Homogenized Korteweg–de Vries and Boussinesq models for anisotropic propagation of solitary waves over a structured bathymetry

Abstract

We derive effective Boussinesq and Korteweg–de Vries equations governing nonlinear wave propagation over a structured bathymetry using a three-scale homogenization approach. The model captures the anisotropic effects induced by the bathymetry, leading to significant modifications in soliton dynamics. Homogenized parameters, determined from elementary cell problems, reveal strong directional dependencies in wave speed and dispersion. Our results provide new insights into nonlinear wave propagation in structured shallow-water environments, and consequently motivate further fundamental and applied studies in wave hydrodynamics and coastal engineering.

1. Introduction

The seminal model proposed by Boussinesq (1871) is a reduced-order model for wave propagation in a shallow-water regime over a flat bottom, characterized by weak nonlinearity and weak dispersion. Over the years, numerous extensions of this model have been developed to account for increasingly complex physical phenomena; comprehensive reviews are provided in Madsen & Schäffer (1999) and Brocchini (2013). A significant body of work has focused on the development of reduced Boussinesq-type models for wave propagation over slowly varying seabed bathymetry, where smooth variations occur on a scale comparable to the wavelength (Madsen & Sørensen 1992; Beji & Nadaoka 1996; Agnon, Madsen & Schäffer 1999; Madsen et al. 2002, 2006; Klopman, Van Groesen & Dingemans 2010). The limitations of this regime become evident when the characteristic length of seabed variations is comparable to the water depth, such as in the case of a sharp step or a rapidly oscillating microstructure. In such cases, it becomes necessary to account for evanescent fields around obstacles and potential velocity singularities in the presence of irregular geometries or sudden changes in depth, see for instance Tuck (2005).

Due to the complexity of the nonlinear water-wave problem, most of the literature in this field focuses on the linear regime, which has seen a resurgence of interest due to recent developments in metamaterials for water-wave control. These metamaterials often involve rapidly varying bathymetries composed of rigid vertical plates placed on the seabed (Porter & Newman 2014). Such structures have been experimentally applied to achieve wave manipulation and cloaking (Farhat et al. 2008; Berraquero et al. 2013; Xu et al. 2015; Li et al. 2018; Hua et al. 2022). Note that, while most studies on water-wave metamaterials focus on linear regimes, a few recent works have explored nonlinear effects, albeit in configurations based on local resonators rather than varying bathymetry (Bobinski et al. 2015; Monsalve et al. 2019; Lorenzo et al. 2023).

Homogenization theories provide a robust framework for understanding the macroscopic behaviour of waves over structured bathymetries. For example, Maurel et al. (2017) and Marangos & Porter (2021) demonstrated that the homogenized model of a rapidly oscillating bathymetry accurately reproduces the anisotropy of experimentally determined dispersion relations. Moreover, Maurel, Pham & Marigo (2019) extended this approach to finite-sized bathymetries by deriving homogenized transmission conditions.

That said, the application of homogenization techniques in the nonlinear regime has received significantly less attention. The pioneering work of Rosales & Papanicolaou (1983) introduced an effective Korteweg–de Vries (KdV)-type equation for waves propagating over a periodic microstructure in a two-dimensional setting, using two-scale homogenization techniques. However, it took over two decades for this result to be extended to a three-dimensional and bidirectional framework by Craig et al. (2005), who derived a homogenized Boussinesq equation for rapidly varying seabeds using Hamiltonian perturbation theory. Note that in this study, the three-dimensional model assumes a scale separation between the two horizontal directions, as it aims to model an effective Kadomtsev–Petviashvili equation. Simultaneously, Nachbin (2003) proposed a theoretical framework based on conformal mapping techniques to derive one-dimensional Boussinesq-type systems with variable coefficients that depend on terrain-following coordinates. This framework was later generalized to two-dimensional propagation by Andrade & Nachbin (2018). Building on these developments, Garnier, Kraenkel & Nachbin (2007) derived a one-dimensional effective Boussinesq equation with explicit expressions for the homogenized coefficients in the case of sinusoidal bathymetry, recovering the results of Rosales & Papanicolaou (1983). We also note the work of Quezada de Luna & Ketcheson (2021) and Ketcheson, Lóczi & Russo (2025), which adopt a two-step asymptotic approach. In these works, the authors start from a reduced Saint-Venant-type equation, where the vertical structure of the flow is incorporated through an effective water depth that depends on the horizontal coordinates. Homogenization is then applied to this reduced equation. This method assumes slowly varying bathymetry and, in the periodic case, that the period is large compared with the water depth. Despite these advances, extending homogenization techniques to fully three-dimensional settings without horizontal scale separation remains an open challenge.

In this paper, we revisit the problem considered in Rosales & Papanicolaou (1983) in a three-dimensional setting, focusing on a rapidly oscillating bathymetry that is invariant in one-direction, see figure 1. To this end, we build on the three-scale homogenization framework introduced in Monsalve, Pham & Maurel (2024) for the analysis of step-like topographies in the weakly dispersive and weakly nonlinear regime. This approach involves three scales: the macroscopic scale of the incident wavelength; the mesoscopic scale associated with the varying bathymetry (the so-called periodic elementary cell); the microscopic scale corresponding to the water surface elevation. These scales are connected through asymptotic matching conditions. A key feature of the method is the interpretation of the surface elevation as a boundary layer phenomenon at the mesoscopic scale as well as the accountability of the singularity of the velocity near sharp edges of the bathymetry. We apply this methodology to the case of a periodic bathymetry, recovering the two-dimensional results of Rosales & Papanicolaou (1983). By generalizing the model to three dimensions, without assuming horizontal scale separation, we provide an explicit prediction of the anisotropic effects induced by irregular seabeds on the propagation of KdV solitons.

Figure 1. Nonlinear wave propagating obliquely over a structured bathymetry.

This article is organized as follows. Section 2 recalls the classical Boussinesq model over a flat bottom and the stationary soliton-type solution of the associated KdV equation. We then introduce the homogenized version of the Boussinesq model for rapidly oscillating bathymetry, involving effective parameters deduced from elementary problems on a unit cell. Due to the anisotropy of the effective medium, the homogenized model leads to an effective KdV equation that depends on the propagation direction of the soliton. Section 2 is devoted to the construction of the homogenized model using a three-scale asymptotic analysis. Finally, in § 3, we analyse the physical properties of the model, with a focus on soliton sensitivity to bathymetric features such as spacing, obstacle height and incidence angle.

2. Anisotropic Boussinesq and KdV models for structured bathymetry

We consider nonlinear propagation of water waves over a structured bathymetry characterized by a periodic alternation of depths between $h$ and $\xi h$ , $0\leqslant \xi \lt 1$ and periodicity $p$ relative to the water depth, see figure 2. Assuming an inviscid, incompressible fluid in irrotational motion, the velocity $\boldsymbol u$ and the associated velocity potential $\varphi$ satisfy

(2.1) \begin{equation} {\textrm {div}}^{\text{3d}}\, {\boldsymbol u}({\boldsymbol r},z,t)=0, \quad {\boldsymbol u}({\boldsymbol r},z,t)={\boldsymbol{\nabla }^{\text{3d}}}\varphi ({\boldsymbol r},z,t), \end{equation}

where $g$ is the gravitational constant, $t$ is time, ${\boldsymbol r}=(x,y)$ represents the coordinate in the horizontal plane and $z$ the vertical coordinate. These equations (2.1) are complemented by the nonlinear dynamic and kinematic boundary conditions at the free surface, i.e.

(2.2) \begin{equation} \displaystyle \frac {\partial \varphi }{\partial t}+\frac {{\boldsymbol u}\boldsymbol{\cdot }{\boldsymbol u}}{2}+g\eta =0,\quad \displaystyle u_z=\frac {\partial \eta }{\partial t}+{\boldsymbol u}\boldsymbol{\cdot }{\boldsymbol{\nabla }} \eta \quad \text{ at } z=\eta ({\boldsymbol r},t), \end{equation}

(we have defined ${\boldsymbol{\nabla }} f=\partial _xf\,\boldsymbol{e}_x+\partial _y\,f\boldsymbol{e}_y$ and ${\boldsymbol{\nabla }^{\text{3d}}} f={\boldsymbol{\nabla }} f+\partial _zf \boldsymbol{e}_z$ ) and condition of vanishing normal velocity

(2.3) \begin{equation} {\boldsymbol u}\boldsymbol{\cdot }{\boldsymbol n}=0 \quad \text{on the rigid walls.} \end{equation}

In the remainder of this section, we first recall the forms of the classical Boussinesq and KdV equations for constant bathymetry. We then present the main contribution of this work, which is the generalization of these equations to the case of a rapidly varying periodic bathymetry.

2.1. The classical Boussinesq and KdV models for a flat bathymetry

For weakly nonlinear and fairly long waves compared with a constant reference water depth $h$ , the classical Boussinesq equations are deduced from the dimensional reduction in the vertical coordinate of (2.1)–(2.3). The resulting, effective, system governing the free surface elevation $\eta ({\boldsymbol r},t)$ and the horizontal velocity $\boldsymbol{u}({\boldsymbol r},t)=u_x({\boldsymbol r},t)\boldsymbol{e}_x+u_y({\boldsymbol r},t)\boldsymbol{e}_y$ at the free surface under such hypothesis reads as

(2.4) \begin{equation} \left \{\begin{array}{l} \displaystyle \displaystyle \frac {\partial \eta }{\partial t} +{\textrm {div}}\left [ ({h}+\eta )\boldsymbol{u} \right ] +\frac {{h}^3}{3}\,{\textrm {div}}\left [ \Delta \boldsymbol{u} \right ]=0,\\[12pt] \displaystyle \frac {\partial \boldsymbol{u}}{\partial t}+g\boldsymbol{\nabla }\eta +\frac {1}{2}\boldsymbol{\nabla }\boldsymbol (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u})=0. \end{array}\right . \end{equation}

Figure 2. (a) Nonlinear wave propagation over a bathymetry that is periodic in $x$ and invariant in $y$ . (b) The effective wave parameters are derived from elementary problems set in the two-dimensional unit cell $\varOmega$ using rescaled coordinates $x_{{m}}=x/{h}$ and ${z_{{m}}}=z/{h}$ , and we set ${S}=|\varOmega |$ the unit cell area.

From the Boussinesq equations, a KdV equation can be obtained (see (3.102) and (3.107) in Whitham (2011)) by considering solitary wave solutions of the form $\eta ({\boldsymbol r},t)=\eta _{\scriptscriptstyle 0} f(x-u t)$ , corresponding to one-directional wave of amplitude $\eta _{\scriptscriptstyle 0}$ propagating in the $\boldsymbol{e}_x$ direction with celerity $c$ , leading to

(2.5) \begin{equation} \frac {\partial \eta }{\partial t}+{c}\left (1+\frac {3}{2}\frac {\eta }{{h}}\right )\frac {\partial \eta }{\partial x}+\frac {{c}{h}^2}{6} \frac {\partial ^3 \eta }{\partial x ^3}=0,\quad {c}=\sqrt {g{h}}, \end{equation}

resulting in the form of a soliton of speed $u$

(2.6) \begin{equation} \eta (x,t)=\eta _{\scriptscriptstyle 0}\text{sech}^2\left \{ \left (\frac {3\eta _{\scriptscriptstyle 0}}{4{h}^3}\right )^{1/2}\left (x-u t\right )\right \}, \quad u={c}\left (1+\frac {\eta _{\scriptscriptstyle 0}}{2{h}}\right ). \end{equation}

As expected for a constant depth, the classical Boussinesq and KdV models are isotropic and thus insensitive to the direction of wave propagation.

2.2. Effective model for a periodic bathymetry

In the case of a fast-varying periodic bathymetry, see figure 2, the dimensional reduction of the governing equations leads to an anisotropic version of the Boussinesq and KdV equations previously given by (2.4) and (2.5) in the case of a constant bathymetry. The effective problem is set on the surface elevation $\overline {\eta }({\boldsymbol r},t)$ and velocity $\overline {{\boldsymbol u}}({\boldsymbol r},t)$ at the free surface, corresponding to averaged versions of the actual field over a periodic cell (hence homogenizing fast oscillations associated with the presence of the periodic microstructure at the bottom).

2.2.1. Anisotropic Boussinesq model

The main result of this study is the derivation of the effective Boussinesq equations for a structured bathymetry, which read as

(2.7) \begin{align} \begin{cases} \displaystyle \displaystyle \frac {\partial \overline {\eta }}{\partial t} +{\textrm {div}}\left [ (\alpha _x\, \overline {h}+n_x \overline {\eta })\overline {u}_x\boldsymbol{e}_x+(\overline {h}+\overline {\eta })\overline {u}_y\boldsymbol{e}_y \right ]\\[10pt] \displaystyle \qquad + \,\overline {h}{h}^2\,{\textrm {div}}\left [ \left (d_{xx}\frac {\partial ^2 \overline {u}_x}{\partial x ^2}+d_{xy}\frac {\partial ^2 \overline {u}_x}{\partial y ^2}\right )\boldsymbol{e}_x + \left (d_{yx}\frac {\partial ^2 \overline {u}_y}{\partial x ^2}+d_{yy}\frac {\partial ^2 \overline {u}_y}{\partial y ^2}\right )\boldsymbol{e}_y \right ]=0,\\[19pt] \displaystyle \frac {\partial \overline {{\boldsymbol u}}}{\partial t}+g\boldsymbol{\nabla } \overline {\eta }+\frac {1}{2}\boldsymbol{\nabla }\left (n_x\overline {u}_x^2+\overline {u}_y^2\right )=0, \end{cases} \end{align}

involving the largest depth $h$ and average water depth $\overline {h}$ , a non-dimensional parameter $\alpha _x$ responsible of the anisotropy in the linear non-dispersive regime, four non-dimensional parameters $(d_{xx},d_{xy},d_{yx},d_{yy})$ accounting for the dispersive effects and one non-dimensional parameter $n_x$ encapsulating the effect of the anisotropic microstructure in the nonlinear regime. These parameters are given by the resolution of linear static problems, of the Poisson-type, which depend only on the geometry of the periodic unit cell, see § 2.2.3. For the case of a flat bathymetry, see tables 1 and 2, these coefficients become explicit and we recover the classical Boussinesq equations given in (2.4). It should be noted that the Boussinesq system in (2.7) corresponds, in the flat-bed limit of our model, to the so-called Kaup system or the integrable version of the Boussinesq system due to Krishnan (1982). This system has been studied in Chen (1998) and Bona, Chen & Saut (2002) as one of the possible forms, obtained via the Boussinesq trick and modulo a choice of the effective velocity reference level (resulting in a family of so-called a–b–c–d systems). For our model, applying the Boussinesq trick poses no particular difficulty. However, in the presence of a variable bathymetry, the choice of the effective velocity level becomes more delicate and requires revisiting the asymptotic derivation. In particular, it would be necessary to perform the asymptotic analysis for different reference levels to establish the equivalent of the anisotropic a–b–c–d systems and then carry out a corresponding stability study. For now, we highlight this as an important direction for future work.

2.2.2. Anisotropic KdV equation and soliton propagation

A key feature of the effective Boussinesq equations (2.7) is their ability to admit anisotropic soliton-type solutions in the homogenized sense, when the influence of the microstructured bathymetry has been spatially averaged. Considering a one-directional effective solution of the form $\overline {\eta }({\boldsymbol r},t)=\eta _{\scriptscriptstyle 0} f(s-u_\theta t)$ , representing a wave of amplitude $\eta _{\scriptscriptstyle 0}$ propagating along the direction $\boldsymbol{e}_\theta$ , which makes an angle $\theta$ with $\boldsymbol{e}_x$ with celerity $u_\theta$ . With $s={\boldsymbol r}\boldsymbol{\cdot }\boldsymbol{e}_\theta$ the local coordinate along the propagation direction, we find the following effective KdV equation:

(2.8) \begin{equation} \frac {\partial \overline {\eta }}{\partial t}+c_\theta \left (1+\frac {3}{2}\frac {\overline {\eta }}{h_\theta }\right )\frac {\partial \overline {\eta }}{\partial s}+\gamma _\theta \frac {c_\theta h_\theta ^2}{6} \frac {\partial ^3 \overline {\eta }}{\partial s ^3}=0,\quad c_\theta =\sqrt {gH_\theta }, \end{equation}

where the anisotropic coefficients $(h_\theta ,H_\theta ,\gamma _\theta )$ are defined as

(2.9) \begin{equation} \left \{\begin{array}{l} \displaystyle h_\theta =\frac {\alpha _x\cos ^2\theta +\sin ^2\theta }{n_x\cos ^2\theta +\sin ^2\theta }\;\overline {h}, \quad H_\theta = \big(\alpha _x\cos ^2\theta +\sin ^2\theta \big)\;\overline {h},\\[15pt] \displaystyle \gamma _\theta =3\frac {\overline {h}}{H_\theta }\left (\frac {{h}}{h_\theta }\right )^2\big(d_{xx}\cos ^4\theta +(d_{xy}+d_{yx})\cos ^2\theta \sin ^2\theta +d_{yy}\sin ^4\theta \big), \end{array}\right . \end{equation}

which yields a $\theta$ -dependent soliton given by

(2.10) \begin{equation} \overline {\eta }(s,t)=\eta _{\scriptscriptstyle 0}\;\text{sech}^2\left \{ \left (\frac {3\eta _{\scriptscriptstyle 0}}{4\gamma _\theta h_\theta ^3}\right )^{1/2}\left (s-u_\theta t\right )\right \}, \quad u_\theta =c_\theta \left (1+\frac {\eta _{\scriptscriptstyle 0}}{2h_\theta }\right ). \end{equation}

A detailed analysis of the property of this $\theta$ -dependent soliton is provided in § 4. Note that in the case of a flat bathymetry (see table 2), one recovers $c_\theta ={c}$ , $h_\theta ={h}$ and $\gamma _\theta =1$ , thus retrieving the classical KdV (2.5) and the isotropic soliton solution (2.6). In the following, we consider the solution (2.10). We emphasize, however, that the stability of the standard KdV equation has been demonstrated by Benjamin (1972) and Bona (1975). With a proper renormalization of (2.8), these stability results directly apply to our anisotropic KdV equation, independently of the propagation direction defined by $\theta$ .

2.2.3. Effective parameters and elementary problems

The non-dimensional parameters $\alpha _x$ , $n_x$ and $d_{ab}$ , $\{a,b\}\in \{x,y\}$ , governing the Boussinesq (2.7), are computed by solving a sequence of linear, static, Poisson-type problems defined over the periodic unit cell $\varOmega$ with local coordinates ${\boldsymbol r}_{{m}}=(x_{{m}},{z_{{m}}})=(x,z)/{h}$ . There are five such elementary problems, denoted $(Q^{\text{(0)}},Q^{\text{(1)}}_x,Q^{\text{(1)}}_y,Q^{\text{(2)}}_x,Q^{\text{(2)}}_y)$ , each satisfying a generic form

(2.11) \begin{equation} \left \{\begin{array}{l} \displaystyle {\textrm {div}}_{\hspace {-.02cm}{m}}\left ({\boldsymbol{\nabla }}_{\hspace {-.05cm}{m}} Q+F\boldsymbol{e}_x\right )=G, \quad \text{in } \varOmega ,\\[12pt] \displaystyle \left ({\boldsymbol{\nabla }}_{\hspace {-.05cm}{m}} Q+F\boldsymbol{e}_x\right )\boldsymbol{\cdot }{\boldsymbol{n}}=0, \quad \text{ on the walls,}\\[10pt] \displaystyle \dfrac {\partial Q}{\partial {z_{{m}}}}(x_{{m}},0)=H, \quad \displaystyle \int _{{z_{{m}}}=0}Q(x_{{m}},0)\,\textrm {d}x_{{m}}=0,\\[15pt] Q,\left ({\boldsymbol{\nabla }}_{\hspace {-.05cm}{m}} Q+F\boldsymbol{e}_x\right ) \text{ $p$-periodic with respect to } x_{{m}}. \end{array}\right . \end{equation}

The source terms $F({\boldsymbol r}_{{m}})$ , $G({\boldsymbol r}_{{m}})$ and $H(x_{{m}})$ vary across the five problems and are given in table 1. Since some elementary problems rely on the solutions of others, they must be solved in a sequential order, starting with $Q^{\text{(0)}}$ , then computing $Q^{\text{(1)}}_x$ and $Q^{\text{(1)}}_y$ , and finally $Q^{\text{(2)}}_x$ and $Q^{\text{(2)}}_y$ .

Table 1. The five elementary problems follow the same Poisson-type cell problem (2.11). For each, we provide the corresponding bulk sources $F({\boldsymbol r}_{{m}})$ , $G({\boldsymbol r}_{{m}})$ and surface source $H(x_{{m}})$ .

Table 2. Explicit solutions for the elementary problems in the flat bathymetry case (with ${S}=p$ and $\overline {h}={h}$ ), and the corresponding effective parameters from (2.12)–(2.13) .

The effective, dimensionless, parameters $\alpha _x$ and $n_x$ are given by

(2.12) \begin{equation} \alpha _x=1+\frac {1}{{S}}\int _\varOmega \frac {\partial Q^{\text{(0)}}}{\partial x_{{m}}}({\boldsymbol r}_{{m}})\,\textrm {d}{\boldsymbol r}_{{m}}, \quad \displaystyle n_x=1+\frac {1}{p}\int _{{z_{{m}}}=0} \left (\frac {\partial Q^{\text{(0)}}}{\partial x_{{m}}}(x_{{m}},0)\right )^2\,\textrm {d}x_{{m}}, \end{equation}

where $S$ denotes the area of the unit cell $\varOmega$ , while the parameters $(d_{xx},d_{xy},d_{yx},d_{yy})$ are defined by

(2.13) \begin{equation} \begin{cases} \displaystyle d_{yx}=\frac {1}{{S}}\int _\varOmega Q^{\text{(1)}}_x({\boldsymbol r}_{{m}})\,\textrm {d}{\boldsymbol r}_{{m}},\quad \displaystyle d_{yy}=\frac {1}{{S}}\int _\varOmega Q^{\text{(1)}}_y({\boldsymbol r}_{{m}})\,\textrm {d}{\boldsymbol r}_{{m}},\\[15pt] \displaystyle d_{xx}=d_{yx}+\frac {1}{{S}}\int _\varOmega \frac {\partial Q^{\text{(2)}}_x}{\partial x_{{m}}}\,\textrm {d}{\boldsymbol r}_{{m}},\quad \displaystyle d_{xy}=d_{yy}+\frac {1}{{S}}\int _\varOmega \frac {\partial Q^{\text{(2)}}_y}{\partial x_{{m}}}\,\textrm {d}{\boldsymbol r}_{{m}}. \end{cases} \end{equation}

For a flat bathymetry, the elementary problems admit explicit solutions, leading to classical values for the effective parameters of the standard Boussinesq model (2.4), as summarized in table 2.

Figure 3. Non-dimensional scales of the problem, from (3.1)–(3.2). (a) The microscopic scale is associated with the wave amplitude $ka=\alpha \delta ^3$ (the macroscopic scale is that of the non-dimensional wavelength $O(1)$ ); (b) the intermediate, mesoscopic, scales are associated with the water depth $kh=\delta$ and array spacing $p\delta =O(\delta )$ . Panels (a ii) and (b ii) show the corresponding systems of coordinates in the microscopic domain and in the mesoscopic domain, within the unit cell, see (3.6).

3. Asymptotic derivation of the homogenized Boussinesq model

3.1. Summary of the asymptotic procedure

To capture in a single asymptotic framework the contributions due to nonlinearities and dispersion, we use classical scalings for characteristic wavelength $1/k$ , water depth $h$ and wave amplitude $a$ . Specifically, by introducing the dimensionless nonlinearity and shallowness parameters

(3.1) \begin{equation} \varepsilon =\frac {{a}}{{h}}\ll 1, \quad \delta =k{h}\ll 1. \end{equation}

Moreover, we consider the following scaling:

(3.2) \begin{equation} \varepsilon =\alpha \delta ^2,\quad \text{with }\,\, \alpha =O(1), \end{equation}

which ensures a non-trivial interplay between nonlinearities and dispersion as well as keeping a single small parameter $\delta$ in the study. In addition, the structured bathymetry has periodicity and step height of the order of $h$ . In the following sections, we will conduct the asymptotic analysis in the different regions shown in figure 3 and, for simplicity, we will consider the non-dimensional form of (2.1)–(2.3). Introducing leading-order dispersion relation $\omega =\sqrt {g{h}}\, k$ , we use the following transformations:

(3.3) \begin{equation} t\to \omega t, \quad (x,y,z)\to k (x,y,z), \quad \eta \to \eta /{a},\quad {\boldsymbol u} \to {\boldsymbol u}/(\omega {a}), \quad \varphi \to k \varphi /(\omega {a}), \end{equation}

which provide the non-dimensional form of (2.1). The incompressibility ( $\text{Inc}$ ), and irrotationality conditions ( $\text{Rot}$ ) become

(3.4) \begin{equation} \begin{array}{l} (\text{Inc}) :\hspace {.2cm} \displaystyle {\textrm {div}}^{\text{3d}}\,{\boldsymbol u}=0,\quad (\text{Rot}):\hspace {.2cm} \displaystyle {\boldsymbol u}={\boldsymbol{\nabla }^{\text{3d}}} \varphi , \end{array} \end{equation}

while the dynamic ( $\text{DC}$ ) and kinematic ( $\text{KC}$ ) conditions at the free surface, described by (2.2) and (2.3), as well as the vanishing normal velocity on the rigid walls formed by the seabed at and the walls ( $\text{RC}$ ) read

(3.5) \begin{align}& (\text{DC}):\hspace {.2cm} \frac {\partial \varphi }{\partial t}+\varepsilon \delta \;\frac {{\boldsymbol u}\boldsymbol{\cdot }{\boldsymbol u}}{2}+\delta ^{-1}\eta =0, \quad (\text{KC}):\hspace {.2cm} u_z=\frac {\partial \eta }{\partial t}+\varepsilon \delta \, {\boldsymbol u}\boldsymbol{\cdot }{\boldsymbol{\nabla }} \eta ,\quad \text{at}\; z=\varepsilon \delta \eta ({\boldsymbol r},t),\nonumber\\[5pt]& (\text{RC}):\hspace {.2cm} {\boldsymbol u}\boldsymbol{\cdot }{\boldsymbol{n}}=0\quad \text{on the walls}. \end{align}

In the above equations, we distinguish three different scales. We first have the macroscopic scale which concerns the horizontal variations that occur at the wavelength scale. Next, we have the mesoscopic scale, which is associated with the scale of the water depth and the bathymetry. Finally, we have the microscopic scale which governs the scale of the wave amplitude $\varepsilon \delta$ in the vicinity of the free surface. Note that the macroscopic and mesoscopic scales are classic for linear shallow-water waves, while the microscale is linked to the (weakly) nonlinear context. This microscopic region corresponds to a zoom revealing the finite wave amplitude, very small compared with the water depth (strictly zero in a linear analysis), and its study will provide the boundary condition for the mesoscopic region. Based on these considerations, we introduce rescaled coordinates capable of taking these different scales into account vertically and horizontally,

(3.6) \begin{equation} \text{macro: } {\boldsymbol r}=(x,y),\quad \text{meso: } {\boldsymbol r}_{{m}}=(x_{{m}},{z_{{m}}})= \left (\frac {x}{\delta },\frac {z}{\delta }\right ), \quad \text{micro: } z_{\mu }=\frac {{z_{{m}}}-\varepsilon \eta }{\varepsilon \eta }. \end{equation}

At the mesoscopic scale in the region of the bathymetry or at the microscopic scale near the free surface, the procedure consists in using the appropriate spatial coordinates (3.6) to get a unitary domain that depends neither on $\delta$ or $\varepsilon$ (it has $O(1)$ dimensions) while the dependence on the small parameter appears explicitly in a new system of equations. The two regions are then connected by asymptotically matching their solutions.

3.2. Setting of the multiscale problem

3.2.1. The mesoscopic domain

The mesoscopic domain $\varOmega$ consists of a two-dimensional periodic cell of vertical extend $\delta$ and horizontal unitary extend. It includes the microstructural bathymetry at the sea bottom but it excludes the free surface of vertical extend $\varepsilon \delta \ll \delta$ , see figure 3. At this scale, the fields follow a two-scale expansion with respect to the macroscopic scale $\boldsymbol r$ and mesoscopic scale ${\boldsymbol r}_{{m}}$ ,

(3.7) \begin{equation} {\boldsymbol u}=\sum _{n\geqslant 0} \delta ^n {\boldsymbol u}^{(n)}({\boldsymbol r},{\boldsymbol r}_{{m}},t),\quad \varphi =\sum _{n\geqslant 0} \delta ^n \varphi ^{(n)}({\boldsymbol r},{\boldsymbol r}_{{m}},t), \end{equation}

with ${\boldsymbol u}^{(n)}$ and $\varphi ^{(n)}$ being $p$ -periodic with respect to $x_{{m}}$ . The differential operator is thus redefined as

(3.8) \begin{equation} {\boldsymbol{\nabla }^{\text{3d}}}\to {\boldsymbol{\nabla }}+\frac {1}{\delta }{\boldsymbol{\nabla }}_{\hspace {-.05cm}{m}}, \quad {\boldsymbol{\nabla }}_{\hspace {-.05cm}{m}} =\boldsymbol{e}_x\frac {\partial }{\partial x_{{m}}}+\boldsymbol{e}_z\frac {\partial }{\partial {z_{{m}}}}, \end{equation}

and it follows that (3.4)–(3.5) in the mesoscopic domain read

(3.9) \begin{align} & {}^{{m}}(\text{Inc})^{}:\hspace {.2cm} {\textrm {div}}{\boldsymbol u}+\frac {1}{\delta }{\textrm {div}}_{\hspace {-.02cm}{m}} {\boldsymbol u}=0,\quad {}^{{m}}(\text{Rot})^{}_{}:\hspace {.2cm} {\boldsymbol u}={\boldsymbol{\nabla }} \varphi +\frac {1}{\delta }{\boldsymbol{\nabla }}_{\hspace {-.05cm}{m}}{\varphi },\nonumber\\[5pt]& {}^{{m}}(\text{RC})^{}:\hspace {.2cm} {\boldsymbol u}\boldsymbol{\cdot }{\boldsymbol {n}}=0\quad \text{on the walls}. \end{align}

We insist on the fact that the mesoscopic region excludes the free surface: the values of the fields at ${z_{{m}}}=0$ are unknowns at this stage and their determination will result from a matched asymptotic analysis with the microscopic domain introduced in the next section.

3.2.2. The microscopic domain

The microscopic domain is the region located at the free surface which results from a zoom made at a point $x_{{m}}$ near the boundary ${z_{{m}}}=0$ in the mesoscopic domain. It corresponds to a water column of infinite vertical extent where nonlinearities are governed by the dynamic/kinematic boundary conditions at the free surface. Here, we have to account for the variations of the fields on the macroscopic scale over $x$ and $y$ , on the mesoscopic scale over $x_{{m}}$ but now over the microscopic scale $z_{\mu }$ . Accordingly, we have a three-scale expansion as follows:

(3.10) \begin{equation} \begin{array}{ll} \displaystyle {\boldsymbol u}=\sum _{n\geqslant 0} \delta ^n {\boldsymbol v}^{(n)}({\boldsymbol r},x_{{m}},z_{\mu },t),\quad \displaystyle \varphi =\sum _{n\geqslant 0} \delta ^n \psi ^{(n)}({\boldsymbol r},x_{{m}},z_{\mu },t), \quad & \displaystyle \eta = \sum _{n\geqslant 1} \delta ^n \eta ^{(n)}({\boldsymbol r},x_{{m}},t). \end{array} \end{equation}

The differential operator reads at this scale

(3.11) \begin{equation} {\boldsymbol{\nabla }^{\text{3d}}} \to {\boldsymbol{\nabla }}+\frac {\boldsymbol{e}_x}{\delta }\frac {\partial }{\partial x_{{m}}}+\frac {\boldsymbol{e}_z}{\alpha \delta ^3\eta }\frac {\partial }{\partial z_{\mu }}. \end{equation}

It follows that (3.4)–(3.5) in the microscopic domain become

(3.12) \begin{align}\begin{split}& {}^{\mu}(\text{Inc}):\hspace {.2cm} {\textrm {div}}{\boldsymbol u}+\frac {1}{\delta }\frac {\partial u_x}{\partial x_{{m}}}+\frac {1}{\alpha \delta ^3\eta }\frac {\partial u_z}{\partial z_{\mu }}=0,\\[5pt]& {}^{\mu}(\text{Rot})^{}_{}:\hspace {.2cm} u_x=\frac {\partial \varphi }{\partial x}+\frac {1}{\delta }\frac {\partial \varphi }{\partial x_{{m}}},\quad u_y=\frac {\partial \varphi }{\partial y},\quad u_z=\frac {1}{\alpha \delta ^3\eta }\frac {\partial \varphi }{\partial z_{\mu }}, \end{split} \end{align}

along with

(3.13) \begin{align}\begin{split}& {}^{\mu}(\text{DC})^{}:\hspace {.2cm} \frac {\partial \varphi }{\partial t}+\delta ^{-1}\eta +\alpha \delta ^3 \;\frac {{\boldsymbol u}\boldsymbol{\cdot }{\boldsymbol u}}{2}=0,\\[5pt]& {}^{\mu}(\text{KC})^{}:\hspace {.2cm} u_z=\frac {\partial \eta }{\partial t}+\alpha \delta ^3 \left ({\boldsymbol u}\boldsymbol{\cdot }{\boldsymbol{\nabla }} \eta +\frac {1}{\delta }u_x\frac {\partial \eta }{\partial x_{{m}}}\right ),\quad \text{at}\; z_{\mu }=0. \end{split} \end{align}

Note that we have anticipated from ${}^{\mu}(\text{DC})^{}$ and ${}^{\mu}(\text{KC})^{}$ that $\varphi =O(1)$ implies $\eta$ and $u_z$ are $O(\delta )$ , hence $v_z^{\text{(0)}}=\eta ^{\text{(0)}}=0$ in (3.10). These equations have to be complemented with appropriate conditions where $z_{\mu }\to -\infty$ which correspond to matching conditions with the mesoscopic region as we will see in the next subsubsection.

3.2.3. Matching the solutions in the mesoregions and microregions

The missing conditions that we have identified in the mesodomain and microdomain are provided simultaneously by matching the solutions in an intermediate region when ${z_{{m}}}\rightarrow 0$ and $z_{\mu }\to -\infty$ , see figure 3. Specifically, for $\varphi$ , we need to have

(3.14) \begin{equation} \varphi ^{\text{(0)}}({\boldsymbol r},{\boldsymbol r}_{{m}},t)+\delta \varphi ^{\text{(1)}}({\boldsymbol r},{\boldsymbol r}_{{m}},t)+\boldsymbol{\cdots} \sim \psi ^{\text{(0)}}({\boldsymbol r},x_{{m}},z_{\mu },t)+ \delta \psi ^{\text{(1)}}({\boldsymbol r},x_{{m}},z_{\mu },t)+\boldsymbol{\cdots} . \end{equation}

The matching procedure is the same for the other fields. Taking into account (3.6) and using the Taylor expansion

(3.15) \begin{equation} f({\boldsymbol r},x_{{m}},{z_{{m}}},t)=f({\boldsymbol r},x_{{m}},0,t)+\alpha \delta ^2 (z_{ \mu }+1)\eta ({\boldsymbol r},x_{{m}},t)\frac {\partial f}{\partial {z_{{m}}}}({\boldsymbol r},x_{{m}},0,t) +\boldsymbol{\cdots} , \end{equation}

we identify in (3.14) the terms of the same power of $\delta$ , which gives

(3.16) \begin{equation} \left \{\begin{array}{l} \displaystyle \varphi ^{(n)}({\boldsymbol r},x_{{m}},0,t)=\lim _{z_{\mu }\to -\infty } \psi ^{(n)}({\boldsymbol r},x_{{m}},z_{\mu },t), \qquad n=0,1,2,\\[12pt] \displaystyle \varphi ^{\text{(3)}}({\boldsymbol r},x_{{m}},0,t)=\lim _{z_{\mu }\to -\infty } \bigg( \psi ^{\text{(3)}}({\boldsymbol r},x_{{m}},z_{\mu },t)\\[9pt] \qquad\qquad\qquad\qquad\qquad\qquad -\alpha (z_{\mu }+1)\eta ^{\text{(1)}}({\boldsymbol r},x_{{m}},t)\dfrac {\partial \varphi ^{\text{(0)}}}{\partial {z_{{m}}}}({\boldsymbol r},x_{{m}},0,t) \bigg). \end{array}\right . \end{equation}

The same relations are obtained by replacing $\varphi ^{(n)}$ by $u_x^{(n)}$ or $u_z^{(n)}$ and $\psi ^{(n)}$ by $v_x^{(n)}$ or $v_z^{(n)}$ .

3.2.4. The macroscopic fields

We aim at obtaining equations on macroscopic fields which depend only on $\boldsymbol r$ and $t$ . The fields are defined by averaging the following fields at the free surface $X_{{m}}$ :

(3.17) \begin{align} \begin{split} \overline {\varphi }^{(n)}({\boldsymbol r},t) & =\frac {1}{p}\int _{X_{{m}}}\varphi ^{(n)}({\boldsymbol r},x_{{m}},0,t)\,\textrm {d}x_{{m}},\hspace {.2cm} \overline {{\boldsymbol u}}^{(n+1)}({\boldsymbol r},t)={\boldsymbol{\nabla }}\overline {\varphi }^{(n)}({\boldsymbol r},t),\\[6pt] \overline {\eta }^{(n)}({\boldsymbol r},t) & =\frac {1}{p}\int _{X_{{m}}}\eta ^{(n)}({\boldsymbol r},x_{{m}},t)\,\textrm {d}x_{{m}}, \end{split} \end{align}

and as an intermediate mean field we shall also use the average horizontal velocity in the unit cell

(3.18) \begin{equation} \displaystyle {\boldsymbol U}^{(n)}({\boldsymbol r},t)=\frac {1}{{S}}\int _\varOmega \Big (u_x^{(n)}({\boldsymbol r},{\boldsymbol r}_{{m}},t)\boldsymbol{e}_x+u_y^{(n)}({\boldsymbol r},{\boldsymbol r}_{{m}},t)\boldsymbol{e}_y\Big)\,\textrm {d}{\boldsymbol r}_{{m}}. \end{equation}

3.3. Resolution at orders zero and one: non-dispersive linear behaviour

Proposition 3.1. The effective wave equation reads at the dominant orders $n=0, 1$ ,

\begin{align*} \kern60pt\left \{\begin{array}{l} \displaystyle\frac {\partial \overline {\varphi }^{(n)}}{\partial t}+\overline {\eta }^{(n+1)}=0, \quad \displaystyle {\mathrm{div}}\, {\boldsymbol U}^{(n)}+\kappa \frac {\partial \overline {\eta }^{(n+1)}}{\partial t}=0,\qquad\qquad \qquad (3.19a) \\[12pt]\displaystyle {\boldsymbol U}^{(n)}=\alpha _x\frac {\partial \overline {\varphi }^{(n)}}{\partial x}\boldsymbol{e}_x+\frac {\partial \overline {\varphi }^{(n)}}{\partial y}\boldsymbol{e}_y,\qquad \qquad \qquad\quad\qquad \qquad\qquad\quad\!\!\! (3.19b) \end{array}\right. \end{align*}

with $\kappa = ({p}/{{S}})$ and $\alpha _x$ given by (2.12).

Remark. At these dominant orders, the effect of the microstructure is already visible since the wave equation is anisotropic, however, the effects of dispersion and nonlinearities are not captured. In the absence of microstructure, i.e. a flat bathymetry, see table 2, we simply obtain the linear isotropic wave equation in shallow water.

3.3.1. Proof of (3.19) for $n=0$

From ${}^{{m}}(\text{Rot})^{-1}_{}$ in (3.9), we have ${\boldsymbol{\nabla }}_{\hspace {-.05cm}{m}}\varphi ^{\text{(0)}}=0$ , and from ${}^{\mu}(\text{Rot})^{-1}_{x}$ and ${}^{\mu}(\text{Rot})^{-4}_{z}$ in (3.12), $\partial _{x_{{m}}}\psi ^{\text{(0)}}=\partial _{z_{ \mu }}\psi ^{\text{(0)}}=0$ , hence $\psi ^{\text{(0)}}=\psi ^{\text{(0)}}({\boldsymbol r},t)$ . Owing to the matching (3.16), we deduce

(3.20) \begin{equation} \psi ^{\text{(0)}}({\boldsymbol r},t)= \varphi ^{\text{(0)}}({\boldsymbol r},t). \end{equation}

Next, from ${}^{\mu}(\text{Inc})^{-3}$ in (3.12), $\partial _{z_{\mu }}v_z^{\text{(1)}}=0$ . Accounting for ${}^{\mu}(\text{KC})^{1}$ and ${}^{\mu}(\text{DC})^{0}$ as well as the matching on $u_z^{\text{(1)}}$ from (3.16), we obtain

(3.21) \begin{equation} v_z^{\text{(1)}}({\boldsymbol r},t)=-\frac {\partial ^2 \psi ^{\text{(0)}}}{\partial t ^2}({\boldsymbol r},t)=\frac {\partial \eta ^{\text{(1)}}}{\partial t}({\boldsymbol r},t)=u_z^{\text{(1)}}({\boldsymbol r},x_{{m}},{z_{{m}}}=0,t), \end{equation}

which tells us that $u_z^{\text{(1)}}$ does not depend on $x_{{m}}$ . By integrating (3.21) over $X_{{m}}$ together with (3.20), we obtain the first relation in (3.19a ) for $n=0$ . Now, we use ${}^{{m}}(\text{Inc})^{0}$ in (3.9), namely ${\textrm {div}}{\boldsymbol u}^{\text{(0)}}+{\textrm {div}}_{\hspace {-.02cm}{m}}{\boldsymbol u}^{\text{(1)}}=0$ , that we integrate over $\varOmega$ . By doing so and accounting for the boundary condition ${}^{{m}}(\text{RC})^{1}$ in (3.9) on the walls and for the $p$ -periodicity of ${\boldsymbol u}^{\text{(1)}}$ with respect to $x_{{m}}$ , we obtain

(3.22) \begin{equation} {S}{\textrm {div}}{\boldsymbol U}^{\text{(0)}}+\int _{X_{{m}}}u_z^{\text{(1)}}({\boldsymbol r},x_{{m}},{z_{{m}}}=0,t)\,\textrm {d}x_{{m}}=0, \end{equation}

where the second term is deduced from (3.21) and we obtain the second relation in (3.19a ) for $n=0$ .

The relation (3.19b ) is obtained by considering the problem satisfied by $({\boldsymbol u}^{\text{(0)}},\varphi ^{\text{(1)}})$ in the ${\boldsymbol r}_{{m}}\in \varOmega$ coordinate. From ${}^{{m}}(\text{Inc})^{-1}$ , ${}^{{m}}(\text{Rot})^{0}_{}$ , ${}^{{m}}(\text{RC})^{0}$ in (3.9) along with the matching condition (3.16) which provides $u_z^{\text{(0)}}({\boldsymbol r},x_{{m}},{z_{{m}}}=0)=0$ (since $v_z^{\text{(0)}}=0$ ), the problem reads

(3.23) \begin{equation} \left \{\begin{array}{l} \displaystyle {\textrm {div}}_{\hspace {-.02cm}{m}}{\boldsymbol u}^{\text{(0)}}=0,\quad {\boldsymbol u}^{\text{(0)}}={\boldsymbol{\nabla }}_{\hspace {-.05cm}\textit{m}}\varphi ^{\text{(1)}}+{\boldsymbol{\nabla }}\overline {\varphi }^{\text{(0)}}({\boldsymbol r},t),\quad \text{in } \varOmega ,\\[6pt] {\boldsymbol u}^{\text{(0)}}\boldsymbol{\cdot }{\boldsymbol{n}}=0 \text{ on the walls,}\quad u_z^{\text{(0)}}=0 \text{ at } {z_{{m}}}=0,\\[6pt] {\boldsymbol u}^{\text{(0)}},\varphi ^{\text{(1)}} \text{ $p$-periodic with respect to } x_{{m}}. \end{array}\right . \end{equation}

The above problem is linear with respect to the two components of ${\boldsymbol{\nabla }}\overline {\varphi }^{\text{(0)}}$ and noticing that ${\boldsymbol r}_{{m}}\boldsymbol{\cdot }\boldsymbol{e}_y=0$ , the solution can be written as

(3.24) \begin{equation} \varphi ^{\text{(1)}}=\overline {\varphi }^{\text{(1)}}+\frac {\partial \overline {\varphi }^{\text{(0)}}}{\partial x}Q^{\text{(0)}}, \end{equation}

with $\overline {\varphi }^{\text{(1)}}$ defined in (3.17) and $Q^{\text{(0)}}({\boldsymbol r}_{{m}})$ solution to the elementary problem announced in (2.11) and table 1.

Integrating ${\boldsymbol u}^{\text{(0)}}={\boldsymbol{\nabla }}_{\hspace {-.05cm}\textit{m}}\varphi ^{\text{(1)}}+{\boldsymbol{\nabla }}\overline {\varphi }^{\text{(0)}}$ (with (3.24)) over $\varOmega$ , see (3.17), we obtain (3.19b ) with $\alpha _x$ given by (2.12).

3.3.2. Proof of (3.19) for $n=1$

We proceed as for $n=0$ . From ${}^{\mu}(\text{Rot})^{-3}_{z}$ in (3.12), $\partial _{z_{\mu }}\psi ^{\text{(1)}}=0$ , and from ${}^{\mu}(\text{Rot})^{0}_{x,y}$ in (3.12), ${}^{{m}}(\text{Rot})^{0}_{x,y}$ in (3.9) and (3.20) (along with the matchings (3.16)), we have

(3.25) \begin{equation} \psi ^{\text{(1)}}({\boldsymbol r},x_{{m}},t)=\varphi ^{\text{(1)}}({\boldsymbol r},x_{{m}},0,t), \quad v_a^{\text{(0)}}({\boldsymbol r},x_{{m}},t)=u_a^{\text{(0)}}({\boldsymbol r},x_{{m}},0,t), \; a=x,y. \end{equation}

Next, from ${}^{\mu}(\text{Inc})^{-2}$ in (3.12), $\partial _{z_{\mu }}v_z^{\text{(2)}}=0$ hence from ${}^{\mu}(\text{KC})^{2}$ and ${}^{\mu}(\text{DC})^{1}$ , along with the matching $u_z^{\text{(2)}}$ in (3.16), we obtain

(3.26) \begin{equation} v_z^{\text{(2)}}({\boldsymbol r},x_{{m}},t)=-\frac {\partial ^2 \varphi ^{\text{(1)}}}{\partial t ^2}({\boldsymbol r},x_{{m}},0,t)=\frac {\partial \eta ^{\text{(2)}}}{\partial t}({\boldsymbol r},x_{{m}},t)=u_z^{\text{(2)}}({\boldsymbol r},x_{{m}},0,t), \end{equation}

whose integration over $X_{{m}}$ provides the first relation in (3.19a ) for $n=1$ . Now, we use ${}^{{m}}(\text{Inc})^{1}$ in (3.9), namely ${\textrm {div}}{\boldsymbol u}^{\text{(1)}}+{\textrm {div}}_{\hspace {-.02cm}{m}}{\boldsymbol u}^{\text{(2)}}=0$ , that we integrate over $\varOmega$ . Proceeding as previously, we obtain

(3.27) \begin{equation} {S}{\textrm {div}}{\boldsymbol U}^{\text{(1)}}+\int _{X_{{m}}}u_z^{\text{(2)}}({\boldsymbol r},x_{{ m}},0,t)\,\textrm {d}x_{{m}}=0, \end{equation}

and using (3.26), we obtain the second relation in (3.19a ) for $n=1$ .

We now consider the problem satisfied by $({\boldsymbol u}^{\text{(1)}},\varphi ^{\text{(2)}})$ in ${\boldsymbol r}_{{m}}\in \varOmega$ coordinate. From ${}^{{m}}(\text{Inc})^{0}$ , ${}^{{m}}(\text{Rot})^{1}_{}$ , ${}^{{m}}(\text{RC})^{1}$ in (3.9) and (3.21), the problem reads

(3.28) \begin{equation} \left \{\begin{array}{l} \displaystyle {\textrm {div}}_{\hspace {-.02cm}{m}}{\boldsymbol u}^{\text{(1)}}=-{\textrm {div}}{\boldsymbol u}^{\text{(0)}},\quad {\boldsymbol u}^{\text{(1)}}={\boldsymbol{\nabla }}_{\hspace {-.05cm}{m}}\varphi ^{\text{(2)}}+{\boldsymbol{\nabla }}\varphi ^{\text{(1)}},\quad \text{in } \varOmega ,\\[10pt] \displaystyle {\boldsymbol u}^{\text{(1)}}\boldsymbol{\cdot }{\boldsymbol{n}}=0 \text{ on the walls,}\quad u_z^{\text{(1)}}=-\frac {\partial ^2 \overline {\varphi }^{\text{(0)}}}{\partial t ^2} \text{ at } {z_{{m}}}=0,\\[10pt] {\boldsymbol u}^{\text{(1)}},\varphi ^{\text{(2)}} p\text{-periodic with respect to } x_{{m}}. \end{array}\right . \end{equation}

We know $({\boldsymbol u}^{\text{(0)}},\varphi ^{\text{(1)}})$ from (3.24) and we use that $\partial _{tt}\overline {\varphi }^{\text{(0)}}=1/\kappa (\alpha _x\partial _{xx}\overline {\varphi }^{\text{(0)}}+\partial _{yy}\overline {\varphi }^{\text{(0)}})$ from (3.19a )–(3.19b ). The solution at this order can be express as a linear combination of elementary solutions pondered by macroscopic sources. Accordingly, we set

(3.29) \begin{equation} \varphi ^{\text{(2)}}= \overline {\varphi }^{\text{(2)}}+\frac {\partial \overline {\varphi }^{\text{(1)}}}{\partial x}Q^{\text{(0)}}+\frac {\partial ^2 \overline {\varphi }^{\text{(0)}}}{\partial x ^2}Q^{\text{(1)}}_x+\frac {\partial ^2 \overline {\varphi }^{\text{(0)}}}{\partial y ^2}Q^{\text{(1)}}_y, \end{equation}

(and ${\boldsymbol u}^{\text{(1)}}$ given in Appendix A.1) with $\overline {\varphi }^{\text{(2)}}$ defined in (3.17), $Q^{\text{(0)}}({\boldsymbol r}_{{m}})$ , $Q^{\text{(1)}}_{x}({\boldsymbol r}_{{m}})$ and $Q^{\text{(1)}}_{y}({\boldsymbol r}_{{m}})$ solutions to the elementary problems given in (2.11) and table 1.

We notice that $Q^{\text{(1)}}_{x}$ and $Q^{\text{(1)}}_{y}$ are even functions of $x_{{m}}$ (with $Q^{\text{(0)}}$ an odd function of $x_{{m}}$ ) due to the symmetry $x_{{m}}\mapsto -x_{{m}}$ of the obstacle. Accordingly, integrating ${\boldsymbol u}^{\text{(1)}}={\boldsymbol{\nabla }}_{\hspace {-.05cm}\textit{m}}\varphi ^{\text{(2)}}+{\boldsymbol{\nabla }}\varphi ^{\text{(1)}}$ over $\varOmega$ (with (3.24) and (3.29), see also Appendix A.1 for the explicit expression of ${\boldsymbol u}^{\text{(1)}}$ ), we obtain (3.19b ) at order $n=1$ .

3.4. Resolution at order two: dispersive contributions

Proposition 3.2. The effective wave equation reads at the order $n=2$ ,

\begin{align*} \kern70pt\left \{\begin{array}{ll} \displaystyle\frac {\partial \overline {\varphi }^{\text{(2)}}}{\partial t}+\overline {\eta }^{\text{(3)}}=0,\quad \displaystyle {\mathrm{div}}\, {\boldsymbol U}^{\text{(2)}}+\kappa \frac {\partial \overline {\eta }^{\text{(3)}}}{\partial t}=0,\qquad \qquad\qquad\qquad\!\!\! (3.30a)\\[16pt] \displaystyle U_x^{\text{(2)}}=\alpha _x\frac {\partial \overline {\varphi }^{\text{(2)}}}{\partial x}+d_{xx}\frac {\partial ^3 \overline {\varphi }^{\text{(0)}}}{\partial x ^3} +d_{xy}\frac {\partial ^3 \overline {\varphi }^{\text{(0)}}}{\partial x \partial y ^2},\qquad \qquad\qquad\qquad\!\!\! (3.30b)\\[18pt] \displaystyle U_y^{\text{(2)}}=\frac {\partial \overline {\varphi }^{\text{(2)}}}{\partial y}+d_{yx}\frac {\partial ^3 \overline {\varphi }^{\text{(0)}}}{\partial x ^2\partial y} +d_{yy}\frac {\partial ^3 \overline {\varphi }^{\text{(0)}}}{\partial y ^3},\qquad \qquad\qquad\qquad\ \ (3.30c)\\ \end{array}\right . \end{align*}

where $(d_{xx},d_{xy},d_{yx},d_{yy})$ are defined in (2.13).

Remark. Anisotropic dispersive terms appear at this order while nonlinear contributions remains absent. In the absence of microstructure, i.e. a flat bathymetry, see table 2, all the dispersive coefficients $d_{ij}$ are equal to $1/3$ .

3.4.1. Proof of (3.30)

We start with the equivalent at order two of relations (3.25) and (3.26). First, from ${}^{\mu}(\text{Rot})^{-2}_{z}$ in (3.12), we have $\partial _{z_{\mu }}\psi ^{\text{(2)}}=0$ , and from ${}^{\mu}(\text{Rot})^{1}_{x,y}$ in (3.12), ${}^{{m}}(\text{Rot})^{1}_{x,y}$ in (3.9) and (3.25) (along with the matchings (3.16)), we have

(3.31) \begin{equation} \psi ^{\text{(2)}}({\boldsymbol r},x_{{m}},t)=\varphi ^{\text{(2)}}({\boldsymbol r},x_{{m}},0,t), \quad v_a^{\text{(1)}}({\boldsymbol r},x_{{m}},t)=u_a^{\text{(1)}}({\boldsymbol r},x_{{m}},0,t), \; a=x,y. \end{equation}

Next, from ${}^{\mu}(\text{Inc})^{-1}$ in (3.12) and (3.25), we have $\partial _{z_{ \mu }}v_z^{\text{(3)}}=-\alpha \eta ^{\text{(1)}}\partial _{x_{{m}}}{u_x^{\text{(0)}}}_{|{z_{{m}}}=0}$ . We obtain $v_z^{\text{(3)}}$ by integration and using that $v_z^{\text{(3)}}=\partial _t\eta ^{\text{(3)}}$ at $z_{\mu }=0$ from ${}^{\mu}(\text{KC})^{3}$ in (3.13) along with (3.21). With $\eta ^{\text{(3)}}$ satisfying ${}^{\mu}(\text{DC})^{2}$ , we eventually obtain

(3.32) \begin{align} \begin{split}v_z^{\text{(3)}}({\boldsymbol r},x_{{m}},z_{\mu },t)& =\frac {\partial \eta ^{\text{(3)}}}{\partial t}({\boldsymbol r},x_{{m}},t)-z_{\mu }\,\alpha \frac {\partial }{\partial x_{{m}}}\big(\eta ^{\text{(1)}}({\boldsymbol r},t)u_x^{\text{(0)}}({\boldsymbol r},x_{{m}},0,t)\big), \\[5pt] \eta ^{\text{(3)}}({\boldsymbol r},x_{{m}},t)& + \frac {\partial \varphi ^{\text{(2)}}}{\partial t}({\boldsymbol r},x_{{m}},0,t)=0, \end{split} \end{align}

and taking the average over $X_{{m}}$ of the second relation above provides the first relation in (3.30a ). Using further the matching (3.16) on $u_z^{\text{(3)}}$ , we also obtain

(3.33) \begin{equation} u_z^{\text{(3)}}({\boldsymbol r},x_{{m}},0,t)=\frac {\partial \eta ^{\text{(3)}}}{\partial t}({\boldsymbol r},x_{{m}},t)+\alpha \frac {\partial }{\partial x_{{m}}}\left (\eta ^{\text{(1)}}({\boldsymbol r},t)u_x^{\text{(0)}}({\boldsymbol r},x_{{m}},0,t)\right ). \end{equation}

Now, we use ${}^{{m}}(\text{Inc})^{2}$ in (3.9), namely ${\textrm {div}}{\boldsymbol u}^{\text{(2)}}+{\textrm {div}}_{\hspace {-.02cm}{m}}{\boldsymbol u}^{\text{(3)}}=0$ , that we integrate over $\varOmega$ . By doing so and accounting for the boundary condition ${}^{{m}}(\text{RC})^{3}$ in (3.9), we obtain

(3.34) \begin{equation} {S}{\textrm {div}}{\boldsymbol U}^{\text{(2)}}+\int _{X_{{m}}}u_z^{\text{(3)}}({\boldsymbol r},x_{{m}},0,t)\,\textrm {d}x_{{m}}=0. \end{equation}

Reporting (3.33) in the last integral and using the $p$ -periodicity of $\eta ^{\text{(1)}}u_x^{\text{(0)}}$ with respect to $x_{{m}}$ , we obtain the second relation in (3.30a ).

We consider the problem satisfied by $({\boldsymbol u}^{\text{(2)}},\varphi ^{\text{(3)}})$ in ${\boldsymbol r}_{{m}}\in \varOmega$ coordinate. From ${}^{{m}}(\text{Inc})^{1}$ , ${}^{{m}}(\text{Rot})^{2}_{}$ , ${}^{{m}}(\text{RC})^{2}$ in (3.9) and with $u^{\text{(2)}}_z$ at ${z_{{m}}}=0$ given by (3.26), the problem reads

(3.35) \begin{equation} \left \{\begin{array}{l} \displaystyle {\textrm {div}}_{\hspace {-.02cm}{m}}{\boldsymbol u}^{\text{(2)}}=-{\textrm {div}}{\boldsymbol u}^{\text{(1)}},\quad {\boldsymbol u}^{\text{(2)}}={\boldsymbol{\nabla }}_{\hspace {-.05cm}\textit{m}}\varphi ^{\text{(3)}}+{\boldsymbol{\nabla }}\varphi ^{\text{(2)}},\quad \text{in } \varOmega ,\\ \displaystyle {\boldsymbol u}^{\text{(2)}}\boldsymbol{\cdot }{\boldsymbol{n}}=0 \text{ on the walls,}\quad u_z^{\text{(2)}}=-\frac {\partial ^2 \varphi ^{\text{(1)}}}{\partial t ^2} \text{ at } {z_{{m}}}=0,\\ {\boldsymbol u}^{\text{(2)}},\varphi ^{\text{(3)}} \text{ $p$-periodic with respect to } x_{{m}}. \end{array}\right . \end{equation}

We know $({\boldsymbol u}^{\text{(1)}},\varphi ^{\text{(2)}})$ from (3.28) and (3.29), and $\varphi ^{\text{(1)}}$ from (3.24), and this tells us that the problem is linear with respect to the two components of ${\boldsymbol{\nabla }}\overline {\varphi }^{\text{(2)}}$ , $\partial _{xx}\overline {\varphi }^{\text{(1)}}$ and $\partial _{yy}\overline {\varphi }^{\text{(1)}}$ (as at order one) as well as $\partial _{xxx}\overline {\varphi }^{\text{(0)}}$ and $\partial _{xyy}\overline {\varphi }^{\text{(0)}}$ (see Appendix A.2 for the detailed derivation). Accordingly, we set

(3.36) \begin{align} \displaystyle \varphi ^{\text{(3)}}= \displaystyle \overline {\varphi }^{\text{(3)}}+\frac {\partial \overline {\varphi }^{\text{(2)}}}{\partial x}Q^{\text{(0)}}+\frac {\partial ^2 \overline {\varphi }^{\text{(1)}}}{\partial x ^2}Q^{\text{(1)}}_x+\frac {\partial ^2 \overline {\varphi }^{\text{(1)}}}{\partial y ^2}Q^{\text{(1)}}_y+\frac {\partial ^3 \overline {\varphi }^{\text{(0)}}}{\partial x ^3}Q^{\text{(2)}}_x+\frac {\partial ^3 \overline {\varphi }^{\text{(0)}}}{\partial x \partial y ^2}Q^{\text{(2)}}_y, \end{align}

with ${\boldsymbol u}^{\text{(2)}}$ (given in Appendix A.1), $\overline {\varphi }^{\text{(3)}}$ defined in (3.17) with the two additional functions $Q^{\text{(2)}}_{x}$ and $Q^{\text{(2)}}_{y}$ appearing at this order, solution of the static elementary problems defined in (2.11) and table 1.

We note that $Q^{\text{(2)}}_{x}$ and $Q^{\text{(2)}}_{y}$ are odd functions of $x_{{m}}$ (with $Q^{\text{(1)}}_{x,y}$ and $Q^{\text{(0)}}$ being even and odd functions of $x_{{m}}$ , respectively). Accordingly, integrating ${\boldsymbol u}^{\text{(2)}}={\boldsymbol{\nabla }}_{\hspace {-.05cm}\textit{ m}}\varphi ^{\text{(3)}}+{\boldsymbol{\nabla }}\varphi ^{\text{(2)}}$ in (3.35) over $\varOmega$ with (3.29) and (3.36) (see also Appendix A.1 for the full expression of ${\boldsymbol u}^{\text{(2)}}$ ), we obtain the form of ${\boldsymbol U}^{\text{(2)}}$ announced in (3.30b ). Eventually, in the absence of microstructure, with $Q^{\text{(0)}}=0$ , $Q^{\text{(1)}}_{x}$ and $Q^{\text{(1)}}_{y}$ independent of $x_{{m}}$ , we recover that $Q^{\text{(2)}}_{x}=Q^{\text{(2)}}_{y}=0$ , as announced in table 2.

3.5. Resolution at order three: nonlinear contributions

Proposition 3.3. The effective wave equation reads at order $n=3$ ,

\begin{align*} \begin{cases} \displaystyle\frac {\partial \overline {\varphi }^{\text{(3)}}}{\partial t}+\overline {\eta }^{\text{(4)}}+\frac {\alpha }{2}\left ( n_x\left (\frac {\partial \overline {\varphi }^{\text{(0)}}}{\partial x}\right )^2+\left (\frac {\partial \overline {\varphi }^{\text{(0)}}}{\partial y}\right )^2 \right )=0,\\[10pt] \displaystyle {\mathrm{div}}\, {\boldsymbol U}^{\text{(3)}}+\kappa \left (\frac {\partial \overline {\eta }^{\text{(4)}}}{\partial t}+\alpha {\textrm {div}}(\overline {\eta }^{\text{(1)}}{\boldsymbol{\nabla }}\overline {\varphi }^{\text{(0)}})\right )=0, \qquad \qquad \qquad \qquad\qquad\qquad (3.37a) \\[15pt] \displaystyle U_x^{\text{(3)}}=\alpha _x\frac {\partial \overline {\varphi }^{\text{(3)}}}{\partial x}+d_{xx}\frac {\partial ^3 \overline {\varphi }^{\text{(1)}}}{\partial x ^3} +d_{xy}\frac {\partial ^3 \overline {\varphi }^{\text{(1)}}}{\partial x \partial y ^2} -\alpha \kappa (1-n_x) \,\frac {\partial \overline {\varphi }^{\text{(0)}}}{\partial x}{\overline {\eta }^{\text{(1)}}},\qquad \quad (3.37b)\\[14pt] \displaystyle U_y^{\text{(3)}}=\frac {\partial \overline {\varphi }^{\text{(3)}}}{\partial y}+d_{yx}\frac {\partial ^3 \overline {\varphi }^{\text{(1)}}}{\partial x ^2\partial y} +d_{yy}\frac {\partial ^3 \overline {\varphi }^{\text{(1)}}}{\partial y ^3} , \qquad \qquad \qquad \qquad\qquad\qquad\qquad (3.37c)\end{cases} \end{align*}

where $n_x$ and $\alpha _x$ are given by (2.12), and $(d_{xx},d_{xy},d_{yx},d_{yy})$ are defined in (2.13).

Remark. A nonlinear anisotropic contribution appears in (3.37a ) at this order. Such contribution becomes isotropic in the absence of microstructure (with $\kappa =n_x=1$ ). Note that the nonlinear contribution in (3.37b ) vanishes also for $n_x=1$ in the absence of a microstructure.

3.5.1. Proof of (3.37)

From ${}^{\mu}(\text{Rot})^{-1}_{z}$ in (3.12), $\partial _{z_{\mu }}\psi ^{\text{(3)}}=0$ , hence from the matching (3.16) along with $\partial _{{z_{{m}}}}\varphi ^{\text{(0)}}=0$ from (3.20), we deduce

(3.38) \begin{equation} \psi ^{\text{(3)}}({\boldsymbol r},x_{{m}},t)=\varphi ^{\text{(3)}}({\boldsymbol r},x_{{m}},0,t). \end{equation}

Using the result above along with ${}^{\mu}(\text{DC})^{3}$ in (3.13) and (3.25), we have

(3.39) \begin{equation} \frac {\partial \varphi ^{\text{(3)}}}{\partial t}({\boldsymbol r},x_{{m}},0,t)+\eta ^{\text{(4)}}({\boldsymbol r},x_{{m}},t)+\frac {\alpha }{2}\Big((u_x^{\text{(0)}})^2+(u_y^{\text{(0)}})^2\Big)_{|{z_{{m}}}=0}=0, \end{equation}

with $(u_x^{\text{(0)}},u_y^{\text{(0)}})$ given in (3.23) together with (3.24). Integrating (3.39) over $X_{{m}}$ and accounting for the periodicity of $Q^{\text{(0)}}$ with respect to $x_{{m}}$ provides the first relation in (3.37a ).

We now consider $v_z^{\text{(4)}}$ which is given by ${}^{\mu}(\text{Inc})^{0}$ in (3.12). Using (3.25) and (3.31) along with ${}^{{m}}(\text{Inc})^{-1,0}$ in (3.9), we obtain $\partial _{z_{\mu }}{v_z^{\text{(4)}}}=\alpha \partial _{{z_{{m}}}}(\eta ^{\text{(1)}}u_z^{\text{(1)}} +\eta ^{\text{(2)}}u_z^{\text{(0)}})_{|{z_{{m}}}=0}$ . By integration, and using ${}^{\mu}(\text{KC})^{4}$ in (3.13) with $\eta ^{\text{(1)}}$ independent of $x_{{m}}$ from (3.21), we obtain

(3.40) \begin{equation} \begin{array}{l} v_z^{\text{(4)}}=\displaystyle z_{\mu }\, \alpha \frac {\partial }{\partial {z_{{m}}}}(\eta ^{\text{(1)}}u_z^{\text{(1)}} +\eta ^{\text{(2)}}u_z^{\text{(0)}})_{|{z_{{m}}}=0} +{v_z^{\text{(4)}}}_{|z_{\mu }=0},\\[10pt] \displaystyle {v_z^{\text{(4)}}}_{|z_{\mu }=0}=\frac {\partial \eta ^{\text{(4)}}}{\partial t}+\alpha \left (u_x^{\text{(0)}}\frac {\partial \eta ^{\text{(2)}}}{\partial x_{{m}}}+ {\boldsymbol u}^{\text{(0)}}\boldsymbol{\cdot }{\boldsymbol{\nabla }}\eta ^{\text{(1)}}\right )_{|{z_{{m}}}=0}. \end{array} \end{equation}

Now, we use ${}^{{m}}(\text{Inc})^{3}$ in (3.9), namely ${\textrm {div}}{\boldsymbol u}^{\text{(3)}}+{\textrm {div}}_{\hspace {-.02cm}{m}}{\boldsymbol u}^{\text{(4)}}=0$ , that we integrate over $\varOmega$ . By doing so and accounting for the boundary condition ${}^{{m}}(\text{RC})^{4}$ in (3.9), we obtain

(3.41) \begin{equation} {S}{\textrm {div}}{\boldsymbol U}^{\text{(3)}}+\int _{X_{{m}}}u_z^{\text{(4)}}({\boldsymbol r},x_{{m}},0,t)\,\textrm {d}x_{{m}}=0. \end{equation}

To deduce ${u_z^{\text{(4)}}}_{|{z_{{m}}}=0}$ , we extend the matching condition (3.16) at the fourth order, namely

(3.42) \begin{equation} {u_z^{\text{(4)}}}_{|{z_{{m}}}=0}=\lim _{z_{\mu }\to -\infty } \left ( v_z^{\text{(4)}}-\alpha (z_{\mu }+1) \frac {\partial }{\partial {z_{{m}}}}\big(\eta ^{\text{(1)}}u_z^{\text{(1)}}+\eta ^{\text{(2)}}u_z^{\text{(0)}} \big)_{|{z_{{m}}}=0} \right). \end{equation}

Using (3.40), and again ${}^{{m}}(\text{Inc})^{0,1}$ , we obtain

(3.43) \begin{equation} {u_z^{\text{(4)}}}_{|{z_{{m}}}=0}= \frac {\partial \eta ^{\text{(4)}}}{\partial t}+ \alpha {\textrm {div}}(\eta ^{\text{(1)}}{\boldsymbol u}^{\text{(0)}})_{|{z_{{m}}}=0} +\alpha \frac {\partial }{\partial x_{{m}}} \big(\eta ^{\text{(2)}} {u_x^{\text{(0)}}}+\eta ^{\text{(1)}} {u_x^{\text{(1)}}}\big)_{|{z_{{m}}}=0} \end{equation}

(we also used the property that $\eta ^{\text{(1)}}$ does not depend on $x_{{m}}$ from (3.21)). Inserting the above result in (3.41), the integral of the third term of (3.43) vanishes by periodicity, as previously, and with $\eta ^{\text{(1)}}=\overline {\eta }^{\text{(1)}}$ from (3.21), we obtain the second relation in (3.37a ).

We now consider the problem satisfied by $({\boldsymbol u}^{\text{(3)}},\varphi ^{\text{(4)}})$ in ${\boldsymbol r}_{{m}}\in \varOmega$ coordinate. From ${}^{{m}}(\text{Inc})^{2}$ , ${}^{{m}}(\text{Rot})^{3}_{}$ , ${}^{{m}}(\text{RC})^{3}$ in (3.9) and with $u^{\text{(3)}}_z$ at ${z_{{m}}}=0$ given by (3.32), the problem reads

(3.44) \begin{equation} \left \{\begin{array}{l} \displaystyle {\textrm {div}}_{\hspace {-.02cm}{m}}{\boldsymbol u}^{\text{(3)}}=-{\textrm {div}}{\boldsymbol u}^{\text{(2)}},\quad {\boldsymbol u}^{\text{(3)}}={\boldsymbol{\nabla }}_{\hspace {-.05cm}\textit{m}}\varphi ^{\text{(4)}}+{\boldsymbol{\nabla }}\varphi ^{\text{(3)}},\quad \text{in } \varOmega ,\\[6pt]\displaystyle {\boldsymbol u}^{\text{(3)}}\boldsymbol{\cdot }{\boldsymbol{n}}=0 \text{ on the walls,}\quad \\[6pt] \displaystyle u_z^{\text{(3)}}=\frac {\partial \eta ^{\text{(3)}}}{\partial t}({\boldsymbol r},x_{{m}},t)+\alpha \frac {\partial }{\partial x_{{m}}}\left (\eta ^{\text{(1)}}({\boldsymbol r},t)u_x^{\text{(0)}}({\boldsymbol r},x_{{m}},0,t)\right ) \text{ at } {z_{{m}}}=0,\\[6pt]{\boldsymbol u}^{\text{(3)}},\varphi ^{\text{(4)}} \text{ $p$-periodic with respect to } x_{{m}}. \end{array}\right . \end{equation}

We know the expressions for $({\boldsymbol u}^{\text{(2)}},\varphi ^{\text{(3)}})$ from (3.35) and (3.36), $\eta ^{\text{(3)}}$ from (3.32) and (3.29), $\eta ^{\text{(1)}}$ from (3.20) and (3.21) and $u_x^{\text{(0)}}$ from (3.23). This static problem can be decomposed as a linear combination of the macroscopic components $\partial _{x}\overline {\varphi }^{\text{(3)}}$ , $(\partial _{xx}\overline {\varphi }^{\text{(2)}},\partial _{yy}\overline {\varphi }^{\text{(2)}})$ and $(\partial _{xxx}\overline {\varphi }^{\text{(1)}},\partial _{xyy}\overline {\varphi }^{\text{(1)}})$ in a similar way tothe previous order, see (3.36). The additional macroscopic components compared with the previous order are the fourth-order derivatives $(\partial _{xxxx}\overline {\varphi }^{\text{(0)}},\partial _{xxyy}\overline {\varphi }^{\text{(0)}},\partial _{yyyy}\overline {\varphi }^{\text{(0)}})$ as well as the nonlinear term $\alpha {\overline {\eta }^{\text{(1)}}}({\partial \overline {\varphi }^{\text{(0)}}}/{\partial x})$ (see Appendix A.2). Accordingly, we set

(3.45) \begin{align} \begin{split} \varphi ^{\text{(4)}} & = \overline {\varphi }^{\text{(4)}}+\frac {\partial \overline {\varphi }^{\text{(3)}}}{\partial x}Q^{\text{(0)}}+\frac {\partial ^2 \overline {\varphi }^{\text{(2)}}}{\partial x ^2}Q^{\text{(1)}}_x+\frac {\partial ^2 \overline {\varphi }^{\text{(2)}}}{\partial y ^2}Q^{\text{(1)}}_y+\frac {\partial ^3 \overline {\varphi }^{\text{(1)}}}{\partial x ^3}Q^{\text{(2)}}_x+\frac {\partial ^3 \overline {\varphi }^{\text{(1)}}}{\partial x \partial y ^2}Q^{\text{(2)}}_y,\\[3pt] & \quad +\frac {\partial ^4\overline {\varphi }^{\text{(0)}}}{\partial x^4}Q^{\text{(3)}}_x +\frac {\partial ^4\overline {\varphi }^{\text{(0)}}}{\partial x^2\partial y^2}Q^{\text{(3)}}_{xy} +\frac {\partial ^4\overline {\varphi }^{\text{(0)}}}{\partial y^4}Q^{\text{(3)}}_y -\alpha {\overline {\eta }^{\text{(1)}}}\frac {\partial \overline {\varphi }^{\text{(0)}}}{\partial x}Q^{\text{(3)}} \end{split} \end{align}

(and ${\boldsymbol u}^{\text{(3)}}$ given in Appendix A.1) with $\overline {\varphi }^{\text{(4)}}$ defined in (3.17). The elementary problems satisfied by $(Q^{\text{(3)}}_{x},Q^{\text{(3)}}_{y},Q^{\text{(3)}}_{xy})$ are collected in Appendix B.1 since we shall not use them. Indeed we only need to know that they are even functions of $x_{{m}}$ ; the problem satisfied by $Q^{\text{(3)}}$ reads

(3.46) \begin{equation} \left \{\begin{array}{l} \displaystyle {\textrm {div}}_{\hspace {-.02cm}{m}} \big({\boldsymbol{\nabla }}_{\hspace {-.05cm}\textit{m}}Q^{\text{(3)}}\big)= 0 \quad \text{in } \varOmega ,\\[5pt] \displaystyle {\boldsymbol{\nabla }}_{\hspace {-.05cm}{m}}Q^{\text{(3)}} \boldsymbol{\cdot }{\boldsymbol{n}}=0 \text{ on the walls,}\\[7pt] \displaystyle \frac {\partial Q^{\text{(3)}}}{\partial {z_{{m}}}}=-\frac {\partial ^2 }{\partial x_{{m}} ^2}Q^{\text{(0)}}(x_{{m}},0), \quad \int _{X_{{m}}}Q^{\text{(3)}}\,\textrm {d}x_{{m}}=0, \quad \text{ at } {z_{{m}}}=0,\\[10pt] Q^{\text{(3)}}, {\boldsymbol{\nabla }}_{\hspace {-.05cm}{m}}Q^{\text{(3)}} p\text{-periodic with respect to } x_{{m}}, \end{array}\right . \end{equation}

and we notice that $Q^{\text{(3)}}$ is an odd function of $x_{{m}}$ , and in the absence of microstructure $Q^{\text{(0)}}=0$ , hence $Q^{\text{(3)}}=0$ . Integrating ${\boldsymbol u}^{\text{(3)}}={\boldsymbol{\nabla }}_{\hspace {-.05cm}\textit{m}}\varphi ^{\text{(4)}}+{\boldsymbol{\nabla }}_{\hspace {-.05cm}{m}}\varphi ^{\text{(3)}}$ in (3.44) over $\varOmega$ with (3.36) and (3.45) (see (A1) in Appendix A.1 for the detailed form of ${\boldsymbol u}^{\text{(3)}}$ ) and eliminating the contributions from odd functions, we obtain the form announced in (3.37b ) for $U_y^{\text{(3)}}$ and

(3.47) \begin{equation} U_x^{\text{(3)}}=\alpha _x\frac {\partial \overline {\varphi }^{\text{(3)}}}{\partial x}+d_{xx}\frac {\partial ^3 \overline {\varphi }^{\text{(1)}}}{\partial x ^3} +d_{xy}\frac {\partial ^3 \overline {\varphi }^{\text{(1)}}}{\partial x \partial y ^2} -\alpha d \,{\overline {\eta }^{\text{(1)}}}\frac {\partial \overline {\varphi }^{\text{(0)}}}{\partial x}, \end{equation}

with $d= {1}/{{S}}\int _\varOmega ({\partial Q^{\text{(3)}}}/{\partial x_{{m}}}) \,\textrm {d}{\boldsymbol r}_{{m}}$ , and it is shown in Appendix B.2 that $d=\kappa (1-n_x)$ .

3.6. Final model

The final model is obtained by gathering the results from order zero to order three, namely (3.19), (3.30) and (3.37), and considering in a first step the unique fields ${\boldsymbol U}={\boldsymbol U}^{\text{(0)}}+\delta {\boldsymbol U}^{\text{(1)}}+\boldsymbol{\cdots} +\delta ^3{\boldsymbol U}^{\text{(3)}}$ , $\overline {\eta }=\delta \overline {\eta }^{\text{(1)}}+\boldsymbol{\cdots} +\delta ^4\overline {\eta }^{\text{(4)}}$ . We obtain

\begin{align} \left \{\begin{array}{l} \displaystyle {\textrm {div}} {\boldsymbol U}+ \frac {\kappa }{\delta }\frac {\partial \overline {\eta }}{\partial t}+\alpha \delta ^2 \kappa {\textrm {div}}(\overline {\eta }\; {\boldsymbol{\nabla }}\,\overline {\varphi })=0,\kern164pt (3.48a)\\[6pt] \displaystyle \frac {\partial \overline {\varphi }}{\partial t}+\frac {1}{\delta }\overline {\eta }+\frac {\alpha \delta ^3}{2}\left ( n_x \left (\frac {\partial \overline {\varphi }}{\partial x}\right )^2+\left (\frac {\partial \overline {\varphi }}{\partial y}\right )^2 \right )=0, \\[10pt] \displaystyle U_x=\alpha _x\frac {\partial \overline {\varphi }}{\partial x}+\delta ^2d_{xx}\frac {\partial ^3 \overline {\varphi }}{\partial x ^3} +\delta ^2d_{xy}\frac {\partial ^3 \overline {\varphi }}{\partial x \partial y ^2} - \alpha \delta ^2 \kappa (1-n_x)\,{\overline {\eta }}\frac {\partial \overline {\varphi }}{\partial x},\qquad \qquad\quad (3.48b)\\[10pt] \displaystyle U_y=\frac {\partial \overline {\varphi }}{\partial y}+\delta ^2d_{yx}\frac {\partial ^3 \overline {\varphi }}{\partial x ^2\partial y} +\delta ^2d_{yy}\frac {\partial ^3 \overline {\varphi }}{\partial y ^3}.\end{array}\right.\nonumber \end{align}

By using the velocity at the free surface, $\overline {{\boldsymbol u}}= {\boldsymbol{\nabla }} \overline {\varphi }$ , instead of $\boldsymbol U$ , and eliminating $\overline {\varphi }$ , we obtain

(3.49) \begin{align} \displaystyle \frac {\partial \overline {\eta }}{\partial t}+\delta {\textrm {div}} \left [\left (\frac {\alpha _x}{\kappa }+ \varepsilon n_x \overline {\eta }\right ) \overline {u}_x\boldsymbol{e}_x + \left (\frac {1}{\kappa } + \varepsilon \overline {\eta }\right ) \overline {u}_y\boldsymbol{e}_y\right ] \nonumber\\[3pt]\displaystyle +\frac {\delta ^3}{\kappa } {\textrm {div}}\left [ \left (d_{xx}\frac {\partial ^2 \overline {u}_x}{\partial x ^2}+d_{xy}\frac {\partial ^2 \overline {u}_x}{\partial y ^2}\right )\boldsymbol{e}_x+ \left (d_{yx}\frac {\partial ^2 \overline {u}_y}{\partial x ^2} +d_{yy}\frac {\partial ^2 \overline {u}_y}{\partial y ^2}\right )\boldsymbol{e}_y \right ]=0, \\[-28pt] \nonumber \end{align}

(3.50) \begin{align} \frac {\partial \overline {\eta }}{\partial x}+\delta \frac {\partial \overline {u}_x}{\partial t}+\frac {\varepsilon \delta ^2}{2}\frac {\partial }{\partial x}\left (n_x \overline {u}_x^2+\overline {u}_y^2\right )=0, \quad \frac {\partial \overline {\eta }}{\partial y}+ \delta \frac {\partial \overline {u}_y}{\partial t}+\frac {\varepsilon \delta ^2}{2}\frac {\partial }{\partial y}\left (n_x \overline {u}_x^2+ \overline {u}_y^2\right )=0, \\[-6pt] \nonumber \end{align}

whose dimensional form is given in (2.7), with $\overline {h}={h}/\kappa$ being the mean water depth.

4. Analysis of the effective model

We investigate in this section the physical properties of the effective Boussinesq (2.7) and the KdV (2.8), as well as the soliton solution (2.10). To simplify the analysis, we consider a bathymetry composed of plates with vanishing thickness. By doing so, the average water depth $\overline {h}={h}$ remains fixed, leaving us with two dimensionless geometrical parameters: the array periodicity rescaled by the water depth column, $p$ , and the rescaled reduced water depth $\xi$ ( $\xi =1$ corresponds to a flat bathymetry). We will see that, in the considered limit of thin plates, the Boussinesq and KdV equations are significantly modified for wave propagation along the $x$ -direction, while they remain unchanged for propagation along the $y$ -direction. The case of oblique propagation continuously varies between these two extremes. This configuration also highlights a potential issue when applying homogenization directly to a reduced problem, written solely in horizontal variables and involving the depth profile $h(x)$ (Quezada de Luna & Ketcheson 2021; Ketcheson et al. 2025). In such an approach, if one simply considers a laminated depth profile $h(x)$ that alternates between $h$ and $\xi {h}$ over plates of vanishing thickness, one would incorrectly conclude using homogenization that these plates have no effect on wave propagation – which is not physically accurate in the $x$ -direction.

4.1. Analysis of the anisotropic Boussinesq (2.7)

According to (2.7), the effect of the bathymetry is encapsulated in the parameters $(\alpha _x,n_x,d_{xx})$ . For a solution depending only on $x$ , (2.7) simplifies to

(4.1) \begin{equation} \left \{\begin{array}{l} \displaystyle \displaystyle \frac {\partial \overline {\eta }}{\partial t} + \alpha _x\, {h}\frac {\partial \overline {u}_x}{\partial x} +n_x \frac {\partial }{\partial x}\overline {\eta }\overline {u}_x +{h}^3 d_{xx}\frac {\partial ^3 \overline {u}_x}{\partial x ^3}=0,\\[10pt] \displaystyle \frac {\partial \overline {u}_x}{\partial t}+g\frac {\partial \overline {\eta }}{\partial x}+n_x\overline {u}_x\frac {\partial \overline {u}_x}{\partial x}=0, \end{array}\right . \end{equation}

and for a flat bathymetry, we have $\alpha _x=n_x=1$ et $d_{xx}=1/3$ . In contrast, for a solution depending only on $y$ , (2.7) simplifies to

(4.2) \begin{equation} \left \{\begin{array}{l} \displaystyle \displaystyle \frac {\partial \overline {\eta }}{\partial t} + {h}\frac {\partial \overline {u}_y}{\partial y} + \frac {\partial }{\partial y}\overline {\eta } \overline {u}_y +{h}^3 d_{yy}\frac {\partial ^3 \overline {u}_y}{\partial y ^3}=0,\\[10pt] \displaystyle \frac {\partial \overline {u}_y}{\partial t}+g\frac {\partial \overline {\eta }}{\partial y}+\overline {u}_y\frac {\partial \overline {u}_y}{\partial y}=0, \end{array}\right . \end{equation}

which differs from the one-dimensional version of (2.4) for a flat bathymetry only in the value of $d_{yy}$ and its deviation from $1/3$ . The coefficients $(d_{xy}, d_{yx})$ appear only for oblique propagation in the $(x,y)$ plane.

Figures 4 and 5 show the variations of the six non-dimensional effective parameters, namely $(\alpha _x,n_x)$ and $(d_{xx},d_{xy},d_{yx},d_{yy})$ , which enter effective Boussinesq equations (2.7), see also (4.1) and (4.2). These parameters, defined by (2.12) and (2.13), were obtained by numerically solving, with COMSOL Multiphysics, the potential flow problems given in (2.11) together with table 1. In the numerical set-up, we considered an elementary unit cell (in dimensionless coordinates) of width $p$ and height 1, containing a plate of thickness $10^{-3}$ and height $(1-\xi )$ .

Figure 4. Parameters $(\alpha _x,n_x)$ entering the effective Boussinesq (2.7): (a) $Q^{\text{(0)}}(x_{{m}},{z_{{m}}})$ solution to (2.11) and table 1 ( $\xi =0.5$ , $p=0.25$ ); (b) ( $\alpha _x,n_x$ ) against the rescaled periodicity $p$ for increasing rescaled water depth over the plate $\xi =0.1$ to $0.5$ , with a step of 0.1).

Figure 5. Parameters $(d_{xx},d_{xy},d_{yx},d_{yy})$ entering the effective Boussinesq (2.7) against the rescaled periodicity $p$ for rescaled water depth over the plate $\xi =0.1$ to $0.5$ with a step of 0.1. The insets show the solutions to the elementary problems in the unit cell for $p=0.25$ and $\xi =0.5$ .

The deviation of $\alpha _x$ from 1 directly reflects the anisotropy of the medium, which applies both in the linear and nonlinear regimes: it quantifies the ratio of effective water depths in the $x$ - and $y$ -directions. In Maurel et al. (2017), it was shown that anisotropy is maximized for infinitely thin plates, which is the case we consider here, and for small spacing $p$ . We recover this result with $\alpha _x \to \xi$ in the limit $p \to 0$ (see table 3), meaning that the effective height in the $x$ -direction tends to $(\xi {h})$ , the smallest water depth imposed by the bathymetry. This result is intuitive since, when the spacing between plates is very small, the flow cannot penetrate the dense vertical layer region, which thus behaves like a flat bottom at depth $\xi {h}$ . This intuitive picture is in line with the analysis presented in Nachbin (2003), see in particular figures 5.1 and 6.2 of that reference. The effective height in the $y$ -direction remains $h$ , the largest water depth imposed by the bathymetry. For infinitely thin plates, this is also intuitive since the flow in the $y$ -direction does not ‘see’ the plates regardless of the periodicity $p$ , at least in the inviscid regime, and thus, only experiences the water column of depth $h$ .

Table 3. Estimates of the effective parameters for flat bathymetry in the limit of small periodicity to water depth ratio $p$ .

In the considered range of plate heights, variations in the nonlinearity parameter $n_x$ are relatively small, deviating only slightly from the value 1 corresponding to a flat bathymetry. Turning to the dispersion parameters, as mentioned earlier and as expected for infinitely thin plates, $d_{yy} \approx {1}/{3}$ regardless of $p$ and $\xi$ . The numerical computation of $d_{yy}$ retrieves this theoretical value with an accuracy of 0.3 %. In comparison, the dispersion parameters $(d_{xx}, d_{xy}, d_{yx})$ decrease significantly as $p$ and $\xi$ decrease. Notably, the effective heights $\alpha _x{h}$ and $(3 d_{xx})^{1/3} {h}$ coincide only for small values of $p$ .

These variations indicate that wave propagation along $y$ is unaffected by the presence of plates. In contrast, when the array is dense ( $p \to 0$ ), wave propagation along $x$ corresponds to propagation over the smallest water depth $\xi {h}$ imposed by the bathymetry, with $\alpha _x {h} = \xi {h}$ , $d_{xx} {h}^3 = ({1}/{3})(\xi {h})^3$ and $n_x=1$ in (4.1).

4.2. Analysis of the anisotropic KdV equation (2.8) and associated soliton solutions

We start by analysing the KdV equation. For wave propagation along $x$ ( $\theta =0$ ), (2.8) combined with (2.9) gives

(4.3) \begin{equation} \frac {\partial \overline {\eta }}{\partial t}+c_0\left (1+\frac {3}{2}\frac {\overline {\eta }}{h_0}\right )\frac {\partial \overline {\eta }}{\partial x}+\gamma _0\frac {c_0h_0^2}{6} \frac {\partial ^3 \overline {\eta }}{\partial x ^3}=0,\quad c_0=\sqrt {gH_0}, \end{equation}

where

(4.4) \begin{equation} \displaystyle H_0=\alpha _x{h}, \quad h_0=\frac {\alpha _x{h}}{n_x}, \quad \displaystyle \gamma _0 =\frac {3d_{xx}}{\alpha _x^3} n_x^2. \end{equation}

This equation differs from (2.5), which holds for a flat bathymetry, in two aspects: the relationship between $c_0$ and $h_0$ , since $H_0\neq h_0$ in general, and the introduction of the coefficient $\gamma _0$ , which deviates from unity. As a result, a potential imbalance between the dispersion and nonlinear terms arises. For wave propagation along $y$ , using $h_{({\pi}/{2})}={h}$ and $\gamma _{({\pi}/{2})}=3d_{yy}\simeq 1$ , we recover the classical KdV equation for a water depth $h$ , as expected.

For the general case of oblique incidence in (2.8), figure 6 displays the variations of the effective depths $(c_\theta /{c})^2=(H_\theta /{h})$ in the linear term and $(h_\theta /{h})$ in the dispersive term as functions of $p$ and $\theta$ . The results are shown for $\xi =0.3$ , but aside from bound values, the trends remain similar for different values of $\xi$ . Both effective depths exhibit similar variations, with $H_\theta \geqslant h_\theta$ . In agreement with the Boussinesq equation analysis, equality is reached for $\theta =0$ and small $p$ ( $H_0=h_0=\xi {h}$ ), as well as for $\theta =\pi /2$ regardless of $p$ ( $H_{({\pi}/{2})}=h_{({\pi }/{2})}={h}$ ). The parameter $\gamma _\theta$ exceeds 1, reaching its maximum at $\theta =0$ for intermediate values of $p$ . As $\xi$ decreases, this maximum increases and occurs at larger values of $p$ . Within the considered range of $\xi$ , the maximum is found around $p\sim 5-6$ and follows an empirical law of the form $1+0.04(1/\xi ^2-1)$ . We refrain from extrapolating the behaviour of the coefficients as $\xi \to 0$ , which would require revisiting the asymptotic analysis. It is worth noting that surface-piercing effects could occur for plates reaching the free surface, which falls beyond the scope of this study.

Figure 6. Anisotropic KdV (2.8)–(2.9): variations as functions of $p$ and $\theta$ of (a) the normalized effective depths $(c_\theta /c)^2=(H_\theta /{h})$ in the linear terms; (b) $h_\theta /h$ in the dispersion; (c) $\gamma _\theta$ in the nonlinear terms. The trends, shown for $\xi =0.3$ are similar for other values of $\xi$ , with $m\simeq 1+0.04(1/\xi ^2-1)$ representing the maximum of $\gamma _\theta$ obtained for $\theta =0$ .

We now examine the soliton solutions in (2.10), characterized by $\ell _\theta = ( {4\gamma _\theta h_\theta ^3}/{3\eta _{\scriptscriptstyle 0}} )^{1/2}$ (its spatial extent) and $u_\theta$ (its velocity). The reference case is given by $\ell = ( {4 {h}^3}/{3\eta _{\scriptscriptstyle 0}} )^{1/2}$ and $u$ in (2.6). Consistent with the KdV analysis, the soliton width is maximal for propagation along $y$ , where $\ell _{({\pi}/{2})}=\ell$ , and minimal for propagation along $x$ , with $\ell _0=\sqrt ({{3d_{xx}}/{n_x}})\ell$ . In the limit of dense plate arrays, we obtain $\ell _0 \sim \xi ^{3/2} \ell$ . To provide a concrete example, we consider a typical configuration with solitons over ${h}=10$ cm and $\eta _{\scriptscriptstyle 0}=0.3{h}$ (with $g=9.8$ m s $^{-2}$ ). Figure 7 presents the variations of these parameters with respect to $p$ for $\theta =0$ and with respect to $\theta$ for $p=1$ . The soliton velocity and exponential tails decrease significantly as $p$ or $\xi$ decreases, to the extent that even for moderate values of $\xi$ , these effects should be experimentally observable. This is illustrated in figure 8, which displays the soliton profiles for $\xi =0.5$ , 0.3 and 0.1 for propagation along $x$ (coloured curves) and along $y$ (grey curves). The profiles are shown at $t=0$ (solid lines) and $t=1$ s (dashed lines). Anisotropy in both soliton shape and velocity is already clearly visible for $\xi =0.5$ and becomes more pronounced as $\xi$ decreases. Note that these profiles correspond to the explicit solutions given in (2.10), rather than to a dynamical solution of the anisotropic KdV (2.8). However, as mentioned earlier, the stability of (2.8) for any propagation direction $\theta$ is guaranteed, since it has the same form as the classical KdV equation; therefore, the results of Benjamin (1972) and Bona (1975) apply directly.

Figure 7. Soliton characteristics: variations of $\ell _\theta$ (spatial extent) and $u_\theta$ (velocity) of solitons, (a–b) as functions of $p$ for $\theta =0$ , and (c–d) as functions of $\theta$ for $p=1$ . The parameters are set to ${h}=0.1$ m and $\eta _{\scriptscriptstyle 0}=0.03$ m.

Figure 8. Soliton profiles $\overline {\eta }(x,t)$ propagating along $x$ ( $\theta =0$ ), shown at $t=0$ (solid coloured lines) and $t=1$ s (dashed coloured lines), obtained from (2.10) for (a) $\xi =0.5$ , (b) $\xi =0.3$ and (c) $\xi =0.1$ . The grey curves represent the reference solution $\eta (x,t)$ from (2.6), which corresponds to propagation along $y$ , i.e. $\overline {\eta }(y,t)$ .

5. Concluding remarks

Following the pioneering work of Rosales & Papanicolaou (1983), we have derived effective two-dimensional Boussinesq and KdV equations encapsulating the effect of a periodic, rapidly varying, bathymetry in the $x$ -direction. These equations reveal strong anisotropy, which has been illustrated in the case where the bathymetry results from an array of thin plates, allowing us to explore the effects of the array periodicity and the plate height. Our deliberate choice of a relatively simple geometry, well-grounded in the metamaterials literature and easily realizable in a laboratory experiment (Berraquero et al. 2013; Maurel et al. 2017), enables us to emphasize effective parameters that remain robust in any configuration, as opposed to terms arising only from symmetry breaking. This simplification provides a clear foundation from which extensions to more complex bathymetries can be pursued within the same theoretical framework, with new terms naturally identified and their influence assessed. It also highlights the necessity of considering the three-dimensional effects of the original problem: any homogenization analysis starting from a reduced formulation based solely on a depth profile $h(x)$ fails to capture the plate-induced effects in the limit of vanishing thickness. Altogether, this emphasizes both the theoretical importance and the experimental feasibility of our approach.

The next step is to perform laboratory experiments, to confirm our results by investigating the ability of soliton-like solutions to propagate according to the derived laws and to explore the limitations of the approach, in particular the influence of vortex structures expected to be generated by the plates. From a numerical perspective, the stability of different formulations of the Boussinesq equation has been discussed in the case of a flat bottom, see for example Chen (1998) and Bona et al. (2002). A similar discussion can be carried out for the homogenized (2.7) we have obtained. Another promising direction would be to compare our homogenized results, for instance (2.10), with direct numerical simulations in the presence of bathymetry. In particular, the recent work of Svärd & Kalisch (2025) is highly relevant, as it establishes a Boussinesq system that is numerically stable even in the case of discontinuous bathymetries, including step configurations. This suggests a natural opportunity to confront the predictions of our homogenized model with those obtained from direct simulations over stepped bathymetries.