# Mechanical balance laws for two-dimensional Boussinesq systems

Chourouk El Hassanieh, Samer Israwi, Henrik Kalisch, Dimitrios Mitsotakis, Amutha Senthilkumar

<figure><img src="image_1.png" /></figure>

## To cite this version:

Chourouk El Hassanieh, Samer Israwi, Henrik Kalisch, Dimitrios Mitsotakis, Amutha Senthilkumar. Mechanical balance laws for two-dimensional Boussinesq systems. 2024. hal-04307057

HAL Id: hal-04307057

https://hal.science/hal-04307057v1

Preprint submitted on 22 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
MECHANICAL BALANCE LAWS FOR TWO DIMENSIONAL BOUSSINESQ
SYSTEMS

CHOUROUK EL HASSANIEH, SAMER ISRAWI, HENRIK KALISCH, DIMITRIOS MITSOTAKIS, AND AMUTHA SENTHILKUMAR

ABSTRACT. Most of the asymptotically derived Boussinesq systems of water wave theory for long waves of small amplitude fail to satisfy exact mechanical conservation laws for mass, momentum and energy. It is thus only fair to consider approximate conservation laws that hold in the context of these systems. Although such approximate mass, momentum and energy conservation laws can be derived, the question of a rigorous mathematical justification still remains unanswered. The aim of this paper is to justify the formally derived mechanical balance laws for weakly nonlinear and weakly dispersive water wave Boussinesq systems. In particular, two asymptotic expansions used for the formal and rigorous derivation of the Boussinesq systems and the same are employed for the derivation and rigorous justification of the balance laws. Numerical validation of the asymptotic orders of approximation is also presented.

# 1. INTRODUCTION

The propagation of surface water waves is described by the Euler equations of fluid mechanics, which are accompanied by dynamic boundary conditions at the free surface and the sea floor, as detailed in [21]. The solutions of the Euler equations also satisfy essential conservation equations, such as those for mass, momentum, and energy. While the Euler equations are a well-justified model both physically and mathematically, as indicated in [17], solving them, whether theoretically or numerically, remains a formidable challenge. This difficulty arises from the fact that the domain in which dependent variables like velocity and pressure are defined changes over time and is bounded by the unknown free surface.

As a result, numerous approximate models have been developed to approximate solutions of the Euler equations. These approximations rely on simplification assumptions about the characteristics of the waves to be described. These assumptions typically give rise to what are known as *wave regimes*. Two significant examples of water wave regimes are the small amplitude and long wave regime and the large amplitude, long wave regime. In both of these regimes, the primary focus is on long water waves, which means that the waves are considered to have a much greater wavelength compared to the depth of the water. In the small amplitude regime, the waves are assumed to have small amplitudes in comparison to the water depth, whereas in the large amplitude regime, there is no restriction on the amplitude of the waves.

In both regimes, formal approximations of the Euler equations are typically referred to as Boussinesq systems, with the exception of the large amplitude and long wave regime, where these systems are often termed Serre-Green-Naghdi equations, as mentioned in [17]. This nomenclature stems from the fact that the initial approximations were first formulated by J. Boussinesq [5] under the small amplitude assumption. Subsequently, the small amplitude assumption was removed by F. Serre [20] and independently by Green and Naghdi in [10], who derived equations describing strongly nonlinear and weakly dispersive water waves. For more comprehensive details, please refer to [17].

The initial approximations of the Euler equations had several theoretical and sometimes practical limitations. Therefore, research efforts were focused on deriving mathematically sound Boussinesq systems that possess favorable nonlinear and dispersive properties. Building upon Boussinesq's work, a whole class of Boussinesq-type systems were derived [3, 4] and were subsequently justified both theoretically and numerically, as documented in [16, 17, 18, 19]. Since these models find applications in the nearshore zone, it is imperative to be able to assess their accuracy concerning mass, momentum, and energy balance laws.

<footnote>Date: January 22, 2024.</footnote>

<footnote n="2000">2000 Mathematics Subject Classification. 35Q35, 76B15, 35B40.</footnote>

<footnote>Key words and phrases. Boussinesq systems, balance laws, asymptotic relations.</footnote>
In the one-dimensional case of the Serre equations, balance laws have been estimated using asymptotic techniques, and it has been established that these balance laws are satisfied by the solutions of the Serre equations with the same level of accuracy as they approximate the Euler equations, as discussed in [14]. Similarly, for one-dimensional Boussinesq systems, similar estimates were obtained, as documented in [1].

In this paper, we focus on a two-dimensional shallow water wave regime modeled by a set of Boussinesq-type equations initially presented in [4]. These two-dimensional Boussinesq Equations can be derived from the Water Waves problem through a straightforward asymptotic expansion based on the potential $\Phi$. They describe the propagation of long waves with small amplitudes on the surface of an ideal fluid. The motion of the free surface and the evolution of the velocity field of an incompressible, inviscid, and irrotational fluid under the influence of gravity is the actual notion of water waves on the surface of an ideal fluid.

After revisiting the derivation of the systems through asymptotic techniques, we conduct an asymptotic analysis of the mass, momentum, and energy balance laws. Theoretical justifications are provided for these balance laws, along with the approximations of the velocity field within the systems under consideration. These justifications are based on error estimates between solutions of the Euler equations and the Boussinesq systems, assuming identical initial data.

This paper is structured as follows: In Section 3, we provide an asymptotic derivation of the two-dimensional Boussinesq systems found in [4, 3]. In Section 4, we offer a mathematically rigorous derivation of the balance laws, expanding significantly on the preliminary work presented in [12, 11, ?]. Additionally, we conduct an analysis of the errors in the balance laws. We introduce essential theorems to enhance the understanding of the relationship between the formal derivation of the balance laws and their theoretical justification. These mechanical balance laws are dependent on the conventional small parameters $\alpha$ and $\beta$, which quantify the nonlinearity and dispersion of the system, respectively. We validate our theoretical findings numerically in Section 5. The numerically computed errors in mass, momentum, and energy conservation align with the theoretically derived estimates for the propagation of two-dimensional water waves. We employ linear analysis to elucidate the behavior of the solutions and their adherence to balance laws. Furthermore, our estimation of the solutions of the linearized Boussinesq equations reveals a distinct behavior of the two-dimensional solutions compared to what is known for one-dimensional problems, as documented in [21].

## 2. NOTATION

We start by introducing some notation. In what follows we denote by $X \in \mathbb{R}^2$ the horizontal variables $X = (x, y)$. We denote by $\alpha$ and $\beta$ the non-linearity and shallowness and parameters respectively given by (3.5). We use the following notations for the gradient and laplacian:

$$ \nabla = (\partial_x, \partial_y)^T, \quad \Delta = \partial_x^2 + \partial_y^2, \quad \nabla_{X,z} = (\partial_x, \partial_y, \partial_z)^T, \quad \nabla_{X,z}^{\beta} = (\sqrt{\beta} \nabla, \partial_z)^T. $$

We denote by $e_z$ the upward normal unit vector in the vertical direction, while $\partial_n u$ is the upward co-normal derivative of $u$. The non-dimensional fluid domain will be defined as

$$ \Omega_t = \{(X, z) \in \mathbb{R}^3, 0 \le z \le 1 + \alpha\eta(t, X)\}. $$

For all $a, b \in \mathbb{R}$ we write $a \vee b = \max\{a, b\}$. We denote by $C(\lambda_1, \lambda_2, \dots)$ a constant depending on the parameters $\lambda_1, \lambda_2, \dots$, and whose dependence on the $\lambda_j$ is always assumed to be non-decreasing.

Let $p$ be a constant with $1 \le p < \infty$ and $L^p(\mathbb{R}^2)$ the space of all Lebesgue-measurable functions $f$ with the standard norm

$$ |f|_{L^p(\mathbb{R}^2)} = \left( \int_{\mathbb{R}^2} |f(X)|^p dX \right)^{1/p} < \infty. $$

The space $L^\infty(\mathbb{R}^2)$ consists of all essentially bounded, Lebesgue-measurable functions $f$ with the norm

$$ |f|_{L^\infty(\mathbb{R}^2)} = \ess \sup |f(X)| < \infty. $$

Let $f: (-1, 0) \mapsto H^s(\mathbb{R}^2)$ be a Bochner measurable function. The space $L^\infty H^s = L^\infty((-1, 0); H^s(\mathbb{R}^2))$ is endowed with the canonical norm:

$$ |f|_{L^\infty H^s} = \ess \sup_{z \in (-1, 0)} |f(X, z)|_{H^s(\mathbb{R}^2)}. $$
For any real constant $s$, $H^s = H^s(\mathbb{R}^2)$ denotes the Sobolev space of all tempered distributions $f$ with the norm $|f|_{H^s} = |\Lambda^s f|_{L^2} < \infty$, where $\Lambda$ is the pseudo-differential operator $\Lambda = (1 - \Delta)^{1/2}$. We denote by $(H^{s,k}, |\cdot|_{H^{s,k}})$ the Banach space defined by

$$H^{s,k} = \bigcap_{j=0}^{k} H^j ((-1,0); H^{s-j}(\mathbb{R}^2)), \quad |u|_{H^{s,k}} = \sum_{j=0}^{k} |\Lambda^{s-j} \partial_z^j u|_{L^2}$$

We denote by $\dot{H}^{s+1}(\mathbb{R}^2)$ the topological vector space

$$\dot{H}^{s+1}(\mathbb{R}^2) = \{f \in L_{\text{loc}}^2(\mathbb{R}^2), \nabla f \in H^s(\mathbb{R}^2) \times H^s(\mathbb{R}^2)\},$$

equipped with the semi-norm $|f|_{\dot{H}^{s+1}(\mathbb{R}^2)} = |\nabla f|_{H^s(\mathbb{R}^2)}$.

### 3. FORMAL DERIVATION OF THE BOUSSINESQ SYSTEMS

We start with the derivation of the $a-b-c-d$ family of Boussinesq systems of [4]. We consider the domain $\{(x, y, z) \in \mathbb{R}^3, -h_0 < z < \eta(x, y, t)\}$ where the parameter $h_0$ represents the undisturbed depth of the fluid and $\eta(x, y, t)$ represents the free surface deviation above its rest position. If $u(x, y, z, t) = (u(x, y, z, t), v(x, y, z, t))$ denotes the horizontal components and $w(x, y, z, t)$ the vertical component of the fluid velocity vector field, the assuming an ideal and irrotational flow we can define the velocity potential $\Phi$ as

$$(3.1) \qquad \mathbf{u} = \nabla \Phi, \quad w = \Phi_z,$$

where $\nabla$ is the two-dimensional gradient operator $(\frac{\partial}{\partial x}, \frac{\partial}{\partial y})$. Then, the Bernoulli equation and the free-surface boundary condition governing the motion of the fluid are formulated in terms of the potential and the surface excursion by

$$(3.2) \qquad \left. \begin{array}{l} \Phi_t + \frac{1}{2}(\Phi_x^2 + \Phi_y^2 + \Phi_z^2) + g\eta = 0 \\ \eta_t + \Phi_x\eta_x + \Phi_y\eta_y - \Phi_z = 0 \end{array} \right\} \quad \text{on } z = \eta(x, y, t),$$

where $g$ is the gravitational acceleration, [21]. The derivation of asymptotic models relies on appropriate scaling. Here, we use non-dimensional variables:

$$(3.3a) \qquad \tilde{x} = \frac{x}{l}, \quad \tilde{y} = \frac{y}{l}, \quad \tilde{z} = \frac{z+h_0}{h_0}, \quad \tilde{t} = \frac{\sqrt{gh_0t}}{l},$$

$$(3.3b) \qquad \tilde{\eta} = \frac{\eta}{A}, \quad \tilde{\Phi} = \frac{h_0}{Al\sqrt{gh_0}}\Phi,$$

where tilde denotes non-dimensional variables, $l$ and $A$ denote a characteristic wavelength and wave amplitude. The governing equations and boundary conditions for the fully dispersive and fully non-linear irrotational water wave problem are given by the Euler equations:

$$(3.4a) \qquad \beta \nabla^2 \tilde{\Phi} + \tilde{\Phi}_{zz} = 0, \quad 0 < \tilde{z} < \alpha \tilde{\eta} + 1,$$

$$(3.4b) \qquad \tilde{\Phi}_{\tilde{z}} = 0, \quad \tilde{z} = 0,$$

$$(3.4c) \qquad \tilde{\Phi}_{\tilde{t}} + \frac{\alpha}{2} \left( \tilde{\Phi}_{\tilde{x}}^2 + \tilde{\Phi}_{\tilde{y}}^2 + \frac{1}{\beta} \tilde{\Phi}_{\tilde{z}}^2 \right) + \tilde{\eta} = 0, \quad \text{on } \tilde{z} = \alpha \tilde{\eta} + 1,$$

$$(3.4d) \qquad \tilde{\eta}_{\tilde{t}} + \alpha \left[ \tilde{\eta}_{\tilde{x}} \tilde{\Phi}_{\tilde{x}} + \tilde{\eta}_{\tilde{y}} \tilde{\Phi}_{\tilde{y}} \right] - \frac{1}{\beta} \tilde{\Phi}_{\tilde{z}} = 0, \quad \text{on } \tilde{z} = \alpha \tilde{\eta} + 1,$$

where $\alpha$ and $\beta$ are measures of nonlinearity and frequency dispersion defined as

$$(3.5) \qquad \alpha = \frac{A}{h_0}, \qquad \beta = \frac{h_0^2}{l^2}.$$

The derivation then follows by considering the following formal asymptotic series expansion of the velocity potential $\tilde{\Phi}$,

$$(3.6) \qquad \tilde{\Phi}(\tilde{x}, \tilde{y}, \tilde{z}, \tilde{t}) = \sum_{n=0}^{\infty} \tilde{z}^n \tilde{\Phi}^n(\tilde{x} \cdot \tilde{y}, \tilde{t}).$$
Substituting (3.6) into (3.4a) and the Neumann boundary condition at bottom (3.4b) we obtain a polynomial in $\tilde{z}$ and requiring the coefficient of each power of $\tilde{z}$ to vanish, we obtain the classical recurrence relation

$$ (3.7) \qquad \tilde{\Phi}^{n+2} = -\beta \frac{\nabla^2 \tilde{\Phi}^n}{(n+1)(n+2)}, \quad n=1,2,3,\dots $$

From (3.1), (3.6) and (3.7) one obtains the expression for velocity potential in the form

$$ (3.8) \qquad \tilde{\Phi}(\tilde{x}, \tilde{y}, \tilde{z}, \tilde{t}) = \sum_{n=0}^{\infty} (-1)^n \frac{\tilde{z}^{2n}}{(2n)!} \beta^n \nabla^{2n} \tilde{f}, $$

where $\tilde{f} = \tilde{\Phi}(\tilde{x}, \tilde{y}, 0, \tilde{t})$ is the potential evaluated at the bed. Substituting the above expression for $\tilde{\Phi}$ into (3.4c) yields the relation

$$ (3.9) \qquad \tilde{f}_{\tilde{t}} - \frac{\beta}{2} \Delta \tilde{f}_{\tilde{t}} + \tilde{\eta} + \frac{\alpha}{2} |\nabla \tilde{f}|^2 = \mathcal{O}(\alpha \beta, \beta^2). $$

Differentiating (3.9) with respect to $x$ and $y$ and expressing the equations in terms of the non-dimensional velocities at the bottom $\tilde{f}_{\tilde{x}} = \tilde{u}$, $\tilde{f}_{\tilde{y}} = \tilde{v}$, we have the following equations

$$ (3.10a) \qquad \tilde{u}_{\tilde{t}} + \tilde{\eta}_{\tilde{x}} - \frac{1}{2} \beta \Delta \tilde{u}_{\tilde{t}} + \alpha (\tilde{u} \tilde{u}_{\tilde{x}} + \tilde{v} \tilde{v}_{\tilde{x}}) = \mathcal{O}(\alpha \beta, \beta^2), $$

$$ (3.10b) \qquad \tilde{v}_{\tilde{t}} + \tilde{\eta}_{\tilde{y}} - \frac{1}{2} \beta \Delta \tilde{v}_{\tilde{t}} + \alpha (\tilde{u} \tilde{u}_{\tilde{y}} + \tilde{v} \tilde{v}_{\tilde{y}}) = \mathcal{O}(\alpha \beta, \beta^2), $$

$$ (3.10c) \qquad \tilde{\eta}_{\tilde{t}} + \tilde{u}_{\tilde{x}} + \tilde{v}_{\tilde{y}} - \frac{1}{6} \beta (\Delta \tilde{u}_{\tilde{x}} + \Delta \tilde{v}_{\tilde{y}}) + \alpha ((\tilde{\eta} \tilde{u})_{\tilde{x}} + (\tilde{\eta} \tilde{v})_{\tilde{y}}) = \mathcal{O}(\alpha \beta, \beta^2). $$

Now we let $\tilde{U}, \tilde{V}$ are the dimensionless velocities at a dimensionless height $\theta$ ($0 \le \theta \le 1$) in the fluid column. A formal use of Taylor's formula shows that

$$ (3.11a) \qquad \tilde{\Phi}_{\tilde{x}}|_{\tilde{z}=\theta} = \tilde{U} = \tilde{u} - \frac{\theta^2}{2} \beta \Delta \tilde{u} + \frac{\theta^4}{24} \beta^2 \Delta^2 \tilde{u} + \mathcal{O}(\beta^3), $$

$$ (3.11b) \qquad \tilde{\Phi}_{\tilde{y}}|_{\tilde{z}=\theta} = \tilde{V} = \tilde{v} - \frac{\theta^2}{2} \beta \Delta \tilde{v} + \frac{\theta^4}{24} \beta^2 \Delta^2 \tilde{v} + \mathcal{O}(\beta^3), $$

as $\beta \to 0$. Solving these equations for $\tilde{u}$ and $\tilde{v}$ we have that

$$ (3.12a) \qquad \tilde{u} = \left( 1 + \frac{\theta^2}{2} \beta \Delta + \frac{5\theta^4}{24} \beta^2 \Delta^2 \right) \tilde{U} + \mathcal{O}(\beta^3), $$

$$ (3.12b) \qquad \tilde{v} = \left( 1 + \frac{\theta^2}{2} \beta \Delta + \frac{5\theta^4}{24} \beta^2 \Delta^2 \right) \tilde{V} + \mathcal{O}(\beta^3). $$

We substitute the relations (3.12a) and (3.12b) into (3.10a), (3.10b) and (3.10c) to obtain

$$ (3.13a) \qquad \tilde{U}_{\tilde{t}} + \tilde{\eta}_{\tilde{x}} + \frac{\beta}{2}(\theta^2 - 1)\Delta\tilde{U}_{\tilde{t}} + \alpha (\tilde{U}\tilde{U}_{\tilde{x}} + \tilde{V}\tilde{V}_{\tilde{x}}) = \mathcal{O}(\alpha\beta, \beta^2), $$

$$ (3.13b) \qquad \tilde{V}_{\tilde{t}} + \tilde{\eta}_{\tilde{y}} + \frac{\beta}{2}(\theta^2 - 1)\Delta\tilde{V}_{\tilde{t}} + \alpha (\tilde{U}\tilde{U}_{\tilde{y}} + \tilde{V}\tilde{V}_{\tilde{y}}) = \mathcal{O}(\alpha\beta, \beta^2), $$

$$ (3.13c) \qquad \tilde{\eta}_{\tilde{t}} + \tilde{U}_{\tilde{x}} + \tilde{V}_{\tilde{y}} + \frac{\beta}{2} \left[ \theta^2 - \frac{1}{3} \right] (\Delta \tilde{U}_{\tilde{x}} + \Delta \tilde{V}_{\tilde{y}}) + \alpha ((\tilde{\eta} \tilde{U})_{\tilde{x}} + (\tilde{\eta} \tilde{V})_{\tilde{y}}) = \mathcal{O}(\alpha \beta, \beta^2). $$

Observing that

$$ (3.14) \qquad \tilde{U}_{\tilde{t}} + \tilde{\eta}_{\tilde{x}} = \mathcal{O}(\alpha, \beta), $$

$$ (3.15) \qquad \tilde{V}_{\tilde{t}} + \tilde{\eta}_{\tilde{y}} = \mathcal{O}(\alpha, \beta), $$

$$ (3.16) \qquad \tilde{\eta}_{\tilde{t}} + \tilde{U}_{\tilde{x}} + \tilde{V}_{\tilde{y}} = \mathcal{O}(\alpha, \beta), $$
we can write for all $\lambda, \mu \in \mathbb{R}$

$$ 
\begin{aligned}
\beta \Delta [\tilde{U}_{\tilde{x}} + \tilde{V}_{\tilde{y}}] &= \lambda\beta \Delta [\tilde{U}_{\tilde{x}} + \tilde{V}_{\tilde{y}}] + (1-\lambda)\beta \Delta [\tilde{U}_{\tilde{x}} + \tilde{V}_{\tilde{y}}] \\
&= \lambda\beta \Delta [\tilde{U}_{\tilde{x}} + \tilde{V}_{\tilde{y}}] - (1-\lambda)\beta \Delta \tilde{\eta}_{\tilde{t}} + \mathcal{O}(\alpha\beta, \beta^2),
\end{aligned}
$$

and also

$$ 
\begin{aligned}
\beta \Delta \tilde{U}_{\tilde{t}} &= \mu\beta \Delta \tilde{U}_{\tilde{t}} + (1-\mu)\beta \Delta \tilde{U}_{\tilde{t}} \\
&= \mu\beta \Delta \tilde{U}_{\tilde{t}} - (1-\mu)\beta \Delta \tilde{\eta}_{\tilde{x}} + \mathcal{O}(\alpha\beta, \beta^2),
\end{aligned}
$$

and similarly,

$$ \beta \Delta \tilde{V}_{\tilde{t}} = \mu\beta \Delta \tilde{V}_{\tilde{t}} - (1-\mu)\beta \Delta \tilde{\eta}_{\tilde{y}} + \mathcal{O}(\alpha\beta, \beta^2). $$

Substitution of these relations into (3.13a)-(3.13c) leads to the following general $a-b-c-d$ Boussinesq system

$$ (3.17a) \qquad \tilde{U}_{\tilde{t}} + \tilde{\eta}_{\tilde{x}} + \alpha (\tilde{U}\tilde{U}_{\tilde{x}} + \tilde{V}\tilde{V}_{\tilde{x}}) + \beta a \Delta \tilde{\eta}_{\tilde{x}} - \beta b \Delta \tilde{U}_{\tilde{t}} = \mathcal{O}(\alpha\beta, \beta^2), $$

$$ (3.17b) \qquad \tilde{V}_{\tilde{t}} + \tilde{\eta}_{\tilde{y}} + \alpha (\tilde{U}\tilde{U}_{\tilde{y}} + \tilde{V}\tilde{V}_{\tilde{y}}) + \beta a \Delta \tilde{\eta}_{\tilde{y}} - \beta b \Delta \tilde{V}_{\tilde{t}} = \mathcal{O}(\alpha\beta, \beta^2), $$

$$ (3.17c) \qquad \tilde{\eta}_{\tilde{t}} + \tilde{U}_{\tilde{x}} + \tilde{V}_{\tilde{y}} + \alpha ((\tilde{\eta}\tilde{U})_{\tilde{x}} + (\tilde{\eta}\tilde{V})_{\tilde{y}}) + \beta c \Delta (\tilde{U}_{\tilde{x}} + \tilde{V}_{\tilde{y}}) - \beta d \Delta \tilde{\eta}_{\tilde{t}} = \mathcal{O}(\alpha\beta, \beta^2). $$

where

$$ (3.18) \qquad \begin{aligned} a &= \frac{1}{2}(1 - \theta^2)\mu, & b &= \frac{1}{2}(1 - \theta^2)(1 - \mu), \\ c &= \frac{1}{2}\left(\theta^2 - \frac{1}{3}\right)\lambda, & d &= \frac{1}{2}\left(\theta^2 - \frac{1}{3}\right)(1 - \lambda). \end{aligned} $$

After neglecting the high-order terms and writing the variables in dimensional form system (3.17a)-(3.17c) is written as

$$ (3.19a) \qquad U_t + g\eta_x + (UU_x + VV_x) + gh_0^2 a \Delta \eta_x - h_0^2 b \Delta U_t = 0, $$

$$ (3.19b) \qquad V_t + g\eta_y + (UU_y + VV_y) + gh_0^2 a \Delta \eta_y - h_0^2 b \Delta V_t = 0, $$

$$ (3.19c) \qquad \eta_t + h_0 (U_x + V_y) + ((\eta U)_x + (\eta V)_y) + h_0^3 c \Delta (U_x + V_y) - h_0^2 d \Delta \eta_t = 0. $$

**Remark 1.** One should note that by introducing another parameter $\nu$, the system (3.17a),(3.17b),(3.17c) can be generalized to the $a-b-a_1-b_1-c-d$ system

$$ \tilde{U}_{\tilde{t}} + \tilde{\eta}_{\tilde{x}} + \alpha (\tilde{U}\tilde{U}_{\tilde{x}} + \tilde{V}\tilde{V}_{\tilde{x}}) + \beta a \Delta \tilde{\eta}_{\tilde{x}} - \beta b \Delta \tilde{U}_{\tilde{t}} = \mathcal{O}(\alpha\beta, \beta^2), $$

$$ \tilde{V}_{\tilde{t}} + \tilde{\eta}_{\tilde{y}} + \alpha (\tilde{U}\tilde{U}_{\tilde{y}} + \tilde{V}\tilde{V}_{\tilde{y}}) + \beta a_1 \Delta \tilde{\eta}_{\tilde{y}} - \beta b_1 \Delta \tilde{V}_{\tilde{t}} = \mathcal{O}(\alpha\beta, \beta^2), $$

$$ \tilde{\eta}_{\tilde{t}} + \tilde{U}_{\tilde{x}} + \tilde{V}_{\tilde{y}} + \alpha ((\tilde{\eta}\tilde{U})_{\tilde{x}} + (\tilde{\eta}\tilde{V})_{\tilde{y}}) + \beta c \Delta (\tilde{U}_{\tilde{x}} + \tilde{V}_{\tilde{y}}) - \beta d \Delta \tilde{\eta}_{\tilde{t}} = \mathcal{O}(\alpha\beta, \beta^2), $$

where

$$ a_1 = \frac{1}{2}(1 - \theta^2)\nu, \quad b_1 = \frac{1}{2}(1 - \theta^2)(1 - \nu). $$

Choosing distinct parameters $\mu$ and $\nu$ could be of interest when studying waves in the nearshore, where the dominant wave direction is approximately normal to the shoreline. However, in the present work, we will stick to the four-parameter system (3.17a),(3.17b),(3.17c) for the sake of tidiness.

In order to compute the associated mass, momentum and energy densities and fluxes, we need expressions for the velocities and pressure. The velocity field can be easily computed using (3.11a)-(3.11b). The expression for the pressure is obtained from Bernoulli's equation,

$$ (3.21) \qquad \Phi_t + \frac{1}{2} |\nabla_{X,z} \Phi|^2 = -\frac{P}{\rho} - gz + C. $$
We can find the constant C by evaluating the previous equation at the free surface $z = \eta$. Specifically, we find

$$ (3.22) \qquad C = \frac{P_{atm}}{\rho}, $$

where $P_{atm}$ refers to the atmospheric pressure. We introduce the dynamic pressure with the equation

$$ (3.23) \qquad P' = P - P_{atm} + \rho gz, $$

which can be scaled using a typical wave amplitude by $\rho g A \tilde{P}' = P'$, then

$$ \begin{aligned} \tilde{P}' &= -\tilde{\Phi}_{\tilde{t}} - \frac{1}{2}\alpha(\tilde{\Phi}_{\tilde{x}}^2 + \tilde{\Phi}_{\tilde{y}}^2) - \frac{1}{2}\frac{\alpha}{\beta}(\tilde{\Phi}_{\tilde{z}}^2) \\ &= -\tilde{f}_{\tilde{t}} + \beta\frac{\tilde{z}^2}{2}\Delta\tilde{f}_{\tilde{t}} - \frac{\alpha}{2}[\tilde{f}_{\tilde{x}}^2 + \tilde{f}_{\tilde{y}}^2] + \mathcal{O}(\alpha\beta, \beta^2). \end{aligned} $$

If we use the relation (3.9) and recall the relation $\tilde{f}_{\tilde{x}\tilde{x}\tilde{t}} + \tilde{f}_{\tilde{y}\tilde{y}\tilde{t}} = \tilde{U}_{\tilde{x}\tilde{t}} + \tilde{V}_{\tilde{y}\tilde{t}} + \mathcal{O}(\beta)$, then the scaled dynamic pressure becomes

$$ (3.24) \qquad \tilde{P}' = \tilde{\eta} + \frac{1}{2}\beta(\tilde{z}^2 - 1)[\tilde{U}_{\tilde{x}\tilde{t}} + \tilde{V}_{\tilde{y}\tilde{t}}] + \mathcal{O}(\alpha\beta, \beta^2). $$

The total pressure in-terms of dimensional variables is then given by

$$ (3.25) \qquad P = P_{atm} - \rho g (z - \eta) + \frac{\rho}{2} [(z + h_0)^2 - h_0^2] (U_{xt} + V_{yt}) + \mathcal{O}(\alpha\beta, \beta^2). $$

The following section is devoted to a mathematically rigorous approach to understanding the validity of the two-dimensional Boussinesq system as an approximation of the water-wave problem represented by the Euler equations (3.4) with a particular focus on justifying approximate mass, momentum and energy balance laws similar to those presented in [1], [2], [14].

## 4. RIGOROUS APPROACH TO BALANCE LAWS

As we showed in the previous section, appropriate assumptions on the respective magnitude of the parameters $\alpha$ and $\beta$, lead to the derivation of (simpler) asymptotic models from the Euler equations. The Stokes number

$$ S = \frac{\alpha}{\beta}, $$

is introduced in order to quantify the applicability of the equation to a particular regime of surface water waves. For the Boussinesq regime, the Stokes number is usually considered to be of order 1. Here, for the sake simplicity, we assume that the Stokes number is equal to 1 ($\alpha = \beta$), so that we can work with a single small parameter $\alpha$ or $\beta$. Throughout this section we will denote by $\mathcal{O}(\beta^n)$, with $n \in \mathbb{N}$ any family of functions $(f^\beta)_{\beta \in [0,1]}$ such that $\frac{1}{\beta^n}f^\beta$ remains bounded in $L^\infty([0, \frac{T}{\beta}], H^r)$, for all $\beta \ll 1$ and for possibly different values of $r$. It is noted that in the sequel we consider the Zakharov-Craig-Sulem formulation of the Euler equations. Specifically, we consider the Euler equations in terms of the Dirichlet-Neumann operator as follows :

$$ (4.1) \qquad \left\{ \begin{array}{l} \tilde{\eta}_{\tilde{t}} - \frac{1}{\beta} \mathcal{G}_{\beta}[\beta \tilde{\eta}] \psi = 0, \\ \psi_{\tilde{t}} + \tilde{\eta} + \frac{\beta}{2} |\nabla \psi|^2 - \frac{[\mathcal{G}_{\beta}[\beta \tilde{\eta}] \psi + \beta^2 \nabla \tilde{\eta} \cdot \nabla \psi]^2}{2(1 + \beta^3 |\nabla \tilde{\eta}|^2)} = 0. \end{array} \right. $$

Given a solution of this system, we reconstruct the potential $\tilde{\Phi}$ by solving the Laplace equation (4.2) below. More precisely, we introduce the trace of the velocity potential at the free surface, defined as

$$ \psi = \tilde{\Phi}_{|\tilde{z}=1+\beta\tilde{n}}, $$

and the Dirichlet-Neumann operator $\mathcal{G}_{\beta}[\beta \tilde{\eta}]$ as

$$ \mathcal{G}_{\beta}[\beta \tilde{\eta}] \psi = \partial_{\tilde{z}} \tilde{\Phi}_{|\tilde{z}=1+\beta \tilde{n}}, $$
with $\tilde{\Phi}$ solving the boundary value problem

$$ (4.2) \qquad \left\{ \begin{array}{l} \beta \partial_x^2 \tilde{\Phi} + \beta \partial_y^2 \tilde{\Phi} + \partial_z^2 \tilde{\Phi} = 0, \\ \partial_n \tilde{\Phi}_{|_{\tilde{z}=0}} = 0, \\ \tilde{\Phi}_{|_{\tilde{z}=1+\beta\tilde{\eta}}} = \psi. \end{array} \right. $$

We look for an asymptotic expansion of $\tilde{\Phi}$ of the form

$$ (4.3) \qquad \tilde{\Phi}^{app} = \sum_{j=0}^{N} \beta^j \tilde{\Phi}_j. $$

Plugging this expression into the boundary value problem (4.2) one can cancel the residual up to the order $O(\beta^{N+1})$ provided that

$$ (4.4) \qquad \partial_{\tilde{z}}^2 \tilde{\Phi}_j = -\partial_x^2 \tilde{\Phi}_{j-1} - \partial_y^2 \tilde{\Phi}_{j-1}, \quad j = 0, \dots, N, $$

(with the convention that $\tilde{\Phi}_{-1} = 0$), together with the boundary conditions

$$ (4.5) \qquad \left\{ \begin{array}{l} \tilde{\Phi}_{j|_{\tilde{z}=1+\beta\tilde{\eta}}} = \delta_{0,j}\psi, \\ \partial_{\tilde{z}}\tilde{\Phi}_{j|_{\tilde{z}=0}} = 0, \end{array} \right. \quad j = 0, \dots, N, $$

(where $\delta_{0,j} = 1$ if $j = 0$ and $0$ otherwise). Solving equation (4.4) with the boundary conditions (4.5) as in [16] one finds

$$ \tilde{\Phi}_0 = \psi, $$

$$ \tilde{\Phi}_1 = -\frac{1}{2}\tilde{z}^2\Delta\psi + \frac{1}{2}\Delta\psi + \beta\tilde{\eta}\Delta\psi + \frac{1}{2}\beta^2\tilde{\eta}^2\Delta\psi $$

$$ \tilde{\Phi}_2 = \frac{1}{24}\tilde{z}^4\Delta^2\psi - \frac{1}{4}\tilde{z}^2\Delta^2\psi + \frac{5}{24}\Delta^2\psi + \frac{5}{6}\beta\tilde{\eta}\Delta^2\psi - \frac{1}{2}\beta\tilde{z}^2\Delta\tilde{\eta}\Delta\psi + \frac{1}{2}\beta\Delta\tilde{\eta}\Delta\psi \\ - \beta\tilde{z}^2\nabla\tilde{\eta} \cdot \nabla\Delta\psi + \beta\nabla\tilde{\eta} \cdot \nabla\Delta\psi - \frac{1}{2}\beta\tilde{z}^2\tilde{\eta}\Delta^2\psi + O(\beta^2). $$

**4.1. Preliminary Results.** We first present results related to the formal approximations of the velocity potential. These results are necessary ingredients not only for the justification of the derivation of the Boussinesq systems but also of the justification of the mechanical balance laws that we will derive in the next section.

We denote by $\tilde{\Phi}^{Asym}$, the formally defined potential in Section 3, and we begin by decomposing $\tilde{\Phi}^{Asym}$ into two parts; a finite part $P^{series}$ and a remainder term $R^{series}$,

$$ \tilde{\Phi}^{Asym} := \sum_{n=0}^{\infty} (-1)^n \frac{\tilde{z}^{2n}}{(2n)!} \beta^n \nabla^{2n} f = P^{series} + R^{series}, $$

where

$$ (4.6) \qquad P^{series} := f - \beta \frac{\tilde{z}^2}{2} \Delta f + \beta^2 \frac{\tilde{z}^4}{24} \Delta^2 f, $$

and

$$ f := \tilde{\Phi}|_{\tilde{z}=0}. $$

We can consider the following choice of $f$:

$$ (4.7) \qquad f = \psi + \beta \frac{1}{2} \Delta \psi + \beta^2 \frac{5}{24} \Delta^2 \psi + \beta^2 \tilde{\eta} \Delta \psi. $$

**Remark 2.** Note that our choice for $f$ turns out to be well justified (up to order $\beta^3$) using the following definition:

$$ (4.8) \qquad \tilde{\Phi}^{app}|_{\tilde{z}=0} = \psi + \beta \frac{1}{2} \Delta \psi + \beta^2 \frac{5}{24} \Delta^2 \psi + \beta^2 \tilde{\eta} \Delta \psi + O(\beta^3). $$
Using the new notation $\tilde{u} = f_{\tilde{x}}$ and $\tilde{v} = f_{\tilde{y}}$, we easily verify the following expressions:

$$ (4.9) \qquad P_{\tilde{x}}^{\text{series}} = \tilde{u} - \frac{\tilde{z}^2}{2} \beta \Delta \tilde{u} + \frac{\tilde{z}^4}{24} \beta^2 \Delta^2 \tilde{u}, $$

and

$$ (4.10) \qquad P_{\tilde{y}}^{\text{series}} = \tilde{v} - \frac{\tilde{z}^2}{2} \beta \Delta \tilde{v} + \frac{\tilde{z}^4}{24} \beta^2 \Delta^2 \tilde{v}. $$

**Theorem 1.** Let $(\eta^{\text{Euler}}, \tilde{\Phi})$ be a regular solution of the Euler system (3.4), such that $(\eta^{\text{Euler}}, \nabla\psi) \in H^s(\mathbb{R}^2) \times H^s(\mathbb{R}^2)$ with $s$ large enough. Then, for $0 < \tilde{t} < T/\beta$ we have,

$$ |P^{\text{series}} - \tilde{\Phi}^{\text{app}}|_{L^{\infty}(\Omega_t)} \leq C\beta^3, $$

where $C$ is uniform with respect to the parameter $\beta$.

*Proof.* We denote by $\tilde{\eta} = \eta^{\text{Euler}}$, and we make use of the following identity

$$ \begin{aligned} P^{\text{series}} &= f - \beta \frac{\tilde{z}^2}{2} \Delta f + \beta^2 \frac{\tilde{z}^4}{24} \Delta^2 f \\ &= \psi + \beta \frac{1}{2} \Delta \psi + \beta^2 \frac{5}{24} \Delta^2 \psi + \beta^2 \tilde{\eta} \Delta \psi - \frac{1}{2} \beta \tilde{z}^2 \Delta \psi - \frac{1}{4} \beta^2 \tilde{z}^2 \Delta^2 \psi + \frac{1}{24} \beta^2 \tilde{z}^4 \Delta^2 \psi + \zeta \\ &= \tilde{\Phi}^{\text{app}} + h, \end{aligned} $$

where $h = \beta^3\zeta$ and $\zeta$ is a function of $\tilde{z}$, $\tilde{\eta}$ and the derivatives of $\psi$ resulting from the asymptotic expansion of the potential given by $\tilde{\Phi}^{\text{app}}$. Looking for solutions to the water-wave equations with $\psi$ in some Sobolev space would be too restrictive. Therefore, we can assume that the derivatives of $\psi$ belong to some Sobolev space $H^s(\mathbb{R}^2)$ with $s$ large enough so that $H^s(\mathbb{R}^2) \hookrightarrow L^\infty(\mathbb{R}^2)$. Using the fact that $\tilde{z}$ is bounded above by $1 + \beta\tilde{\eta}$ and $(\eta^{\text{Euler}}, \nabla\psi) \in H^s(\mathbb{R}^2) \times H^s(\mathbb{R}^2)$, we deduce that $h = \mathcal{O}(\beta^3)$, and hence the result follows. $\square$

Since we aim in finding a rigorous justification of the formally derived mechanical balance laws for two-dimensional Boussinesq systems, we must first estimate approximations between the velocities $\nabla_{X,z} P^{\text{series}}$ and $\nabla_{X,z} \tilde{\Phi}$ as well as for the pressure $\tilde{P}'$.

**4.2. Flattening the Domain.** The aim of this subsection is to present the method of "flattening the domain" of [17] so as to normalize the domain of the differential equations. Consider the boundary-value problem

$$ (4.11) \qquad \left\{ \begin{array}{l} \beta \partial_{\tilde{x}}^2 \omega + \beta \partial_{\tilde{y}}^2 \omega + \partial_{\tilde{z}}^2 \omega = R \quad \text{in } \Omega_t, \\ \omega|_{\tilde{z}=1+\beta\tilde{\eta}} = 0, \\ \partial_n \omega|_{\tilde{z}=0} = 0, \end{array} \right. $$

where $\omega$ and $R$ are smooth functions. Our goal is to give estimates on $\omega$ and $\omega_t$ that justify the results obtained in theorems 3 and 4. To do that we transform the above system into an elliptic boundary value problem on a fixed domain. We assume that the water height is always bounded from below. In other words,

$$ (4.12) \qquad \exists h_{\min} > 0, \quad \forall X = (x, y) \in \mathbb{R}^2, \quad 1 + \beta \tilde{\eta} \geq h_{\min}. $$

We transform the variable domain $\Omega_t$ into a flat strip $S = \mathbb{R}^2 \times (-1, 0)$ by introducing the following diffeomorphism<fnref n="1" /> (see Figure 1):

$$ \Sigma: \quad S \to \Omega_t \\ (X, z) \to \Sigma(X, z) = (X, (1 + \beta \tilde{\eta})z + 1 + \beta \tilde{\eta}) $$

<footnote n="1">¹We choose here the most obvious diffeomorphism but there are other choices of regularizing diffeomorphisms that are useful for obtaining optimal regularity estimates</footnote>
<figure><img src="image_3.png" /><figcaption>FIGURE 1. The transformation Σ in a two-dimensional setting</figcaption></figure>

Then, the new variables $\mathbf{w} = \omega \circ \Sigma$ and $\mathbf{R} = R \circ \Sigma$ satisfy the following elliptic boundary-value problem on the fixed domain $S$:

$$ (4.13) \qquad \begin{cases} \frac{1}{1 + \partial_{\tilde{z}} \sigma} \nabla_{X,z}^{\beta} \cdot P(\Sigma) \nabla_{X,z}^{\beta} \mathbf{w} = \mathbf{R} \quad \text{in } (-1,0) \times \mathbb{R}^2 \\ \mathbf{w}|_{\tilde{z}=0} = 0 \\ e_z \cdot P(\Sigma) \nabla_{X,z}^{\beta} \mathbf{w}|_{\tilde{z}=-1} = 0 \end{cases} $$

Here we use the notation

$$ \sigma(X, \tilde{z}) = \beta \tilde{\eta} \tilde{z} + 1 + \beta \tilde{\eta}, \quad P(\Sigma) = I + Q(\Sigma), \quad Q(\Sigma) = \begin{pmatrix} \partial_{\tilde{z}} \sigma I_2 & -\sqrt{\beta} \nabla \sigma \\ -\sqrt{\beta} \nabla \sigma T & \frac{-\partial_{\tilde{z}} \sigma + \beta |\nabla \sigma|^2}{1 + \partial_{\tilde{z}} \sigma} \end{pmatrix}. $$

**Theorem 2.** Let $\mathbf{R} \in H^{s,0}$, $\mathbf{w} \in H^{s+1,1}$, $s \ge 0$ solution of (4.13) such that assumption (4.12) is satisfied, then there exists a constant $C$ such that, for some $t_0 > 2$, the following estimate holds

$$ |\Lambda^s \nabla_{X,z}^{\beta} \mathbf{w}|_{L^2} \le C(h_{min}^{-1}, \beta_{max}, |\tilde{\eta}|_{H^{s+1 \vee t_0+1}}) |\Lambda^s \mathbf{R}|_{L^2}, $$

where $\beta_{max}$ is an upper bound of the shallowness parameter $\beta$.

*Proof.* We briefly state the results, first for the case $s = 0$. Multiplying the first equation of (4.13) by $\mathbf{w}$ and integrating yields

$$ \int_S \nabla_{X,z}^{\beta} \mathbf{w} \cdot P(\Sigma) \nabla_{X,z}^{\beta} \mathbf{w} = \int_S \mathbf{R} \mathbf{w}. $$

Using the coercivity<fnref n="2" /> of $P(\Sigma)$ along with Poincaré inequality one gets the following estimate:

$$ \begin{aligned} |\nabla_{X,z}^{\beta} \mathbf{w}|_{L^2}^2 &\le C(h_{min}^{-1}, \beta_{max}, |\tilde{\eta}|_{H^{t_0+2}}) |\mathbf{R}|_{L^2} |\mathbf{w}|_{L^2}, \\ &\le C |\mathbf{R}|_{L^2} |\nabla_{X,z}^{\beta} \mathbf{w}|_{L^2}, \end{aligned} $$

which proves the theorem for $s = 0$. The result follows by an induction relation for $s > 0$. For details of this argument, the reader may consult [17]. □

**Theorem 3.** Let $(\eta^{\text{Euler}}, \tilde{\Phi})$ be a regular solution of the Euler system such that $(\eta^{\text{Euler}}, \nabla\psi) \in H^s(\mathbb{R}^2) \times H^s(\mathbb{R}^2)$ with $s$ large enough. Assume that the total water depth satisfies (4.12). Then, for $0 < \tilde{t} < T/\beta$ we have,

$$ |\tilde{\Phi}_{\tilde{x}}^{app} - \tilde{\Phi}_{\tilde{x}}|_{L^{\infty}(\Omega_t)} \le C\beta^2, $$

$$ |\tilde{\Phi}_{\tilde{y}}^{app} - \tilde{\Phi}_{\tilde{y}}|_{L^{\infty}(\Omega_t)} \le C\beta^2, $$

$$ |\tilde{\Phi}_{\tilde{z}}^{app} - \tilde{\Phi}_{\tilde{z}}|_{L^{\infty}(\Omega_t)} \le C\beta^2, $$

where $C$ is a constant independent of $\beta$.

<footnote n="2">²The coercivity of $P(\Sigma)$ follows from the choice of the diffeomorphism, and we have $\forall \theta \in \mathbb{R}^3, \forall (X, z) \in S, \exists K > 0/P(\Sigma)(X, z)\theta. \theta \ge K|\theta|^2$</footnote>
Proof. Let $\omega = \tilde{\Phi}^{app} - \tilde{\Phi}$. Since $\tilde{\Phi}$ is a solution of (4.2) then its asymptotic expansion given by $\tilde{\Phi}^{app}$ satisfies

$$ egin{cases} \beta \partial_x^2 \tilde{\Phi}^{app} + \beta \partial_y^2 \tilde{\Phi}^{app} + \partial_z^2 \tilde{\Phi}^{app} = \beta^3 r, \\ \tilde{\Phi}^{app}|_{\tilde{z}=1+\beta\tilde{\eta}} = \psi, \\ \partial_n \tilde{\Phi}^{app}|_{\tilde{z}=0} = 0, \end{cases} $$

where $r$ is a regular function in terms of $\tilde{z}$ and the derivatives of $\psi$. Moreover, it is evident that

$$ \partial_n \omega|_{\tilde{z}=0} = 0 \quad \text{and} \quad \omega|_{\tilde{z}=1+\beta\tilde{\eta}} = 0. $$

Then $\omega$ satisfies the following boundary-value problem:

$$ (4.14) \qquad \left\{ \begin{array}{l} \beta \partial_x^2 \omega + \beta \partial_y^2 \omega + \partial_z^2 \omega = \beta^3 r, \\ \partial_n \omega|_{\tilde{z}=0} = 0, \\ \omega|_{\tilde{z}=1+\beta \tilde{\eta}} = 0. \end{array} \right. $$

Let $\Sigma$ be the diffeomorphism defined in section 4.2. Using $\Sigma$, we transform system (4.14) into a boundary value problem on a flat strip $S = (-1, 0) \times \mathbb{R}^2$ and hence it follows by theorem 2 that there exists a constant $C$ such that:

$$ |\Lambda^s \nabla_{X,z}^\beta \mathbf{w}|_{L^2} \le C \beta^3 |\Lambda^s \mathbf{r}|_{L^2}, $$

where $\mathbf{w} = \omega \circ \Sigma$ and $\mathbf{r} = r \circ \Sigma$. We shall use the fact that $H^{s-1,1}(S) \hookrightarrow L^\infty((-1,0); H^{s-3/2}(\mathbb{R}^2))$ with the following estimate

$$ \begin{aligned} |\nabla_{X,z}^\beta \mathbf{w}|_{H^{s-1,1}(S)} &= |\Lambda^{s-1} \nabla_{X,z}^\beta \mathbf{w}|_{L^2} + |\Lambda^{s-1} \nabla_{X,z}^\beta \partial_z \mathbf{w}|_{L^2} \\ &\le C(h_{min}^{-1}, \beta_{max}, |\tilde{\eta}|_{H^{s+1}}, |Q|_{L^\infty H^s}, |Q|_{H^{s,1}}, |\partial_z Q|_{L^\infty H^s}) |\Lambda^s \nabla_{X,z}^\beta \mathbf{w}|_{L^2}. \end{aligned} $$

Moreover, for $s$ large enough we have $H^{s-3/2}(\mathbb{R}^2) \hookrightarrow L^\infty(\mathbb{R}^2)$, therefore by using the above information we get

$$ \begin{aligned} |\nabla_{X,z}^\beta \mathbf{w}|_{L^\infty(S)} &= \underset{(X,z) \in S}{\operatorname{ess\,sup}} |\nabla_{X,z}^\beta \mathbf{w}| \\ &\le |\nabla_{X,z}^\beta \mathbf{w}|_{L^\infty((-1,0); H^{s-3/2}(\mathbb{R}^2))} \\ &\le C\beta^3. \end{aligned} $$

where $C$ is a constant independent of $\beta$. The last estimate remains true on $\Omega_t$ for $\nabla_{X,z}^\beta \omega$. Since, after changing variables we have

$$ \nabla_{X,z}^\beta \omega = \nabla_{X,z}^\beta (\mathbf{w} \circ \Sigma^{-1}) = J_{\Sigma^{-1}}^T (\nabla_{X,z}^\beta \mathbf{w} \circ \Sigma^{-1}), $$

where the coefficients of the Jacobian matrix $J_{\Sigma^{-1}}^T$ are bounded. $\square$

In the following corollary, we justify rigorously the equations (3.11a) and (3.11b),

**Corollary 1.** Let $(\eta^{Euler}, \tilde{\Phi})$ be a regular solution of the Euler system (3.4) such that $(\eta^{Euler}, \nabla\psi) \in H^s(\mathbb{R}^2) \times H^s(\mathbb{R}^2)$ with $s$ large enough. Then, for $0 < \tilde{t} < T/\beta$ we have,

$$ (4.15) \qquad |\tilde{\Phi}_{\tilde{x}} - P_{\tilde{x}}^{\text{series}}|_{L^{\infty}(\Omega_t)} \le C\beta^2, $$

$$ (4.16) \qquad |\tilde{\Phi}_{\tilde{y}} - P_{\tilde{y}}^{\text{series}}|_{L^{\infty}(\Omega_t)} \le C\beta^2, $$

$$ (4.17) \qquad |\tilde{\Phi}_{\tilde{z}} - P_{\tilde{z}}^{\text{series}}|_{L^{\infty}(\Omega_t)} \le C\beta^2. $$

Where $C$ is uniform with respect to the parameter $\beta$.

Proof. We proceed by using a triangle inequality followed by previous estimates. In fact, from Theorem 1 we have

$$ P_{\tilde{x}}^{\text{series}} - \tilde{\Phi}_{\tilde{x}}^{\text{app}} = h_{\tilde{x}}, $$
where $h_{\tilde{x}} = \beta^3 \zeta_{\tilde{x}}$ and $\zeta$ is a function of $\tilde{z}, \tilde{\eta}$ and the derivatives of $\psi$ regular enough. Therefore, for $s$ large enough such that $H^s(\mathbb{R}^2) \hookrightarrow L^\infty(\mathbb{R}^2)$ with $\tilde{z}$ bounded above by $1 + \beta\tilde{\eta}$, one can write

$$ \begin{aligned} |\tilde{\Phi}_{\tilde{x}} - P_{\tilde{x}}^{\text{series}}|_{L^{\infty}(\Omega_t)} &\le |\tilde{\Phi}_{\tilde{x}} - \tilde{\Phi}_{\tilde{x}}^{\text{app}}|_{L^{\infty}(\Omega_t)} + |\tilde{\Phi}_{\tilde{x}}^{\text{app}} - P_{\tilde{x}}^{\text{series}}|_{L^{\infty}(\Omega_t)} \\ &\le C\beta^2, \end{aligned} $$

where $C$ is uniform with respect to the parameter $\beta$. Similar considerations apply to the second and third inequalities. $\square$

**4.3. Expression for the pressure.** The aim of this section is to find an approximation of the pressure defined in the context of the Euler equations. We begin with a useful remark derived from the Zakharov-Craig-Sulem equations, cf. [18].

**Remark 3.** Let $\tilde{\eta} \in H^{s+1/2}(\mathbb{R}^2) \cap H^{t_0+2}(\mathbb{R}^2)$ with $s \ge 0$, $t_0 > 1$ satisfying (4.12). Then, the following mappings are continuous:

$$ (4.18) \qquad \mathcal{G}_{\beta}[\beta\tilde{\eta}] : \dot{H}^{s+1}(\mathbb{R}^2) \to H^{s-1/2}(\mathbb{R}^2) $$

$$ (4.19) \qquad \psi \mapsto \mathcal{G}_{\beta}[\beta\tilde{\eta}]\psi $$

$$ (4.20) \qquad \nu[\beta\tilde{\eta}] : \dot{H}^{s+1/2}(\mathbb{R}^2) \to H^{s-1/2}(\mathbb{R}^2) $$

$$ (4.21) \qquad \psi \mapsto \frac{[\mathcal{G}_{\beta}[\beta\tilde{\eta}]\psi + \beta^2 \nabla\tilde{\eta} \cdot \nabla\psi]^2}{2(1 + \beta^3 |\nabla\tilde{\eta}|^2)} $$

**Theorem 4.** Let $(\eta^{\text{Euler}}, \tilde{\Phi})$ be a regular solution of the Euler equations, such that $(\eta^{\text{Euler}}, \nabla\psi) \in H^s(\mathbb{R}^2) \times H^s(\mathbb{R}^2)$ with $s$ large enough. Assume that the total water depth satisfies (4.12). Then, for $0 < \tilde{t} < T/\beta$ we have,

$$ |\tilde{\Phi}_{\tilde{t}}^{\text{app}} - \tilde{\Phi}_{\tilde{t}}|_{L^{\infty}(\Omega_t)} \le C\beta^2, $$

where $C$ is independent of $\beta$.

*Proof.* Let $(\tilde{\eta}, \nabla\psi) \in H^s(\mathbb{R}^2) \times H^s(\mathbb{R}^2)$ with $s$ large enough. As a result of Remark 3 we have,

$$ \frac{[\mathcal{G}_{\beta}[\beta\tilde{\eta}]\psi + \beta^2 \nabla\tilde{\eta} \cdot \nabla\psi]^2}{2(1 + \beta^3 |\nabla\tilde{\eta}|^2)} \in H^{s-1/2}(\mathbb{R}^2). $$

In fact, this quantity appears in the expression of $\psi_{\tilde{t}}$ as shown in system (4.1):

$$ \psi_{\tilde{t}} = -\tilde{\eta} - \frac{\beta}{2} |\nabla\psi|^2 + \frac{[\mathcal{G}_{\beta}[\beta\tilde{\eta}]\psi + \beta^2 \nabla\tilde{\eta} \cdot \nabla\psi]^2}{2(1 + \beta^3 |\nabla\tilde{\eta}|^2)}. $$

Hence $\psi_{\tilde{t}} \in H^{s-1/2}(\mathbb{R}^2)$, and for $s$ large enough we have $H^{s-1/2}(\mathbb{R}^2) \hookrightarrow L^\infty(\mathbb{R}^2)$.

We will use the notation $\omega = \tilde{\Phi}^{\text{app}} - \tilde{\Phi}$. Then $\omega$ satisfies the boundary value problem (4.14). In fact, if we denote $\phi = \tilde{\Phi} \circ \Sigma$ then $\phi$ satisfies:

$$ (4.22) \qquad \left\{ \begin{array}{l} \nabla_{X,z}^{\beta} \cdot P(\Sigma) \nabla_{X,z}^{\beta} \phi = 0 \quad \text{in } \mathbb{R}^2 \times (-1, 0) \\ \phi|_{\tilde{z}=0} = \psi \\ e_z \cdot P(\Sigma) \nabla_{X,z}^{\beta} \phi|_{\tilde{z}=-1} = 0 \end{array} \right. $$

We look for an approximate solution to the above system of the form

$$ \phi^{\text{app}} = \sum_{j=0}^{N} \beta^j \phi_j. $$

By replacing $\phi^{\text{app}}$ in (4.22) and canceling higher order terms in $\beta$ one obtains:

$$ (4.23) \qquad \frac{1}{H} \partial_{\tilde{z}}^2 \phi_0 = 0, \qquad \phi_0|_{\tilde{z}=0} = \psi, \qquad \frac{1}{H} \partial_z \phi_0 = 0, $$
and for all $1 \le j \le n$,

$$ (4.24) \qquad \begin{cases} \frac{1}{H} \partial_{\tilde{z}}^2 \phi_j = -A(\nabla, \partial_{\tilde{z}}) \phi_{j-1}, \\ \phi_j|_{\tilde{z}=0} = \psi, \quad \frac{1}{H} \partial_{\tilde{z}} \phi_j = 0. \end{cases} $$

The operator $A(\nabla, \partial_z)$ is given by

$$ A(\nabla, \partial_{\tilde{z}})\bullet = \nabla \cdot (H\nabla\bullet) + \partial_{\tilde{z}}\left(\frac{|\nabla\sigma|^2}{H}\partial_{\tilde{z}}\bullet\right) - \nabla \cdot (\nabla\sigma \partial_{\tilde{z}}\bullet) - \partial_{\tilde{z}}(\nabla\sigma \cdot \nabla\bullet), $$

where $H(t, X) = 1 + \beta\tilde{\eta}(t, X)$ denotes the water depth and $\sigma(t, X, \tilde{z}) = \beta\tilde{\eta}(t, X)\tilde{z} + (1 + \beta\tilde{\eta}(t, X))$. Then, $\phi_0 = \psi$ and hence we can write:

$$ \phi = \psi + h, $$

where $h$ is function of $\tilde{z}$, $\tilde{\eta}$ and the derivatives of $\psi$. Thus, for $w = \omega \circ \Sigma$ one has $w_{\tilde{t}} \in H^{s-1/2}(\mathbb{R}^2)$. (One can choose initially $(\tilde{\eta}, \nabla\psi) \in H^{s+t_0}(\mathbb{R}^2) \times H^{s+t_0}(\mathbb{R}^2)$ with $t_0$ large enough.) Differentiating (4.13) with respect to time, and for the sake of brevity (we refer to Lemma 5.4 in [17]) one readily obtains the following estimate:

$$ |\Lambda^s \nabla_{X,z}^\beta \mathbf{w}_{\tilde{t}}|_{L^2} \le \beta^3 C(h_{\min}^{-1}, \beta_{max}, |\tilde{\eta}|_{H^{s+t_0}}, |\nabla\psi|_{H^{s+t_0}}). $$

Using the embedding $H^{s-1,1} \hookrightarrow L^\infty((-1,0); H^{s-3/2})$ as in Theorem 3 one obtains the result. $\square$

**Theorem 5.** Let $(\eta^{\text{Euler}}, \tilde{\Phi})$ be a regular solution of the Euler system such that $(\eta^{\text{Euler}}, \nabla\psi) \in H^s(\mathbb{R}^2) \times H^s(\mathbb{R}^2)$ with $s$ large enough. Then, for $0 < \tilde{t} < T/\beta$ we have,

$$ |P_{\tilde{t}}^{\text{series}} - \tilde{\Phi}_{\tilde{t}}|_{L^{\infty}(\Omega_t)} \le C\beta^2, $$

where $C$ is uniform with respect to the parameter $\beta$.

*Proof.* The following identity holds (we refer to the proof of Theorem 1 for clarification):

$$ P^{\text{series}} = \tilde{\Phi}^{\text{app}} + h, $$

where $h = \mathcal{O}(\beta^3)$ is a function of $\tilde{z}$, $\tilde{\eta}$ and the derivatives of $\psi$. Hence, we can write

$$ P_{\tilde{t}}^{\text{series}} = \tilde{\Phi}_{\tilde{t}}^{\text{app}} + h_{\tilde{t}}, $$

and we aim to show that $h_{\tilde{t}} \in H^{s-1/2}(\mathbb{R}^2)$ for $s$ large enough. In fact, from the second equation in system (4.1) we have

$$ \psi_{\tilde{t}} = -\tilde{\eta} - \frac{\beta}{2} |\nabla\psi|^2 + \frac{[G_{\beta}[\beta\tilde{\eta}]\psi + \beta^2\nabla\tilde{\eta} \cdot \nabla\psi]^2}{2(1 + \beta^3|\nabla\tilde{\eta}|^2)}. $$

Then, as a consequence of Remark 3, $\psi_{\tilde{t}} \in H^{s-1/2}(\mathbb{R}^2)$. Moreover, from the first equation in system (4.1) we have

$$ \tilde{\eta}_{\tilde{t}} = \frac{1}{\beta} G_{\beta}[\beta \tilde{\eta}] \psi, $$

which implies that $\tilde{\eta}_{\tilde{t}} \in H^{s-1/2}(\mathbb{R}^2)$. Choosing $s$ large enough such that $H^{s-1/2}(\mathbb{R}^2) \hookrightarrow L^\infty(\mathbb{R}^2)$ and for $\tilde{z}$ bounded from above by $1 + \beta\tilde{\eta}$, we can write:

$$ (4.25) \qquad |P_{\tilde{t}}^{\text{series}} - \tilde{\Phi}_{\tilde{t}}^{\text{app}}|_{L^{\infty}(\Omega_t)} \le C\beta^2, $$

where $C$ is independent of $\beta$. $\square$

The approximation of the pressure term $\tilde{P}'$ is given by the following corollary.

**Corollary 2.** Let $(\eta^{\text{Euler}}, \tilde{\Phi})$ be a regular solution of the Euler system such that $(\eta^{\text{Euler}}, \nabla\psi) \in H^s(\mathbb{R}^2) \times H^s(\mathbb{R}^2)$ with $s$ large enough. Then, for $0 < \tilde{t} < T/\beta$ we have,

$$ (4.26) \qquad |\tilde{P}' - \tilde{Q}|_{L^{\infty}(\Omega_t)} \le C(\alpha\beta + \beta^2), $$
where $\tilde{Q}$ is an approximation of the pressure defined by:

$$ (4.27) \qquad \tilde{Q} = -P_{\tilde{t}}^{\text{series}} - \frac{1}{2}\alpha \left( (P_{\tilde{x}}^{\text{series}})^2 + (P_{\tilde{y}}^{\text{series}})^2 \right) - \frac{1}{2}\frac{\alpha}{\beta}(P_{\tilde{z}}^{\text{series}})^2, $$

and $C$ is a constant independent of $\beta$.

*Proof.* In fact, the definition of the pressure $\tilde{P}'$ is given in section 3 by

$$ (4.28) \qquad \tilde{P}' = -\tilde{\Phi}_{\tilde{t}} - \frac{1}{2}\alpha (\tilde{\Phi}_{\tilde{x}}^2 + \tilde{\Phi}_{\tilde{y}}^2) - \frac{1}{2}\frac{\alpha}{\beta}\tilde{\Phi}_{\tilde{z}}^2, $$

and thus one has

$$ (4.29) \qquad \tilde{P}' - \tilde{Q} = -(\tilde{\Phi}_{\tilde{t}} - P_{\tilde{t}}^{\text{series}}) - \frac{1}{2}\alpha (\tilde{\Phi}_{\tilde{x}} - P_{\tilde{x}}^{\text{series}})(\tilde{\Phi}_{\tilde{x}} + P_{\tilde{x}}^{\text{series}}) - \frac{1}{2}\alpha (\tilde{\Phi}_{\tilde{y}} - P_{\tilde{y}}^{\text{series}})(\tilde{\Phi}_{\tilde{y}} + P_{\tilde{y}}^{\text{series}}) $$

$$ (4.30) \qquad -\frac{1}{2}\frac{\alpha}{\beta}(\tilde{\Phi}_{\tilde{z}} - P_{\tilde{z}}^{\text{series}})(\tilde{\Phi}_{\tilde{z}} + P_{\tilde{z}}^{\text{series}}). $$

Given that $\tilde{\Phi}$ is regular enough and choosing $s$ sufficiently large such that $H^s(\mathbb{R}^2) \hookrightarrow L^\infty(\mathbb{R}^2)$ implies that $\nabla_{X,z} P^{\text{series}} \in L^\infty(\Omega_t)$ and consequently

$$ \begin{aligned} |\nabla_{X,z} \tilde{\Phi} + \nabla_{X,z} P^{\text{series}}|_{L^\infty(\Omega_t)} &\le |\nabla_{X,z} \tilde{\Phi}|_{L^\infty(\Omega_t)} + |\nabla_{X,z} P^{\text{series}}|_{L^\infty(\Omega_t)} \\ &\le k, \end{aligned} $$

for some constant $k$. Hence using the previous estimates, one can deduce the result. $\square$

4.4. **Mass balance.** In this section, we establish the mass conservation properties of the general family of the two-dimensional Boussinesq systems (3.17a)–(3.17c), (3.18). The incompressibility of the fluid can be expressed in the form

$$ \rho_t + \frac{\partial}{\partial x}(\rho u) + \frac{\partial}{\partial y}(\rho v) + \frac{\partial}{\partial z}(\rho w) = 0, $$

which after integration and using Leibniz rule yields

$$ \begin{aligned} & \frac{\partial}{\partial t} \int_{-h_0}^{\eta} \rho \, dz - \rho \eta_t + \frac{\partial}{\partial x} \int_{-h_0}^{\eta} \rho u \, dz - \rho u|_{z=\eta} \eta_x \\ & \qquad + \frac{\partial}{\partial y} \int_{-h_0}^{\eta} \rho v \, dz - \rho v|_{z=\eta} \eta_y + \int_{-h_0}^{\eta} \frac{\partial}{\partial z}(\rho w) \, dz = 0. \end{aligned} $$

However, without loss of generality we assume that $\rho = 1$ which simplifies the previous relation. Since the vertical velocity at the bottom $w|_{(-h_0)} = 0$ and the kinematic boundary condition at the free surface of the water is $\eta_t + \Phi_x \eta_x + \Phi_y \eta_y - \Phi_z = 0$ on $z = \eta(x, y, t)$, we have

$$ (4.31) \qquad \frac{\partial}{\partial t} \int_{-h_0}^{\eta} dz + \frac{\partial}{\partial x} \int_{-h_0}^{\eta} u \, dz + \frac{\partial}{\partial y} \int_{-h_0}^{\eta} v \, dz = 0, $$

or equivalently in non-dimensional variables

$$ (4.32) \qquad \frac{\partial}{\partial \tilde{t}} (1 + \alpha \tilde{\eta}) + \frac{\partial}{\partial \tilde{x}} \int_{\tilde{z}=0}^{1+\alpha \tilde{\eta}} \alpha \tilde{\Phi}_{\tilde{x}} \, d\tilde{z} + \frac{\partial}{\partial \tilde{y}} \int_{\tilde{z}=0}^{1+\alpha \tilde{\eta}} \alpha \tilde{\Phi}_{\tilde{y}} \, d\tilde{z} = 0. $$

Substituting the expressions for $\tilde{\Phi}_{\tilde{x}}$ and $\tilde{\Phi}_{\tilde{y}}$ in terms of $\tilde{U}$ and $\tilde{V}$ gives

$$ (4.33) \qquad \begin{aligned} & \frac{\partial}{\partial \tilde{t}} (1 + \alpha \tilde{\eta}) + \frac{\partial}{\partial \tilde{x}} \left[ \tilde{U}(\alpha + \alpha^2 \tilde{\eta}) + \frac{\alpha \beta}{2} \left( \theta^2 - \frac{1}{3} \right) \Delta \tilde{U} \right] \\ & \qquad + \frac{\partial}{\partial \tilde{y}} \left[ \tilde{V}(\alpha + \alpha^2 \tilde{\eta}) + \frac{\alpha \beta}{2} \left( \theta^2 - \frac{1}{3} \right) \Delta \tilde{V} \right] = \mathcal{O}(\alpha \beta^2, \alpha^2 \beta). \end{aligned} $$

Finally, we obtain the differential balance equation

$$ (4.34) \qquad \tilde{\eta}_{\tilde{t}} + \tilde{U}_{\tilde{x}} + \tilde{V}_{\tilde{y}} + \alpha \left[ (\tilde{U}\tilde{\eta})_{\tilde{x}} + (\tilde{V}\tilde{\eta})_{\tilde{y}} \right] + \frac{\beta}{2} \left( \theta^2 - \frac{1}{3} \right) (\Delta \tilde{U}_{\tilde{x}} + \Delta \tilde{V}_{\tilde{y}}) = \mathcal{O}(\alpha \beta, \beta^2). $$
From (4.33) the non-dimensional mass density and the non-dimensional mass fluxes are

$$
\begin{align*}
\tilde{M} &= 1 + \alpha \tilde{\eta}, \\
\tilde{q}_{m_x} &= \tilde{U}(\alpha + \alpha^2 \tilde{\eta}) + \frac{\alpha \beta}{2} \left(\theta^2 - \frac{1}{3}\right) \Delta \tilde{U}, \\
\tilde{q}_{m_y} &= \tilde{V}(\alpha + \alpha^2 \tilde{\eta}) + \frac{\alpha \beta}{2} \left(\theta^2 - \frac{1}{3}\right) \Delta \tilde{V}.
\end{align*}
$$

Thus, the mass balance is

$$
(4.35) \qquad \frac{\partial}{\partial t} \tilde{M} + \frac{\partial}{\partial \tilde{x}} \tilde{q}_{m_x} + \frac{\partial}{\partial \tilde{y}} \tilde{q}_{m_y} = \mathcal{O}(\alpha \beta^2, \beta^3) .
$$

In dimensional variables the quantities in mass balance equation (4.35) are the following:

$$
\begin{align*}
M &= h_0 + \eta, \\
q_{m_x} &= U(h_0 + \eta) + h_0^3 \frac{1}{2} \left(\theta^2 - \frac{1}{3}\right) \Delta U, \\
q_{m_y} &= V(h_0 + \eta) + h_0^3 \frac{1}{2} \left(\theta^2 - \frac{1}{3}\right) \Delta V.
\end{align*}
$$

4.5. Momentum balance. In this section we find an approximate expression for momentum density and flux. for obtaining momentum balance, we first consider the Euler equation (4.36) and the incompressibility condition (4.37)

$$
(4.36) \qquad \vec{u}_t + (\vec{u} \cdot \nabla)\vec{u} + \nabla P = \vec{g},
$$

$$
(4.37) \qquad \nabla \cdot \vec{u} = 0,
$$

where $\vec{u} = (u, v, w)$ represents the velocity field, $P(x, y, z, t)$ represents the pressure and $\vec{g} = (0, 0, -g)$ represents gravitational vector. Writing the Euler equations in terms of the velocity potential $\vec{u} = \nabla \Phi$ we obtain the equations

$$
\Phi_{xt} + (\Phi_x^2)_x + (\Phi_x \Phi_y)_y + (\Phi_x \Phi_z)_z + P_x = 0,
$$

$$
\Phi_{yt} + (\Phi_y^2)_y + (\Phi_y \Phi_x)_x + (\Phi_y \Phi_z)_z + P_y = 0.
$$

Integrating over a fluid column and using the kinematic boundary condition (3.2) yields

$$
\frac{\partial}{\partial t} \int_{-h_0}^{\eta} \Phi_x dz + \frac{\partial}{\partial x} \int_{-h_0}^{\eta} (\Phi_x^2 + P) dz + \frac{\partial}{\partial y} \int_{-h_0}^{\eta} \Phi_x \Phi_y dz = 0,
$$

$$
\frac{\partial}{\partial t} \int_{-h_0}^{\eta} \Phi_y dz + \frac{\partial}{\partial y} \int_{-h_0}^{\eta} (\Phi_y^2 + P) dz + \frac{\partial}{\partial x} \int_{-h_0}^{\eta} \Phi_y \Phi_x dz = 0.
$$

Expressing the above relations in non-dimensional variables leads to the equations

$$
\alpha \frac{\partial}{\partial t} \int_{\tilde{z}=0}^{1+\alpha \tilde{\eta}} \tilde{\Phi}_{\tilde{x}} d\tilde{z} + \frac{\partial}{\partial \tilde{x}} \int_{\tilde{z}=0}^{1+\alpha \tilde{\eta}} \left\{ \alpha^2 (\tilde{\Phi}_{\tilde{x}}^2) + \alpha \tilde{P}' - (\tilde{z}-1) \right\} d\tilde{z} + \alpha^2 \frac{\partial}{\partial \tilde{y}} \int_{\tilde{z}=0}^{1+\alpha \tilde{\eta}} \tilde{\Phi}_{\tilde{x}} \tilde{\Phi}_{\tilde{y}} d\tilde{z} = 0,
$$

$$
\alpha \frac{\partial}{\partial t} \int_{\tilde{z}=0}^{1+\alpha \tilde{\eta}} \tilde{\Phi}_{\tilde{y}} d\tilde{z} + \frac{\partial}{\partial \tilde{y}} \int_{\tilde{z}=0}^{1+\alpha \tilde{\eta}} \left\{ \alpha^2 (\tilde{\Phi}_{\tilde{y}}^2) + \alpha \tilde{P}' - (\tilde{z}-1) \right\} d\tilde{z} + \alpha^2 \frac{\partial}{\partial \tilde{x}} \int_{\tilde{z}=0}^{1+\alpha \tilde{\eta}} \tilde{\Phi}_{\tilde{x}} \tilde{\Phi}_{\tilde{y}} d\tilde{z} = 0.
$$

Substituting non-dimensional velocity potentials $\tilde{\Phi}_{\tilde{x}}$ and $\tilde{\Phi}_{\tilde{y}}$ and $\tilde{P}'$ in terms of $\tilde{U}$ and $\tilde{V}$ gives the momentum balance equations

$$
(4.38) \quad \frac{\partial}{\partial \tilde{t}} \left\{ (1 + \alpha \tilde{\eta}) \tilde{U} + \frac{\beta}{2} \left( \theta^2 - \frac{1}{3} \right) \Delta \tilde{U} \right\} + \frac{\partial}{\partial \tilde{x}} \left\{ \tilde{\eta} + \alpha \tilde{U}^2 + \frac{\alpha}{2} \tilde{\eta}^2 - \frac{1}{3} \beta (\tilde{U}_{\tilde{x}\tilde{t}} + \tilde{V}_{\tilde{y}\tilde{t}}) + \frac{1}{2} \right\} \\
\qquad \qquad + \frac{\partial}{\partial \tilde{y}} (\alpha \tilde{U} \tilde{V}) = \mathcal{O}(\alpha \beta, \beta^2),
$$
$$ (4.39) \qquad \frac{\partial}{\partial \tilde{t}} \left\{ \left(1 + \alpha \tilde{\eta}\right) \tilde{V} + \frac{\beta}{2} \left(\theta^2 - \frac{1}{3}\right) \Delta \tilde{V} \right\} + \frac{\partial}{\partial \tilde{y}} \left\{ \tilde{\eta} + \alpha \tilde{V}^2 + \frac{\alpha}{2} \tilde{\eta}^2 - \frac{1}{3} \beta (\tilde{U}_{\tilde{x}\tilde{t}} + \tilde{V}_{\tilde{y}\tilde{t}}) + \frac{1}{2} \right\} \\ \qquad \qquad \qquad \qquad + \frac{\partial}{\partial \tilde{x}} (\alpha \tilde{U} \tilde{V}) = \mathcal{O}(\alpha \beta, \beta^2). $$

If the terms of order $\mathcal{O}(\alpha\beta, \beta^2)$ are neglected, then the momentum balance equations written in dimensional variables take the following form

$$ \frac{\partial}{\partial t} \left\{ (h_0 + \eta)U + \frac{1}{2} \left(\theta^2 - \frac{1}{3}\right) \Delta U \right\} + \frac{\partial}{\partial x} \left\{ h_0 U^2 + \frac{g}{2} (h_0 + \eta)^2 - \frac{h_0^3}{3} (U_{xt} + V_{yt}) \right\} + \frac{\partial}{\partial y} (h_0 U V) = 0, $$

$$ \frac{\partial}{\partial t} \left\{ (h_0 + \eta)V + \frac{1}{2} \left(\theta^2 - \frac{1}{3}\right) \Delta V \right\} + \frac{\partial}{\partial y} \left\{ h_0 V^2 + \frac{g}{2} (h_0 + \eta)^2 - \frac{h_0^3}{3} (U_{xt} + V_{yt}) \right\} + \frac{\partial}{\partial x} (h_0 U V) = 0. $$

4.6. **Energy balance.** The exact energy balance equation reads

$$ (4.40) \qquad \frac{\partial}{\partial t} \left\{ \frac{1}{2} |\nabla \Phi^2| + gz \right\} + \nabla \cdot \left\{ \left( \frac{1}{2} |\nabla \Phi^2| + gz + P \right) \nabla \Phi \right\} = 0. $$

Integrating over the a water column yields

$$ \frac{\partial}{\partial t} \left\{ \int_{-h_0}^{\eta} \frac{1}{2} |\nabla \Phi|^2 dz + \int_{0}^{\eta} gz dz \right\} + \frac{\partial}{\partial x} \left\{ \int_{-h_0}^{\eta} \left( \frac{1}{2} |\nabla \Phi|^2 + gz + P \right) \Phi_x \right\} dz \\ + \frac{\partial}{\partial y} \left\{ \int_{-h_0}^{\eta} \left( \frac{1}{2} |\nabla \Phi|^2 + gz + P \right) \Phi_y \right\} dz = 0. $$

Using non-dimensional variables the last equation reads

$$ \frac{\partial}{\partial \tilde{t}} \left\{ \int_{\tilde{z}=0}^{1+\alpha\tilde{\eta}} \frac{\alpha^2}{2} \left( \tilde{\Phi}_{\tilde{x}}^2 + \tilde{\Phi}_{\tilde{y}}^2 + \frac{1}{\beta} \tilde{\Phi}_{\tilde{z}}^2 \right) d\tilde{z} + \int_{\tilde{z}=1}^{1+\alpha\tilde{\eta}} (\tilde{z}-1) d\tilde{z} \right\} \\ + \frac{\partial}{\partial \tilde{x}} \int_{\tilde{z}=0}^{1+\alpha\tilde{\eta}} \left\{ \frac{\alpha^3}{2} \left( \tilde{\Phi}_{\tilde{x}}^2 + \tilde{\Phi}_{\tilde{y}}^2 + \frac{1}{\beta} \tilde{\Phi}_{\tilde{z}}^2 \right) + \alpha(\tilde{z}-1) + \alpha^2 \tilde{P}' + \alpha(1-\tilde{z}) \right\} \tilde{\Phi}_{\tilde{x}} d\tilde{z} \\ + \frac{\partial}{\partial \tilde{y}} \int_{\tilde{z}=0}^{1+\alpha\tilde{\eta}} \left\{ \frac{\alpha^3}{2} \left( \tilde{\Phi}_{\tilde{x}}^2 + \tilde{\Phi}_{\tilde{y}}^2 + \frac{1}{\beta} \tilde{\Phi}_{\tilde{z}}^2 \right) + \alpha(\tilde{z}-1) + \alpha^2 \tilde{P}' + \alpha(1-\tilde{z}) \right\} \tilde{\Phi}_{\tilde{y}} d\tilde{z} = 0. $$

Substituting the expressions for $\tilde{\Phi}_{\tilde{x}}$ and $\tilde{\Phi}_{\tilde{y}}$ in terms of $\tilde{U}$ and $\tilde{V}$ and (3.24) leads to the energy balance equation

$$ \frac{\partial}{\partial \tilde{t}} \left[ \frac{1}{2}(\tilde{U}^2 + \tilde{V}^2 + \tilde{\eta}^2) + \frac{\beta}{2} \left(\theta^2 - \frac{1}{3}\right) (\tilde{U} \Delta \tilde{U} + \tilde{V} \Delta \tilde{V}) + \frac{\beta}{6} (\tilde{U}_{\tilde{x}} + \tilde{V}_{\tilde{y}})^2 + \frac{\alpha}{2} \tilde{\eta} (\tilde{U}^2 + \tilde{V}^2) \right] \\ + \frac{\partial}{\partial \tilde{x}} \left[ \frac{\alpha}{2} (\tilde{U}^3 + \tilde{V}^2 \tilde{U}) + \alpha \tilde{\eta}^2 \tilde{U} + \tilde{\eta} \tilde{U} + \frac{\beta}{2} \left(\theta^2 - \frac{1}{3}\right) \tilde{\eta} \Delta \tilde{U} - \frac{\beta}{3} \tilde{U} (\tilde{U}_{\tilde{x}\tilde{t}} + \tilde{V}_{\tilde{y}\tilde{t}}) \right] \\ + \frac{\partial}{\partial \tilde{y}} \left[ \frac{\alpha}{2} (\tilde{V}^3 + \tilde{U}^2 \tilde{V}) + \alpha \tilde{\eta}^2 \tilde{V} + \tilde{\eta} \tilde{V} + \frac{\beta}{2} \left(\theta^2 - \frac{1}{3}\right) \tilde{\eta} \Delta \tilde{V} - \frac{\beta}{3} \tilde{V} (\tilde{U}_{\tilde{x}\tilde{t}} + \tilde{V}_{\tilde{y}\tilde{t}}) \right] = \mathcal{O}(\alpha\beta, \beta^2). $$

Hence, the general form of the energy balance equation is

$$ (4.41) \qquad \frac{\partial}{\partial \tilde{t}} \tilde{E} + \frac{\partial}{\partial \tilde{x}} \tilde{q}_{E_x} + \frac{\partial}{\partial \tilde{y}} \tilde{q}_{E_y} = \mathcal{O}(\alpha\beta, \beta^2). $$

The dimensional form of the quantities are

$$ (4.42) \qquad E = \frac{1}{2}h_0 (U^2 + V^2) + \frac{1}{2}h_0^3 \left(\theta^2 - \frac{1}{3}\right) (U\Delta U + V\Delta V) \\ + \frac{1}{2} (U^2 + V^2) \eta + \frac{1}{6}h_0^3 (U_x + V_y)^2 + \frac{1}{2}g\eta^2, $$

and

$$ (4.43) \qquad q_{E_x} = \frac{1}{2} h_0 (U^3 + U V^2) + g h_0 \eta U + \frac{1}{2} g h_0^3 \left(\theta^2 - \frac{1}{3}\right) \eta \Delta U - \frac{1}{3} h_0^3 U (U_{xt} + V_{yt}) + g U \eta^2, $$
$$ (4.44) \qquad q_{E_y} = \frac{1}{2}h_0(V^3 + U^2V) + gh_0\eta V + \frac{1}{2}gh_0^3\left(\theta^2 - \frac{1}{3}\right)\eta\Delta V - \frac{1}{3}h_0^3V(U_{xt} + V_{yt}) + gV\eta^2. $$

4.7. On the role of the rigorous approach. Here we justify the terms $O(\alpha\beta, \beta^2)$ and $O(\alpha\beta^2, \beta^3)$ in the mechanical balance laws by justifying the formal use of Taylor's formula in (3.11a) and (3.11b). We first verify that the expressions (3.13a), (3.13b), and (3.13c) hold for an explicit definition of the variable $f$ given by (4.7).

We focus on finding a common ground between the expressions used in Section 3 with new similar expressions satisfied by $P_{\tilde{x}}^{series}$ and $P_{\tilde{y}}^{series}$.

In fact, from equations (4.9) and (4.10) one can write up to order $O(\beta^3)$:

$$ (4.45) \qquad P_{\tilde{x}}^{series}|_{\tilde{z}=\theta} = \tilde{U} = \tilde{u} - \frac{\theta^2}{2}\beta\Delta\tilde{u} + \frac{\theta^4}{24}\beta^2\Delta^2\tilde{u}, $$

$$ (4.46) \qquad P_{\tilde{y}}^{series}|_{\tilde{z}=\theta} = \tilde{V} = \tilde{v} - \frac{\theta^2}{2}\beta\Delta\tilde{v} + \frac{\theta^4}{24}\beta^2\Delta^2\tilde{v}. $$

One can verify that $f$ satisfies (3.9) and thus (3.13a), (3.13b), and (3.13c) hold in terms of the variables $\tilde{U}$ and $\tilde{V}$ which are now well justified in terms of $f$ and its derivatives. Therefore, giving a rigorous justification to the formally applied Taylor's formula used in (3.11a) and (3.11b). This would not have been possible without using the fact that the approximations in Section 4.1 hold, specifically Corollary 1 and Theorem 5 that bring together the results and ensure that the use of $\nabla P^{series}$, $P_{\tilde{z}}^{series}$ and $P_{\tilde{t}}^{series}$ to approximate $\nabla\tilde{\Phi}$, $\tilde{\Phi}_z$ and $\tilde{\Phi}_t$ is valid up to order $O(\beta^2)$.

**Theorem 6.** Let $(\eta^{Euler}, \tilde{\Phi})$ be a regular solution of the Euler system such that $(\eta^{Euler}, \nabla\psi) \in H^s(\mathbb{R}^2) \times H^s(\mathbb{R}^2)$ with $s$ large enough. Then, there exists constants $C_1, C_2, C_3, C_4$ independent of $\beta$ such that the following approximate mechanical balance laws are satisfied:

$$ \begin{aligned} & \left| \frac{\partial}{\partial \tilde{t}} (1 + \alpha \tilde{\eta}) + \frac{\partial}{\partial \tilde{x}} \left\{ \tilde{U}(\alpha + \alpha^2 \tilde{\eta}) + \frac{\alpha \beta}{2} \left( \theta^2 - \frac{1}{3} \right) \Delta \tilde{U} \right\} \right. \\ & \qquad \left. + \frac{\partial}{\partial \tilde{y}} \left\{ \tilde{V}(\alpha + \alpha^2 \tilde{\eta}) + \frac{\alpha \beta}{2} \left( \theta^2 - \frac{1}{3} \right) \Delta \tilde{V} \right\} \right|_{L^\infty(\mathbb{R}^2)} \le C_1 (\alpha \beta^2 + \alpha^2 \beta), \\ & \left| \frac{\partial}{\partial \tilde{t}} \left\{ (1 + \alpha \tilde{\eta}) \tilde{U} + \frac{\beta}{2} \left( \theta^2 - \frac{1}{3} \right) \Delta \tilde{U} \right\} + \frac{\partial}{\partial \tilde{x}} \left\{ \tilde{\eta} + \alpha \tilde{U}^2 + \frac{\alpha}{2} \tilde{\eta}^2 - \frac{1}{3} \beta (\tilde{U}_{\tilde{x}\tilde{t}} + \tilde{V}_{\tilde{y}\tilde{t}}) + \frac{1}{2} \right\} \right. \\ & \qquad \left. + \frac{\partial}{\partial \tilde{y}} (\alpha \tilde{U} \tilde{V}) \right|_{L^\infty(\mathbb{R}^2)} \le C_2 (\alpha \beta + \beta^2), \\ & \left| \frac{\partial}{\partial \tilde{t}} \left\{ (1 + \alpha \tilde{\eta}) \tilde{V} + \frac{\beta}{2} \left( \theta^2 - \frac{1}{3} \right) \Delta \tilde{V} \right\} + \frac{\partial}{\partial \tilde{y}} \left\{ \tilde{\eta} + \alpha \tilde{V}^2 + \frac{\alpha}{2} \tilde{\eta}^2 - \frac{1}{3} \beta (\tilde{U}_{\tilde{x}\tilde{t}} + \tilde{V}_{\tilde{y}\tilde{t}}) + \frac{1}{2} \right\} \right. \\ & \qquad \left. + \frac{\partial}{\partial \tilde{x}} (\alpha \tilde{U} \tilde{V}) \right|_{L^\infty(\mathbb{R}^2)} \le C_3 (\alpha \beta + \beta^2), \\ & \left| \frac{\partial}{\partial \tilde{t}} \tilde{E} + \frac{\partial}{\partial \tilde{x}} \tilde{q}_{E_x} + \frac{\partial}{\partial \tilde{y}} \tilde{q}_{E_y} \right|_{L^\infty(\mathbb{R}^2)} \le C_4 (\alpha \beta + \beta^2), \end{aligned} $$

where $\tilde{E}, \tilde{q}_{E_x}$, and $\tilde{q}_{E_y}$ are given by:

$$ \begin{aligned} \tilde{E} &= \frac{1}{2}(\tilde{U}^2 + \tilde{V}^2 + \tilde{\eta}^2) + \frac{\beta}{2}\left(\theta^2 - \frac{1}{3}\right)(\tilde{U}\Delta\tilde{U} + \tilde{V}\Delta\tilde{V}) + \frac{\beta}{6}(\tilde{U}_{\tilde{x}} + \tilde{V}_{\tilde{y}})^2 + \frac{\alpha}{2}\tilde{\eta}(\tilde{U}^2 + \tilde{V}^2), \\ \tilde{q}_{E_x} &= \frac{\alpha}{2}(\tilde{U}^3 + \tilde{V}^2\tilde{U}) + \alpha\tilde{\eta}^2\tilde{U} + \tilde{\eta}\tilde{U} + \frac{\beta}{2}\left(\theta^2 - \frac{1}{3}\right)\tilde{\eta}\Delta\tilde{U} - \frac{\beta}{3}\tilde{U}(\tilde{U}_{\tilde{x}\tilde{t}} + \tilde{V}_{\tilde{y}\tilde{t}}), \\ \tilde{q}_{E_y} &= \frac{\alpha}{2}(\tilde{V}^3 + \tilde{U}^2\tilde{V}) + \alpha\tilde{\eta}^2\tilde{V} + \tilde{\eta}\tilde{V} + \frac{\beta}{2}\left(\theta^2 - \frac{1}{3}\right)\tilde{\eta}\Delta\tilde{V} - \frac{\beta}{3}\tilde{V}(\tilde{U}_{\tilde{x}\tilde{t}} + \tilde{V}_{\tilde{y}\tilde{t}}). \end{aligned} $$

*Proof.* We present a full justification only for the first inequality (approximate mass balance). The others follow by similar substitutions. In fact,

$$ (4.47) \qquad \frac{\partial}{\partial \tilde{t}} \{1 + \alpha \tilde{\eta}\} + \frac{\partial}{\partial \tilde{x}} \int_{\tilde{z}=0}^{1+\alpha \tilde{\eta}} \alpha P_{\tilde{x}}^{series} d\tilde{z} + \frac{\partial}{\partial \tilde{y}} \int_{\tilde{z}=0}^{1+\alpha \tilde{\eta}} \alpha P_{\tilde{y}}^{series} d\tilde{z} = O(\alpha \beta^2). $$
We replace $P_{\tilde{x}}^{series}$ and $P_{\tilde{y}}^{series}$ by expressions (4.45) and (4.46) as follows

$$ \frac{\partial}{\partial \tilde{t}} \{1 + \alpha \tilde{\eta}\} + \frac{\partial}{\partial \tilde{x}} \left\{ \alpha \tilde{u}(1 + \alpha \tilde{\eta}) - \frac{\alpha(1 + \alpha \tilde{\eta})^3}{6} \beta \Delta \tilde{u} + \frac{(1 + \alpha \tilde{\eta})^5}{120} \alpha \beta^2 (\Delta^2 \tilde{u}) \right\} \\ + \frac{\partial}{\partial \tilde{y}} \left\{ \alpha \tilde{v}(1 + \alpha \tilde{\eta}) - \frac{\alpha(1 + \alpha \tilde{\eta})^3}{6} \beta \Delta \tilde{v} + \frac{(1 + \alpha \tilde{\eta})^5}{120} \alpha \beta^2 (\Delta^2 \tilde{v}) \right\} = \mathcal{O}(\alpha \beta^2). $$

Substituting $\tilde{u}$ and $\tilde{v}$ by their expressions in terms of $\tilde{U}$ and $\tilde{V}$ given up to order $\mathcal{O}(\beta^3)$ by (4.45) and (4.46), one can write:

$$ \frac{\partial}{\partial \tilde{t}} \{1 + \alpha \tilde{\eta}\} + \frac{\partial}{\partial \tilde{x}} \left\{ \tilde{U}(\alpha + \alpha^2 \tilde{\eta}) + \frac{\alpha \beta}{2} \left( \theta^2 - \frac{1}{3} \right) \Delta \tilde{U} \right\} + \frac{\partial}{\partial \tilde{y}} \left\{ \tilde{V}(\alpha + \alpha^2 \tilde{\eta}) + \frac{\alpha \beta}{2} \left( \theta^2 - \frac{1}{3} \right) \Delta \tilde{V} \right\} = \Xi. $$

Where,

$$ \Xi = \frac{\partial}{\partial \tilde{x}} \left\{ \frac{\alpha^2 \beta}{6} (\alpha^2 \tilde{\eta}^3 + 3\alpha \tilde{\eta}^2 + 3\tilde{\eta} - 3\theta^2 \tilde{\eta}) \Delta \tilde{U} - \frac{\alpha \beta^2}{120} (\alpha^5 \tilde{\eta}^5 + 5\alpha^4 \tilde{\eta}^4 + 10\alpha^3 \tilde{\eta}^3 - 10\alpha^3 \theta^2 \tilde{\eta}^3 + 10\alpha^2 \tilde{\eta}^2 \right. \\ - 30\alpha^2 \theta^2 \tilde{\eta}^2 + 5\alpha \tilde{\eta} - 30\alpha \theta^2 \tilde{\eta} + 25\alpha \theta^4 \tilde{\eta} - 10\theta^2 + 25\theta^4 + 1) \Delta^2 \tilde{U} \right\} \\ + \frac{\partial}{\partial \tilde{y}} \left\{ \frac{\alpha^2 \beta}{6} (\alpha^2 \tilde{\eta}^3 + 3\alpha \tilde{\eta}^2 + 3\tilde{\eta} - 3\theta^2 \tilde{\eta}) \Delta \tilde{V} - \frac{\alpha \beta^2}{120} (\alpha^5 \tilde{\eta}^5 + 5\alpha^4 \tilde{\eta}^4 + 10\alpha^3 \tilde{\eta}^3 - 10\alpha^3 \theta^2 \tilde{\eta}^3 + 10\alpha^2 \tilde{\eta}^2 \right. \\ - 30\alpha^2 \theta^2 \tilde{\eta}^2 + 5\alpha \tilde{\eta} - 30\alpha \theta^2 \tilde{\eta} + 25\alpha \theta^4 \tilde{\eta} - 10\theta^2 + 25\theta^4 + 1) \Delta^2 \tilde{V} \right\} + \mathcal{O}(\alpha \beta^2). $$

Since $\tilde{U}$ and $\tilde{V}$ can be expressed in terms of $\tilde{\eta}$ and $\tilde{\psi}$, hence taking $s$ large enough such that $H^s(\mathbb{R}^2) \hookrightarrow L^\infty(\mathbb{R}^2)$, we get:

$$ \Xi = \mathcal{O}(\alpha \beta^2, \alpha^2 \beta). $$

Replacing $P_{\tilde{x}}^{series}$ and $P_{\tilde{y}}^{series}$ by expressions (4.45) and (4.46) once again in

$$ \alpha \frac{\partial}{\partial \tilde{t}} \int_{\tilde{z}=0}^{1+\alpha \tilde{\eta}} P_{\tilde{x}}^{series} d\tilde{z} + \frac{\partial}{\partial \tilde{x}} \int_{\tilde{z}=0}^{1+\alpha \tilde{\eta}} \left\{ \alpha^2 (P_{\tilde{x}}^{series})^2 + \alpha \tilde{Q} - (\tilde{z} - 1) \right\} d\tilde{z} + \alpha^2 \frac{\partial}{\partial \tilde{y}} \int_{\tilde{z}=0}^{1+\alpha \tilde{\eta}} P_{\tilde{x}}^{series} P_{\tilde{y}}^{series} d\tilde{z} = \mathcal{O}(\alpha \beta, \beta^2), $$

For the approximate momentum balance we give some details on the pressure term,

$$ \int_{0}^{1+\alpha \tilde{\eta}} \tilde{Q} d\tilde{z} = - \frac{\partial}{\partial \tilde{x}} \int_{0}^{1+\alpha \tilde{\eta}} \left\{ P_{\tilde{t}}^{series} + \frac{1}{2} \alpha \left( (P_{\tilde{x}}^{series})^2 + (P_{\tilde{y}}^{series})^2 \right) - \frac{1}{2} \frac{\alpha}{\beta} (P_{\tilde{z}}^{series})^2 \right\} d\tilde{z} $$

Using the expressions of $P^{series}$ and $f$ given by (4.6) and (3.9) respectively, and then substituting the expressions of $\tilde{u}$ and $\tilde{v}$ in terms of of $\tilde{U}$ and $\tilde{V}$ given by (4.45) and (4.46), one can write

$$ \int_{0}^{1+\alpha \tilde{\eta}} \tilde{Q} d\tilde{z} = \frac{\partial}{\partial \tilde{x}} \left\{ \tilde{\eta} + \alpha \tilde{\eta}^2 - \frac{\beta}{3} (\tilde{U}_{\tilde{x}\tilde{t}} + \tilde{V}_{\tilde{x}\tilde{t}}) + r' \right\} $$

with

$$ r' = \frac{\alpha^2 \beta}{2} \left( \tilde{\eta}^2 + \frac{\alpha}{3} \tilde{\eta}^3 \right) (\tilde{U}_{\tilde{x}\tilde{t}} + \tilde{V}_{\tilde{x}\tilde{t}}) - \frac{\beta^2}{120} (1 + \alpha \tilde{\eta})^5 (\Delta \tilde{U}_{\tilde{x}\tilde{t}} + \Delta \tilde{V}_{\tilde{x}\tilde{t}}) - \frac{\alpha \beta}{6} (1 + \alpha \tilde{\eta})^3 \tilde{U} \Delta \tilde{U} + \frac{\alpha \beta^2}{2} \left( - \frac{\theta^2}{6} (1 + \alpha \tilde{\eta})^3 \right. \\ \left. \left( (\Delta \tilde{U})^2 + (\Delta \tilde{V})^2 + \tilde{U} \Delta^2 \tilde{U} + \tilde{V} \Delta^2 \tilde{V} \right) + \alpha \theta^2 \left( \tilde{\eta}^2 + \frac{\alpha}{3} \tilde{\eta}^3 \right) (\Delta \tilde{U}_{\tilde{x}\tilde{t}} + \Delta \tilde{V}_{\tilde{x}\tilde{t}}) + \frac{(1 + \alpha \tilde{\eta})^5}{12} (\tilde{U} \Delta^2 \tilde{U} + \tilde{V} \Delta^2 \tilde{V}) \right. \\ \left. + \frac{(1 + \alpha \tilde{\eta})^3}{12} \left( (\Delta \tilde{U})^2 + (\Delta \tilde{V})^2 \right) \right) - \frac{\alpha \beta}{6} (1 + \alpha \tilde{\eta})^3 (\tilde{U}_{\tilde{x}} + \tilde{V}_{\tilde{x}})^2 + \frac{\alpha \beta^2}{6} \theta^2 \left( (\Delta \tilde{U}_{\tilde{x}})^2 + (\Delta \tilde{V}_{\tilde{x}})^2 + \tilde{U}_{\tilde{x}} \Delta \tilde{V}_{\tilde{x}} + \tilde{V}_{\tilde{x}} \Delta \tilde{U}_{\tilde{x}} \right) \\ + \frac{\alpha \beta^2}{30} (1 + \alpha \tilde{\eta})^5 (\tilde{U}_{\tilde{x}} + \tilde{V}_{\tilde{x}}) (\Delta \tilde{U}_{\tilde{x}} + \Delta \tilde{V}_{\tilde{x}}) - \frac{\beta^2}{6} \theta^2 (\Delta \tilde{U}_{\tilde{x}\tilde{t}} + \Delta \tilde{V}_{\tilde{x}\tilde{t}}) + \mathcal{O}(\beta^3). $$

The other terms are treated similarly as the terms in the approximate mass conservation equation. $\square$
5. NUMERICAL VALIDATION

In this final section, we present a computational study to gain insight into the approximate conservation
of mass (4.35), momentum (4.38)–(4.39), and energy (4.41) in a practical scenario. We focus on the circular
expansion of water waves originating from an initial localized source.

For this purpose, we numerically solve the Boussinesq system (3.17a)–(3.17c) with $\theta^2 = 9/11$, using periodic boundary conditions within the domain $D = [-20, 20] \times [-20, 20]$. Boussinesq systems corresponding to different values of $\theta$ are asymptotically equivalent to the present system and exhibit similar behavior. Our numerical method employs the standard pseudo-spectral approach in combination with the classical four-stage, fourth-order Runge-Kutta method.

We initiate the simulation with the following initial conditions: $\eta(x, y, 0) = \exp(-(x^2+y^2)/5)$, $u(x, y, 0) = v(x, y, 0) = 0$, for various values of $\alpha = \beta$ as detailed in Table 1. We numerically integrate up to time $T = 10$ using a time step of $\Delta t = 10^{-4}$ and a spatial step of $\Delta x = 0.1$. Throughout the computation, we record values of the discrete quantities (4.35), (4.38)–(4.39), and (4.41).

To discretize the temporal derivatives in these quantities, we employ forward finite differences, while the
spatial derivatives are computed using the fast Fourier transform.

<figure><img src="image_4.png" /><figcaption>FIGURE 2. Generation of expanding water waves from a heap of water. The waves are propagating in a circular pattern with growing radius.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 3. Evolution of the amplitude of the leading expanding wave. The leading wave in the numerical approximation is shown in blue.</figcaption></figure>

The specific initial condition resembles a mound of water, often observed when an object falls into the
water or when a tsunami is triggered by changes in the bottom topography (for more details, see [19]).
<table><caption>TABLE 1. Residuals of mass, momentum and energy conservation laws evaluated using the numerical approximation of a solution of the Boussinesq system.</caption><thead><tr><th>α = β</th><th>Mass</th><th>Momentum</th><th>Energy</th></tr></thead><tbody><tr><td>0.05</td><td>2.21 × 10<sup>-5</sup></td><td>4.16 × 10<sup>-3</sup></td><td>3.76 × 10<sup>-4</sup></td></tr><tr><td>0.10</td><td>1.57 × 10<sup>-4</sup></td><td>8.22 × 10<sup>-3</sup></td><td>1.35 × 10<sup>-3</sup></td></tr><tr><td>0.15</td><td>4.99 × 10<sup>-4</sup></td><td>1.21 × 10<sup>-2</sup></td><td>2.84 × 10<sup>-3</sup></td></tr><tr><td>0.20</td><td>1.12 × 10<sup>-3</sup></td><td>1.59 × 10<sup>-2</sup></td><td>4.76 × 10<sup>-3</sup></td></tr><tr><td>0.25</td><td>2.08 × 10<sup>-3</sup></td><td>1.96 × 10<sup>-2</sup></td><td>7.05 × 10<sup>-3</sup></td></tr><tr><td>0.30</td><td>3.43 × 10<sup>-3</sup></td><td>2.33 × 10<sup>-2</sup></td><td>9.67 × 10<sup>-3</sup></td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>FIGURE 4. Time development of the residuals (absolute errors) of mass, momentum and energy balance laws using the numerical approximation of a solution of the Boussinesq system.</figcaption></figure>

Figure 2 illustrates the initial condition and its evolution at $T = 10$. The initial condition transforms into a series of expanding waves traveling at varying speeds due to dispersion. Over time, the amplitudes of these waves decrease as a result of dispersion's influence. This reduction in amplitude follows a rate proportional to $t^{-1}$ due to the radial spreading of the two-dimensional waves [21]. Figure 3 presents the amplitude of the leading expanding wave, and for a detailed derivation of the decrease rate, we direct interested readers to [21].

Subsequently, for the specific initial data presented in Figure 2, the computed residuals of the mass, momentum, and energy balance laws initially increase but decay over time. The maximum value attained by each conservation law is listed in Table 1. It's worth noting that these maximum values were not reached again, and the residuals appeared to decrease on average. Figure 4 displays the absolute error patterns for the case where $α = β = 0.3$; patterns for other cases are quite similar and thus omitted.

A plausible explanation for the decaying of the residuals is that the expanding waves tend to approach linear waves, as explained in the previous section. These numerical results validate the orders of accuracy of the balance laws, particularly in the common scenario of dispersive waves generated from a general initial condition.
ACKNOWLEDGMENTS

This research was supported by the Research Council of Norway under grant numbers 213474/F20 and 239033/F20. SI would like to thank the Department of Mathematics at the University of Bergen, Norway, and especially Henrik Kalisch for their kind hospitality when work on this paper was begun.

REFERENCES

[1] A. Ali and H. Kalisch. Mechanical balance laws for Boussinesq models of surface water waves. *J Nonlinear Sci*, 22:371–398, 2012.

[2] A. Ali and H. Kalisch. On the formulation of mass, momentum and energy conservation in the KdV equation. *Acta Appl Math*, 133:113–131, 2014.

[3] J. Bona, M. Chen, and J.-C. Saut. Boussinesq equations and other systems for small-amplitude long waves in nonlinear dispersive media. i: Derivation and linear theory. *J.Nonlin. Sci.*, 12(4), 2002.

[4] J. Bona, T. Colin, and D. Lannes. Long wave approximations for water waves. *Arch. Rat. Mech. Anal.*, 178:373–410, 2005.

[5] J. Boussinesq. Théorie des ondes et des remous qui se propagent le long d'un canal rectangulaire horizontal, en commu- niquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond. *J. Math. Pures. Appl.*, 17:55-108, 1872.

[6] V. Dougalis and D. Mitsotakis. Theory and numerical analysis of boussinesq systems: A review. In N. Kampanis, V. Dou- galis, and J. Ekaterinaris, editors, *Effective Computational Methods in Wave Propagation*, pages 63–110. CRC Press, 2008.

[7] V. Dougalis, D. Mitsotakis, and J.-C. Saut. On some Boussinesq systems in two space dimensions: Theory and numerical analysis. *ESAIM: Math. Model. Num. Anal.*, 41:825–854, 2007.

[8] V. Dougalis, D. Mitsotakis, and J.-C. Saut. On initial-boundary value problems for a Boussinesq system of BBM-BBM type in a plane domain. *Discrete Contin. Dyn. Syst*, 23:1191–1204, 2009.

[9] V. Dougalis, D. Mitsotakis, and J.-C. Saut. Boussinesq systems of Bona-Smith type on plane domains: theory and numerical analysis. *J. Sci. Comp.*, 44(2):109–135, 2010.

[10] A.E. Green and P.M. Naghdi. A derivation of equations for wave propagation in water of variable depth. *J. Fluid Mech.*, 78:237–246, 1976.

[11] S. Israwi and H. Kalisch. Approximate conservation laws in the KdV equation. *Physics Letters A*, 383(9):854-858, 2019.

[12] S. Israwi and H. Kalisch. A mathematical justification of the momentum density function associated to the KdV equation. *Comptes Rendus. Mathématique*, 359(1):39-45, 2021.

[13] S Israwi, H. Kalisch, Th. Katsaounis, and D. Mitsotakis. A regularized shallow-water waves system with slip-wall boundary conditions in a basin: theory and numerical analysis. *Nonlinearity*, 35:750–786, 01 2022.

[14] H. Kalisch, Z. Khorsand, and D. Mitsotakis. Mechanical balance laws for fully nonlinear and weakly dispersive water waves. *Physica D: Nonlinear Phenomena*, 333:243–253, 2016.

[15] B. Khorbatly, R. Lteif, S. Israwi, and St. Gerbi. Mathematical modeling and numerical analysis for the higher order Boussinesq system. *ESAIM: M2AN*, 56(2):593–615, 2022.

[16] B. Khorbatly, I. Zaiter, and S. Israwi. Derivation and well-posedness of the extended Green-Naghdi equations for flat bottoms with surface tension. *J. Math. Phys.*, 59(7):071501, 20, 2018.

[17] D. Lannes. *The water waves problem: mathematical analysis and asymptotics*, volume 188. Americal Mathematical Society, Providence, Rhode Island, 2013.

[18] D. Lannes and P. Bonneton. Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation. *Physics of Fluids*, 21:016601, 2009.

[19] D. Mitsotakis. Boussinesq systems in two space dimensions over a variable bottom for the generation and propagation of tsunami waves. *Mathematics and Computers in Simulation*, 80(4):860–873, 2009.

[20] F. Serre. Contribution à l'étude des écoulements permanents et variables dans les canaux. *La Houille blanche*, 8:374-872, 1953.

[21] G.B. Whitham. *Linear and nonlinear waves*, volume 42. John Wiley & Sons, New York, 1974.
**C. E. Hassanieh:** DEPARTMENT OF MATHEMATICS, FACULTY OF SCIENCES 1, LEBANESE UNIVERSITY, BEIRUT, LEBANON
& SORBONNE UNIVERSITÉ, 4 PLACE JUSSIEU, 75005 PARIS, FRANCE & INRIA PARIS

*Email address: chourouk.el-hassanieh@inria.fr*

**S. Israwi:** DEPARTMENT OF MATHEMATICS, FACULTY OF SCIENCES 1, LEBANESE UNIVERSITY, BEIRUT, LEBANON
& ARTS, SCIENCES, AND TECHNOLOGY UNIVERSITY IN LEBANON, CRAMS CENTER OF RESEARCH IN APPLIED MATHEMATICS
AND STATISTICS, FACULTY OF SCIENCES, BEIRUT, LEBANON

*Email address: s_israwi83@hotmail.com*

**H. Kalisch:** DEPARTMENT OF MATHEMATICS, UNIVERSITY OF BERGEN, POSTBOX 7800, 5020 BERGEN, NORWAY

*Email address: henrik.kalisch@uib.no*

**D. Mitsotakis:** VICTORIA UNIVERSITY OF WELLINGTON, SCHOOL OF MATHEMATICS AND STATISTICS, PO BOX 600, WELLINGTON 6140, NEW ZEALAND

*Email address: dimitrios.mitsotakis@vuw.ac.nz*

**A. Senthilkumar:** DEPARTMENT OF MATHEMATICS, UNIVERSITY OF BERGEN, POSTBOX 7800, 5020 BERGEN, NORWAY

*Email address: amutha.senthilkumar@math.uib.no*