The IsobeKakinuma model
Higher order models for surface gravity waves
The socalled IsobeKakinuma model is a hierarchy of partial differential equations modelling the propagation of surface gravity waves (or water waves).
They have been introduced by M. Isobe^{1} (and numerically studied by T. Kakinuma^{2}) using J. C. Luke’s Lagrangian formulation^{3} of the (incompressible, homogeneous) Bernoulli equations with freesurface. Alternatively^{4}, one can derive the equations by using a Galerkin dimension reduction method on a variational formulation of the underlying Laplace problem. With an appropriate function system basis for the vertical profile of the velocity potential, one obtains the following set of equations.
$$ \begin{cases} \sum_{j=0}^N h^{p_j}(\partial_t\phi_j )+g\zeta + \frac{1}2\Big( \big \sum_{j=0}^N h^{p_j} (\nabla\phi_j)\big^2 \\ \qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad +\big( \sum_{j=0}^N p_j h^{p_j1} \phi_j\big)^2 \Big) = 0,\\ h^{p_i}\partial_t\zeta+\sum_{j=0}^N\nabla\cdot\big( \tfrac{h^{p_i+p_j+1}}{p_i+p_j+1} \nabla\phi_j \big)  \sum_{j=0}^N\tfrac{p_ip_j}{p_i+p_j1}h^{p_i+p_j1}\phi_j = 0 \\ \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \qquad \qquad \forall i\in{0,1,\dots,N}. \end{cases} $$ Above, $h$ denotes the height of the fluid layer at a prescribed horizontal location, $x$, and time, $t$. The graph of $\zeta$ describes the free surface, so that $h=d+\zeta$ where the constant $d$ is the depth of the layer at rest. The constant $g$ denotes the vertical gravity acceleration. The $\phi_j$ are artificial unknowns related to the velocity potential. Finally, the $p_j$ are integers and $N$ is the order of the model, roughly speaking measuring the size and hence the complexity of the model. Although it has not yet been proved, it expected that choosing $p_j=2j$ for $j\in{0,1,\dots,N}$, the solution to the IsobeKakinuma model approximates the corresponding solution to the free surface Bernoulli equations with an arbitrary accuracy as $N$ go to infinity. In that way, the IsobeKakinuma model provides a strategy for numerically computing the free surface Bernoulli equations without having to deal with a threedimensional unsteady domain.
One can remark that the IsobeKakinuma model has a peculiar structure. There are $N+1$ evolution equations for the surface deformation, and only one equation prescribing the evolution of (a combination of) the $\phi_j$. However the number of equations coincide with the number of unknowns, and this over/underdetermination can be tackled. In fact, by combining the last $N+1$ equations, we obtain a set of $N$ compatibility conditions, and we can view the IsobeKakinuma model as a set of two evolution equations coupled with a set of constraints, which allows to reconstruct uniquely (or up to a harmless constant) the full set of unknowns. This is reminiscent of the ZakharovCraigSulem formulation^{5} of the freesurface Bernoulli equations, where the equations are written equivalently as two evolution equations for the surface deformation variable and the trace of the velocity potential at the surface, the velocity potential in the bulk of the domain being reconstructed by the aforementioned Laplace problem. One key property of this formulation is that it brings to lights the canonical Hamiltonian structure of the equations. Together with T. Iguchi, we have proved^{6} that the IsobeKakinuma model enjoys an analogous canonical Hamiltonian structure.
However this Hamiltonian structure is not the one which is the most suitable for the study of the initialvalue problem for instance. By using a clever tweak on the energy method on the original formulation of the equations, T. Iguchi and collaborators ^{7} have been able to rigorously justify the IsobeKakinuma model as a higher order approximate model for the free surface Bernoulli equations in the shallow water regime. Roughly speaking, the result is as follows:
 For any sufficiently regular initial data satisfying the aforementioned compatibility conditions and the noncavitation assumption^{8}, there exists a unique solution to the IsobeKakinuma model on the relevant timescale
 This solution is at a distance $O(\mu^{2N+1})$ to the corresponding solution to the free surface Bernoulli equations, where $\mu$ is the dimensionless shallowness parameter measuring the ratio of a typical vertical length to a typical horizontal length.
It should be noted that these result apply to the dimensionless versions of the equations, and are uniform with respect to $\mu\in (0,1]$. This is a key point for the justification of the model in the shallowwater regime, and brings specific difficulties as the limit $\mu\to 0$ is singular. However, the result is not uniform with respect to $N$: in particular the level of regularity (measured through the index of Sobolev norms) that we impose on the initial data to secure the $O(\mu^{2N+1})$ precision is growing with $N$.
I would like to conclude with a brief description of the strategy which allows to obtain such a result. The key point consists in obtaing stability estimates for a linearized version of the IsobeKakinuma equations. It can be checked that the linearized IsobeKakinuma equations (after slight manpulations) reads as follows. $$ \mathcal{A}\partial_t U + \mathcal{L} U = F .$$ Here, $U$ is the set of unknowns, $\mathcal{A}$ and $\mathcal{L}$ are linear operators, and $F$ is a harmless remainder term. In the standard energy method for (say) symmetrizable hyperbolic quasilinear systems, the operators $\mathcal{A}$ is bounded, selfadjoint (for the $L^2$ innerproduct) positive definite, and the operator $\mathcal{L}$ can be decomposed as the sum of a skewHermitian (for the $L^2$ innerproduct) operator and a bounded operator. In that case, by testing the equation against $U$ (that is, assuming a solution exists, considering the identity obtained by taking the $L^2$ innerproduct with $U$ of both sides of the equation), and using the (skew)symmetry properties and CauchySchwarz inequality, we infer the differential inequality $$ \frac12 \frac{{\rm d}}{{\rm d} t}\big(\mathcal{A} U,U \big)_{L^2} \leq C_0  F_{L^2}U_{L^2}+C_1 U_{L^2}^2$$ with $C_0$ and $C_1$ some constants. Using that $\big(\mathcal{A} U,U \big)_{L^2}\approx U_{L^2}^2$ and Gronwall’s Lemma, we infer an a priori control of the quantity $\big(\mathcal{A} U,U )_{L^2}$, and hence of $U_{L^2}$. This control is in terms crucial to obtain stability estimates concerning the nonlinear and linearized equations.
In our case the situation is different as it is $\mathcal{A}$ which is skewHermitian and $\mathcal{L} $ can be decomposed as $\mathcal{L}=\mathcal{A}{\bf u}\cdot\nabla+\mathcal{L}_0$ where $\mathcal{L}_0$ is a bounded, selfadjoint positive definite linear operator. We can hence rewrite the linearized system of equations as $$ \mathcal{A}(\partial_t U +{\bf u}\cdot\nabla U)+ \mathcal{L}_0 U = F .$$ Then the trick consists in testing the equations against $\partial_t U +{\bf u}\cdot\nabla U$. The contribution of $\mathcal{A}$ vanishes due to the skewHermitian property, while the contribution of $\mathcal{L}_0$ allows to control (proceeding roughly speaking as above) the quantity $(\mathcal{L}_0U,U)_{L^2}\approx U_{L^2}^2$. As the standard energy method, this analysis is quite robust, and in particular behaves well with the presence of singular terms with respect to the shoallowness parameter $\mu$. Together with a fine analysis of the system of compatibility constraints (which again relies on symmetry proprties, but this time in view of elliptic estimates), this provides the key ingredient for the the results of T. Iguchi and collaborators stated above.

M. Isobe, A proposal on a nonlinear gentle slope wave equation, Proceedings of Coastal Engineering, Japan Society of Civil Engineers, 41 (1994)
M. Isobe, Timedependent mildslope equations for random waves, Proceedings of 24th International Conference on Coastal Engineering, ASCE (1994) ↩︎ 
T. Kakinuma, [title in Japanese], Proceedings of Coastal Engineering, Japan Society of Civil Engineers, 47 (2000)
T. Kakinuma, A set of fully nonlinear equations for surface and internal gravity waves, Coastal Engineering V: Computer Modelling of Seas and Coastal Regions (2001)
T. Kakinuma, A nonlinear numerical model for surface and internal waves shoaling on a permeable beach, Coastal engineering VI: Computer Modelling and Experimental Measurements of Seas and Coastal Regions(2003) ↩︎ 
J. C. Luke, A variational principle for a fluid with a free surface, J. Fluid Mech., 27 (1967) ↩︎

V. Duchêne, Many Models for Water Waves. Open Math Notes 202109.111309 ↩︎

V. E. Zakharov, Stability of periodic waves of finite amplitude on the surface of a deep fluid, J. Appl. Mech. Tech. Phys., 9 (1968)
W. Craig and C. Sulem, Numerical simulation of gravity waves, J. Comput. Phys., 108 (1993) ↩︎ 
V. Duchêne and T. Iguchi, A Hamiltonian structure of the Isobe–Kakinuma model for water waves, Water Waves, 3 (2021) ↩︎

Y. Murakami and T. Iguchi, Solvability of the initial value problem to a model system for water waves, Kodai Math. J., 38 (2015)
R. Nemoto and T. Iguchi, Solvability of the initial value problem to the Isobe–Kakinuma model for water waves, J. Math. Fluid Mech., 20 (2018)
T. Iguchi, Isobe–Kakinuma model for water waves as a higher order shallow water approximation, J. Differential Equations, 265 (2018)
T. Iguchi, A mathematical justification of the Isobe–Kakinuma model for water waves with and without bottom topography, J. Math. Fluid Mech., 20 (2018) ↩︎ 
the height of water must be everywhere bounded from below by a positive constant. ↩︎