The Isobe-Kakinuma model

The so-called Isobe-Kakinuma 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¹ (and numerically studied by T. Kakinuma²) using J. C. Luke’s Lagrangian formulation³ of the (incompressible, homogeneous) Bernoulli equations with free-surface. Alternatively⁴, 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_j-1} \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_j-1}h^{p_i+p_j-1}\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 Isobe-Kakinuma model approximates the corresponding solution to the free surface Bernoulli equations with an arbitrary accuracy as $N$ go to infinity. In that way, the Isobe-Kakinuma model provides a strategy for numerically computing the free surface Bernoulli equations without having to deal with a three-dimensional unsteady domain.

One can remark that the Isobe-Kakinuma 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/under-determination 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 Isobe-Kakinuma 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 Zakharov-Craig-Sulem formulation⁵ of the free-surface 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⁶ that the Isobe-Kakinuma 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 initial-value problem for instance. By using a clever tweak on the energy method on the original formulation of the equations, T. Iguchi and collaborators ⁷ have been able to rigorously justify the Isobe-Kakinuma 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:

1. For any sufficiently regular initial data satisfying the aforementioned compatibility conditions and the non-cavitation assumption⁸, there exists a unique solution to the Isobe-Kakinuma model on the relevant timescale
2. 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 shallow-water 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 Isobe-Kakinuma equations. It can be checked that the linearized Isobe-Kakinuma 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, self-adjoint (for the $L^2$ inner-product) positive definite, and the operator $\mathcal{L}$ can be decomposed as the sum of a skew-Hermitian (for the $L^2$ inner-product) 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$ inner-product with $U$ of both sides of the equation), and using the (skew)-symmetry properties and Cauchy-Schwarz 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 skew-Hermitian 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, self-adjoint 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 skew-Hermitian 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.