A linear and nonlinear analysis of the shallow water equations and its impact on boundary conditions

TL;DR

This paper compares linear and nonlinear analyses of the shallow water equations, revealing differences in boundary condition requirements and emphasizing the importance of nonlinear energy methods for accurate modeling.

Contribution

It introduces a nonlinear energy analysis approach that refines boundary condition estimates for the shallow water equations, highlighting differences from linear analysis.

Findings

Nonlinear analysis alters the minimal number of boundary conditions needed.

Flow magnitude does not affect boundary condition requirements in the nonlinear case.

Energy and entropy analyses are consistent but provide different insights.

Abstract

We derive boundary conditions and estimates based on the energy and entropy analysis of systems of the nonlinear shallow water equations in two spatial dimensions. It is shown that the energy method provides more details, but is fully consistent with the entropy analysis. The details brought forward by the nonlinear energy analysis allow us to pinpoint where the difference between the linear and nonlinear analysis originate. We find that the result from the linear analysis does not necessarily hold in the nonlinear case. The nonlinear analysis leads in general to a different minimal number of boundary conditions compared with the linear analysis. In particular, and contrary to the linear case, the magnitude of the flow does not influence the number of required boundary conditions.

Equations

$$h_t+(hu{)}_x+(hv{)}_y$$

$$(hu{)}_t+{\left(h{u}^2+\frac{g}{2}{h}^2\right)}_x+(huv{)}_y-fhv$$

$$(hv{)}_t+(huv{)}_x+{\left(h{v}^2+\frac{g}{2}{h}^2\right)}_y+fhu$$

$${\varphi }_t+u{\varphi }_x+v{\varphi }_y+\varphi (u_x+v_y)$$

$$u_t+u{u}_x+v{u}_y+{\varphi }_x-fv$$

$$v_t+u{v}_x+v{v}_y+{\varphi }_y+fu$$

$${\mathbf{q}}_t+\mathcal{A}{\mathbf{q}}_x+\mathcal{B}{\mathbf{q}}_y+\mathcal{C}\mathbf{q}=0,$$

$$\mathcal{A}=\left[\begin{array}{ccc}u& \varphi & 0\\ 1& u& 0\\ 0& 0& u\end{array}\right],\mathcal{B}=\left[\begin{array}{ccc}v& 0& \varphi \\ 0& v& 0\\ 1& 0& v\end{array}\right],\mathcal{C}=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& f\\ 0& -f& 0\end{array}\right].$$

$$ϵ=\frac{\varphi }{2g}\left(u^2+v^2\right)+\frac{{\varphi }^2}{2g}.$$

$$f^{ϵ}=\frac{\varphi u}{2g}\left(u^2+v^2\right)+\frac{{\varphi }^{2}u}{g},\text{and}{g}^{ϵ}=\frac{\varphi v}{2g}\left(u^2+v^2\right)+\frac{{\varphi }^{2}v}{g},$$

$${ϵ}_t+f_x^{ϵ}+g_y^{ϵ}=0.$$

$$s_1=\frac{\partial ϵ}{\partial h},s_2=\frac{\partial ϵ}{\partial (hu)},s_3=\frac{\partial ϵ}{\partial (hv)}.$$

$${\mathbf{q}}_t+\mathsf{D}\mathbf{q}$$

$$\mathsf{L}\mathbf{q}$$

$$\mathbf{q}$$

$$⟨\mathbf{u},\mathbf{v}⟩=\underset{\Omega }{\int }{\mathbf{u}}^T\mathcal{P}\mathbf{v}\mathrm{d}\Omega ,\Vert \mathbf{u}{\Vert }^2=⟨\mathbf{u},\mathbf{u}⟩,$$

$$⟨\mathbf{u},\mathsf{D}\mathbf{u}⟩\ge \lambda \Vert \mathbf{u}{\Vert }^2$$

$$\Vert \mathbf{q}{\Vert }_t^2=2⟨\mathbf{q},{\mathbf{q}}_t⟩=-2⟨\mathbf{q},\mathsf{D}\mathbf{q}⟩\le -2\lambda \Vert \mathbf{q}{\Vert }^2\Rightarrow \Vert \mathbf{q}(\cdot ,T)\Vert \le e^{-\lambda T}\Vert f(\cdot )\Vert .$$

$$\mathcal{S}=\frac{1}{\sqrt{2g}}\left[\begin{array}{ccc}1& 0& 0\\ 0& \sqrt{\varphi }& 0\\ 0& 0& \sqrt{\varphi }\end{array}\right].$$

$$\mathcal{S}\mathcal{A}{\mathcal{S}}^{-1}=\left[\begin{array}{ccc}u& \sqrt{\varphi }& 0\\ \sqrt{\varphi }& u& 0\\ 0& 0& u\end{array}\right],\mathcal{SB}{\mathcal{S}}^{-1}=\left[\begin{array}{ccc}v& 0& \sqrt{\varphi }\\ 0& v& 0\\ \sqrt{\varphi }& 0& v\end{array}\right].$$

$${\left(\mathcal{S}\mathbf{q}\right)}^T\mathcal{S}\mathbf{q}={\mathbf{q}}^T\mathcal{P}\mathbf{q}=\frac{\varphi u^2+\varphi v^2+{\varphi }^2}{2g},$$

$${\mathbf{q}}^T\mathcal{P}{\mathbf{q}}_t+{\mathbf{q}}^T\mathcal{P}\mathcal{A}{\mathbf{q}}_x+{\mathbf{q}}^T\mathcal{P}\mathcal{B}{\mathbf{q}}_y+{\mathbf{q}}^T\mathcal{P}\mathcal{C}\mathbf{q}=0.$$

$$\mathcal{P}\mathcal{A}=\frac{1}{2g}\left[\begin{array}{ccc}u& \varphi & 0\\ \varphi & \varphi u& 0\\ 0& 0& \varphi u\end{array}\right],\mathcal{P}\mathcal{B}=\frac{1}{2g}\left[\begin{array}{ccc}v& 0& \varphi \\ 0& \varphi v& 0\\ \varphi & 0& \varphi v\end{array}\right].$$

$${\mathbf{q}}^T\mathcal{P}{\mathbf{q}}_t$$

$${\mathbf{q}}^T\mathcal{P}\mathcal{A}{\mathbf{q}}_x$$

$${\mathbf{q}}^T\mathcal{P}\mathcal{B}{\mathbf{q}}_y$$

$${\left({\mathbf{q}}^T\mathcal{P}\mathbf{q}\right)}_t+{\left({\mathbf{q}}^T\mathcal{P}\mathcal{A}\mathbf{q}\right)}_x+{\left({\mathbf{q}}^T\mathcal{P}\mathcal{B}\mathbf{q}\right)}_y-{\mathbf{q}}^T\left({\mathcal{P}}_t+{\left(\mathcal{P}\mathcal{A}\right)}_x+{\left(\mathcal{P}\mathcal{B}\right)}_y\right)\mathbf{q}=0.$$

$${\mathcal{P}}_t=\frac{1}{2g}\left[\begin{array}{ccc}0& 0& 0\\ 0& {\varphi }_t& 0\\ 0& 0& {\varphi }_t\end{array}\right],{\left(\mathcal{P}\mathcal{A}\right)}_x$$

$${\left(\mathcal{P}\mathcal{B}\right)}_y$$

$${\mathcal{P}}_t+{\left(\mathcal{P}\mathcal{A}\right)}_x+{\left(\mathcal{P}\mathcal{B}\right)}_y=\frac{1}{2g}\left[\begin{array}{ccc}{u}_x+v_y& {\varphi }_x& {\varphi }_y\\ {\varphi }_x& 0& 0\\ {\varphi }_y& 0& 0\end{array}\right].$$

$${\mathbf{q}}^T\left({\mathcal{P}}_t+{\left(\mathcal{P}\mathcal{A}\right)}_x+{\left(\mathcal{P}\mathcal{B}\right)}_y\right)\mathbf{q}$$

$$=\frac{1}{2g}\left[{\left({\varphi }^{2}u\right)}_x+{\left({\varphi }^{2}v\right)}_y\right].$$

$${\left({\mathbf{q}}^T\mathcal{P}\mathbf{q}\right)}_t+{\left({\mathbf{q}}^T\mathcal{P}\mathcal{A}\mathbf{q}\right)}_x+{\left({\mathbf{q}}^T\mathcal{P}\mathcal{B}\mathbf{q}\right)}_y-\frac{1}{2g}\left[{\left({\varphi }^{2}u\right)}_x+{\left({\varphi }^{2}v\right)}_y\right]=0.$$

$$\frac{d}{dt}\Vert \mathbf{q}{\Vert }_{\mathcal{P}}^2+\frac{1}{2g}\underset{\partial \Omega }{\oint }\left({\mathbf{q}}^T\tilde{\mathcal{A}}\mathbf{q}-\frac{{\varphi }^2}{2g}{u}_n\right)\mathrm{d}S=0,$$

$$\frac{1}{2g}\tilde{\mathcal{A}}=\mathcal{P}\left(\mathcal{A}{n}_x+\mathcal{B}{n}_y\right)=\frac{1}{2g}\left[\begin{array}{ccc}{u}_n& n_x\varphi & n_y\varphi \\ n_x\varphi & \varphi u_n& 0\\ n_y\varphi & 0& \varphi u_n\end{array}\right].$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Shallow water models provide a system of nonlinear hyperbolic equations that are useful to model fluid flows in lakes, rivers, oceans or near coastlines, e.g [1, 2, 3]. Specific information regarding the exact physical case is generally provided by boundary and initial conditions. Proper treatment of the boundary conditions is also an essential aspect of bounding the solution in the continuous analysis. This continuous analysis, in turn, is used as a guide to develop a stable numerical approximation for the shallow water model. For large scale problems, like those found in tsunami modelling on a regional oceanic area, the lateral boundaries of the domain are not physical boundaries. Therefore, artificial “open”-type boundary conditions must be employed while maintaining proper solution estimates.

Due to the nonlinear nature of the shallow water equations (SWEs) the development of a bound of the solution becomes more subtle. In the literature one finds two general avenues in order to obtain estimates of solutions to nonlinear hyperbolic systems of conservation laws:

The energy method advocated by Kreiss [4, 5], Gustafsson [6, 7] and Oliger and Sundström [8]. 2. 2.

The entropy stability theory used by Tadmor [9, 10, 11] and originating in works by Godunov [12], Vol’pert [13], Kružkov [14], Dafermos [15], Lax [16], and Harten [17].

Traditionally, the energy method has been applied to linearized versions of systems of hyperbolic equations in order to develop boundary treatments that lead to energy estimates. In practice, these boundary conditions are needed to develop energy stable numerical approximations that weakly impose boundary information, e.g., through penalty terms [18, 19, 20, 21] or numerical flux functions [22, 23, 24, 25]. The entropy method has often been applied to nonlinear hyperbolic systems on domains with periodic boundary conditions (or infinite domains) in order to obtain entropy conservation. This makes most investigations using entropy analysis similar to the classical von Neumann analysis in the sense that boundaries are often ignored [26, 27, 28, 29]. Entropy conservation is typically an intermediate step and dissipation is added to obtain entropy stability [11]. The use of entropy stability theory to derive boundary conditions has been considered previously, see e.g.[30, 31, 22, 32, 33, 34, 35], but mostly focused on solid walls. The main character of these two approaches are (roughly speaking): The energy analysis provides boundary conditions, but is difficult to apply in the nonlinear case. The entropy analysis is straightforward to apply for nonlinear equations, but its use for development of boundary conditions is limited.

In this work, we will derive stable nonlinear boundary treatments for the SWEs in two spatial dimensions. For this, we will examine the (physical) total energy of the system. This energy is an auxiliary conserved quantity, not explicitly built into the SWEs that can be used to create a nonlinear stability estimate [36, 37, 24] as well as build numerical approximations that respect the evolution of the total energy [38, 36, 39]. An interesting quality of the SWEs and its nonlinear estimate is that the total energy also acts as a mathematical entropy function and fits into the entropy analysis framework [11].

Therefore, the total energy and the analysis of it for the SWEs act as a bridge between the classical energy method and the entropy analysis. So, a second contribution of this work is to examine the connection between the energy and entropy analysis, exemplifying their similarities and differences. We contrast the estimates from the nonlinear analysis to the ones from linear energy analysis and discuss their differences. In particular, we highlight that the linear analysis, in general, indicate a different number of boundary conditions than the nonlinear analysis. We also demonstrate that the nonlinear energy analysis provide additional details compared to the more easily applied entropy analysis, and makes it possible to pinpoint where the difference between the linear and nonlinear analysis reside.

We apply the energy method [7, 4, 5] and derive a bound for the total energy, which for the SWEs is a particular scaled version of the $L^{2}$ norm of the solution, based on a space and time-dependent weight matrix. It will be shown that the energy method is consistent with the mathematical entropy analysis, but also that it provides additional information, with respect to the boundary treatment. Investigations into energy stable boundary conditions for the linearized SWEs are many [40, 41, 42, 43, 8, 44] and leads to well-posed initial boundary value problems [42, 8, 44]. These linear boundary conditions are often applied in the nonlinear case [42, 45]. However, as we will show, the linear boundary treatments do not, in general, provide the correct minimal number of boundary conditions for the nonlinear shallow water problem.

The paper is organized as follows: Various forms of the SWEs are presented in Section 2. A brief review of energy estimates as well as linear or nonlinear boundedness and well-posedness is provided in Section 3. An estimate of the total energy and entropy for the shallow water system is provided in Section 4. In Section 5, we provide details, analysis, and discussion of the general open boundary conditions for the two-dimensional nonlinear SWEs in inflow and outflow regimes. In particular, we discuss similarities and differences between the results from the linear analysis and the nonlinear one. Concluding remarks are drawn in Section 6.

2 Shallow water equations

The SWEs over a flat bottom topography in conservative form [46] are

[TABLE]

which includes the continuity and momentum equations. Here $h(x,y,t)$ is the water height, $u(x,y,t)$ and $v(x,y,t)$ are the fluid velocities in the $x$- and $y$-directions, and $g$ is the gravitational constant. The system of equations (2.1) are derived under the physical requirement that the water height is strictly positive, i.e., $h>0$ [44, 46]. Additionally, we include the influence of Coriolis forces with the parameter $f$ which, for convenience, is assumed to be a constant. In practical applications $f$ is typically a function of latitude [44], which would not affect the subsequent analysis in this work.

In order to apply the energy method it is convenient to work with the equivalent non-conservative form of the governing equations [47, 19, 20, 8] which is

[TABLE]

In (2.2), we formulate the equations in terms of the geopotential $\phi=gh$ to simplify the analysis [8]. Note, that $\phi>0$ according to the physical and mathematical requirements of the problem. Next, we write (2.2) compactly in matrix-vector form by introducing the solution vector $\mathbf{q}=(\phi\,,\,u\,,\,v)^{T}$ as

[TABLE]

where

[TABLE]

The specific total energy is the sum of the kinetic and potential energy [36]

[TABLE]

It is an auxiliary conserved quantity for smooth solutions of (2.1) or (2.2). The specific total energy (2.3) has associated flux functions [38]

[TABLE]

that yield the specific total energy conservation law

[TABLE]

In the entropy analysis, the specific total energy in (2.3) acts as a mathematical entropy function for the SWEs [38]. As such, it is possible to define a new set of entropy variables $\mathbf{s}=(s_{1},\,s_{2},\,s_{3})^{T}$ where

[TABLE]

Multiplying the conservative form of the SWEs (2.1) from the left with $\mathbf{s}^{T}$ yields the auxiliary conservation law of the entropy function (in this case the specific total energy) (2.5). This contraction into entropy space involves the chain rule and certain compatibility conditions between the conservative fluxes from (2.1) and the entropy fluxes from (2.4) [38, 11].

Remark 2.1*.*

The Coriolis force is not present in the specific total energy conservation law (2.5). This agrees with the underlying physics of the problem because the Coriolis terms do not perform work on the fluid [44].

Remark 2.2*.*

The forthcoming nonlinear energy analysis considers the governing equations in non-conservative advective form (2.2). It is also possible to apply the nonlinear energy method to the conservative form (2.1) although it more unwieldy in terms of algebraic manipulations required to obtain an energy statement at the boundary. Further, we will demonstrate in Section 4 that the nonlinear energy analysis of the advective form is equivalent to the entropy analysis of the conservative equations as both arrive at the same final statement of the entropy flux in the normal direction at the boundary. However, we find that the nonlinear energy analysis reveals additional information to specifically see how each equation variable $(h,u,v)^{T}$ contributes to the boundary statement rather than the single expression that results form the entropy contraction. In essence, the entropy contraction “hides” an intermediate form of the boundary information that is revealed through the nonlinear energy analysis. This intermediate form is particularly useful to examine different open boundary configurations.

3 Semi-boundedness and energy estimates for initial boundary value problems

Before proceeding with the nonlinear energy analysis we review a few important concepts relevant for energy boundedness and well-posedness.

Consider the generic initial boundary value problem (IBVP)

[TABLE]

In (3.1), $\mathbf{q}$ is the solution, $\mathsf{D}$ is a (possibly nonlinear) differential operator, $\mathsf{L}$ is the boundary operator, $\Omega$ is the domain, $\partial\Omega$ its boundary and $f$ is initial data.

For two vector functions $\mathbf{u}$ and $\mathbf{v}$ defined on $\Omega$, we introduce the inner product and corresponding norm

[TABLE]

where $\mathcal{P}$ is a positive definite symmetric matrix. Further, we require the concept of semi-boundedness for linear problems [5, 48, 7] and its extended definition to include nonlinear problems [49]. For linear problems we consider a space of functions and in the nonlinear case a set of functions.

Definition 3.1**.**

Let $\mathbb{V}$ be the space (or set) of all differentiable functions satisfying the boundary conditions $\mathsf{L}\mathbf{u}=0$. The differential operator $\mathsf{D}$ in (3.1) is semi-bounded if for all $\mathbf{u}\in\mathbb{V}$

[TABLE]

holds, where $\lambda$ is a constant independent of $\mathbf{u}$.

Definition 3.2**.**

The differential operator $\mathsf{D}$ in (3.1) is maximally semi-bounded if it is semi-bounded in the function space (or function set) $\mathbb{V}$ but not in any space (set) with fewer boundary conditions.

If $\mathsf{D}$ is maximally semi-bounded, we can estimate the solution to (3.1) since

[TABLE]

Remark 3.3*.*

The relation (3.2) for non-positive constants $\lambda$, yields a weak estimate. Homogeneous boundary conditions which lead to $\lambda>0$ provide an estimate also for non-zero boundary data [20, 49].

Next, we relate the concepts of semi-boundedness presented above to well-posedness which we roughly define as follows. The IBVP (3.1) is well posed if i) the solution depends continuously on the given initial data (i.e. a bound like in (3.3) exist), ii) the solution is unique and iii) the solution exists.

For linear IBVPs with smooth initial data and a maximally semi-bounded operator $\mathsf{D}$, well-posedness follows directly. Uniqueness follows by considering the difference between two solutions with the same data. The bound on the original problem bounds the difference and show that it is zero. Existence follows since no over-specification of boundary conditions is done (see [50] for more details on how to determine the minimal number). Hence maximal semi-boundedness is a sufficient condition for well-posedness of linear IBVPs.

For nonlinear IBVPs we consider the analogous definition. We seek to determine the minimal number of boundary conditions such that the differential operator $\mathsf{D}$ is semi-bounded and lead to an energy estimate. As in the linear case, we then say that the operator is maximally semi-bounded. This is largely motivated by the linear theory, since the well-posedness theory for nonlinear hyperbolic systems is incomplete. Maximal semi-boundedness in the nonlinear case lead to an energy estimate, and hence it is a necessary condition for well-posedness. However, since neither uniqueness nor existence follows, it is not a sufficient condition.

4 Energy and entropy analysis

Instead of using the energy (or entropy) conservation law (2.5) directly, we will apply the classical energy method to the non-conservative system (2.2), which first require a suitable symmetrization matrix [19, 8]. Guided by the linear analysis in [8], we select the symmetrizer

[TABLE]

The matrix $\mathcal{S}$ simultaneously symmetrizes the flux matrices $\mathcal{A}$ and $\mathcal{B}$ as

[TABLE]

Further, the symmetrization matrix directly relates the solution energy to the specific total energy (2.3) as

[TABLE]

where $\mathcal{P}=\mathcal{S}^{2}$.

Now we are equipped to apply the classical energy method [48, 7, 4, 20, 21]. We pre-multiply (2.2) by $\mathbf{q}^{T}\mathcal{P}$ to obtain

[TABLE]

From the skew-symmetry of the Coriolis matrix we immediately see that $\mathbf{q}^{T}\mathcal{P}\mathcal{C}\mathbf{q}=0$. The flux matrices are now symmetrized and take the form

[TABLE]

We seek to rewrite (4.2) with complete derivatives, and use the relations

[TABLE]

The expression (4.2) becomes

[TABLE]

4.1 Nonlinear energy analysis

To begin, we compute the derivatives of the matrices in (4.4) to be

[TABLE]

which gives, noting the continuity equation (the first equation) from (2.2),

[TABLE]

It is now straightforward to compute

[TABLE]

Therefore, (4.4) becomes

[TABLE]

We continue by introducing the norm $\|\mathbf{q}\|^{2}_{\mathcal{P}}=\int_{\Omega}\mathbf{q}^{T}\mathcal{P}\mathbf{q}\,\mathrm{d}x\mathrm{d}y$, integrating (4.5) over $\Omega\subset\mathbb{R}^{2}$, and apply Gauss’ theorem to find

[TABLE]

where $u_{n}=un_{x}+vn_{y}$ is the normal velocity to the boundary $\partial\Omega$ and

[TABLE]

To simplify the upcoming analysis, we rotate the solution vector $\mathbf{q}$ as well as the matrix $\widetilde{\mathcal{A}}$ into the normal and tangential velocity directions, using

[TABLE]

such that

[TABLE]

and

[TABLE]

The energy evolution equation (4.6) can now be rewritten as

[TABLE]

Next, we examine the effect of the additional scalar term in (4.7) introduced by the nonlinearity of the problem.

Lemma 4.1**.**

There are three unique ways to rewrite the scalar term in (4.7) as a quadratic form

[TABLE]

where

[TABLE]

Proof.

See A. ∎

The results of Lemma 4.1 allows the reformulation of (4.7) to become

[TABLE]

Remark 4.2*.*

If we use the linear analysis and neglect the matrices $\mathcal{N}_{i}$, the contribution to the energy will be different, which leads to other boundary procedures [8, 44] (that will be discussed in detail in Section 4.3).

To develop an energy estimate we must bound the surface integral term in (4.10). To find such a bound we rotate the boundary matrix $\widehat{\mathcal{A}}-\mathcal{N}_{i}$ into a diagonal form [19, 20, 50, 21]. For each available boundary matrix, we introduce a solution dependent rotation matrix and corresponding vector of linearly independent rotated variables, i.e.,

[TABLE]

Each rotation matrix $\mathcal{R}_{i}$, $i=1,2,3$ transform the boundary matrix into a diagonal form

[TABLE]

such that the boundary integral term in (4.10) becomes

[TABLE]

We collect the rotations of the form (4.11) for the three boundary matrices:

$\widehat{\mathcal{A}}-\mathcal{N}_{1}$

The terms in $\mathcal{R}_{1}$ are $\alpha_{1}=1$, $\beta_{1}=-1/u_{n}$ and $\gamma_{1}=\delta_{1}=1/u_{n}$. The rotated variables are $\mathbf{w}_{1}=\left[\phi,u_{n}^{2}+\phi,u_{n}u_{s}\right]^{T}$ with the diagonal matrix

[TABLE]

$\widehat{\mathcal{A}}-\mathcal{N}_{2}$

For $\mathcal{R}_{2}$, the values $\alpha_{2}=-2$ and $\beta_{2}=\gamma_{2}=\delta_{2}=1/u_{n}$ in (4.11) give the rotated variables $\mathbf{w}_{2}=\left[-\phi/2,u_{n}^{2}+\phi/2,u_{n}u_{s}\right]^{T}$ and diagonal matrix

[TABLE]

$\widehat{\mathcal{A}}-\mathcal{N}_{3}$

The third matrix $\mathcal{R}_{3}$ takes the values

[TABLE]

which yields the rotated variables $\mathbf{w}_{3}=\left[\phi,u_{n}^{2},u_{n}u_{s}\right]^{T}$ and diagonal matrix

[TABLE]

Remark 4.3*.*

The contracted, rotated forms $\mathbf{w}_{i}^{T}\mathcal{D}_{i}\mathbf{w}_{i}$, $i=1,2,3$ are equivalent. However, only the boundary matrix associated with $\mathcal{N}_{1}$ can be rotated such that one can, by inspection, determine the number of positive and negative diagonal values. The reason is that the funtional dependence of the diagonal values in (4.12) can be factored out, clearly displaying the signs of the diagonal values. In contrast, the entries in (4.13) and (4.14) are complicated functions of the solution variables $\phi$ and $u_{n}$, that cannot be factored out. Therefore, determining the number of positive and negative entries in either $\mathcal{D}_{2}$ or $\mathcal{D}_{3}$ requires careful examination of the size and sign of the normal velocity in comparison to the geopotential. In (4.13) and (4.14) one must consider the diagonal values and the rotated variables at the same time. One cannot, as in (4.12), first determine the sign of the diagonal values and next combine the variables to create boundary conditions.

For the upcoming analysis of stable boundary conditions in Section 5 we select the first rotation to have (4.12), and (4.10) becomes

[TABLE]

4.2 Nonlinear entropy analysis

It is interesting to examine how the statement (4.5) compares to the energy and entropy conservation law (2.5). Due to the construction of the symmetrization matrix (4.1), we have the time evolution of the specific total energy

[TABLE]

We look at the energy flux contribution in the $x$-direction to find

[TABLE]

Similarly, in the $y$-direction we have

[TABLE]

So, we find that (4.5) is equivalent to the entropy conservation law (2.5). Integrating (2.5) over a general domain and applying Gauss’ theorem [51] yield

[TABLE]

with the outward pointing normal vector $\vec{n}=(n_{x},n_{y})^{T}$. The energy conservation law (4.16) is the same as (4.15) albeit derived from a different perspective. We expand the boundary term to highlight the similarity to (4.15)

[TABLE]

Remark 4.4*.*

For non-smooth solutions, the specific total energy conservation law (4.16) becomes the energy (or entropy) inequality

[TABLE]

which states that the specific total energy (or entropy) dissipates in the presence of discontinuities [38, 36]. It is the typical entropy condition for conservation laws [52, 53].

Remark 4.5*.*

The analysis above shows that the nonlinear energy and entropy analysis are similar, but that the energy analysis provides more details. In particular, it provides the matrix $\mathcal{P}$ forming the norm, the symmetrized matrices $\mathcal{P}\mathcal{A}$ and $\mathcal{P}\mathcal{B}$ in (4.3) (that depend on space and time but have the same form as in the linear analysis [8]) and an explicit matrix-vector formulation of the boundary term. The energy analysis also shows that the difference between the linear and nonlinear analysis depends only on $\mathcal{N}_{i}$ in (4.8). The entropy analysis on the other hand is more direct and easily applied. Having established the similarities and differences between the energy and entropy analysis, we in the following refer to both when we state energy (or entropy).

4.3 Linear energy analysis

The standard linear energy analysis of the SWE is contained within equation (4.4). This is because the nonlinearity of the problem produces the last term in (4.4) which originates from the space- and time-dependent symmetrizer $\mathcal{S}$ in (4.1). If one linearizes the equations (2.2) around a constant background state, the coefficient matrices $\mathcal{A}$ and $\mathcal{B}$ as well as the symmetrizer $\mathcal{S}$ become constants with the same functional form as in the nonlinear case and the final term in (4.4) vanishes.

Ignoring the nonlinear term $\mathcal{N}_{i}$ and performing such a linearization changes the form of the boundary contribution from (4.10) to

[TABLE]

where $\overline{\phi}$ and $\bar{u}_{n}$ are constant mean states for the geopotential and normal velocity. We use a rotation like (4.11) to bring the linearized matrix term in (4.17) into a diagonal form

[TABLE]

This also gives the linearized rotated variables $\overline{\mathbf{w}}=\overline{\mathcal{R}}^{-1}\widetilde{\mathbf{q}}=[\phi\,,\,\bar{u}_{n}u_{n}+\phi\,,\,\bar{u}_{n}u_{s}]^{T}$. Now the boundary term can be compactly written as

[TABLE]

Next, we manipulate the first term of the boundary integral in (4.18) to get

[TABLE]

Now the sign of the diagonal matrix shifts with the sign of $\bar{u}_{n}^{2}-\overline{\phi}$. When $\bar{u}_{n}^{2}<\overline{\phi}$ the flow is subcritical whereas when $\bar{u}_{n}^{2}>\overline{\phi}$, we have supercritical flow [44]. The matrix in (4.19) is scaled by the value $\overline{\phi}/\bar{u}_{n}$ which further influences the overall sign of each term. So we see that the signs of the matrix in (4.19) progresses through four states depending on the direction and magnitude of the normal velocity [54]:

Supercritical inflow where $\bar{u}_{n}<0$ yields three negative values.

• 2.

Subcritical inflow where $\bar{u}_{n}<0$ yields two negative values.

• 3.

Supercritical outflow where $\bar{u}_{n}>0$ yields one negative value.

• 4.

Subcritical outflow where $\bar{u}_{n}>0$ yields zero negative values.

The impact such sign changes have on the number of boundary conditions is discussed in the next Section. Linear boundary treatments for the SWE were originally developed by Oliger and Sundström [8] and expanded upon elsewhere, e.g., by Ghader and Nordström [42] and Huang and Temam [54].

Remark 4.6*.*

Using lifting approaches, e.g. [19], it is straightforward to apply the results from the forthcoming analysis and implement weak boundary conditions for numerical calculations [18, 19, 20, 21, 22, 23, 24, 25].

5 Energy and entropy stable boundary conditions

The sign of the normal velocity determines whether there are inflow or outflow conditions at the domain boundary. That is, $u_{n}<0$ corresponds to inflow conditions, $u_{n}>0$ to outflow conditions and $u_{n}=0$ to glancing conditions (whose prototypical physical example is a rigid boundary). Energy stable open boundary conditions of Dirichlet type for a general domain must be chosen such that the reformulated boundary term in (4.12) and (4.15)

[TABLE]

is made positive semi-definite. That is, a minimum number of boundary conditions are determined so that

[TABLE]

holds, which means that the operator $\mathcal{D}$ is maximally semi-bounded.

Above we know that the geopotential $\phi$ is strictly positive by the physical requirements, and we have supressed the subscript $1$. In (5.1), $\mathbf{w}^{+}$ contains the outgoing boundary information and $\mathbf{w}^{-}$ the incoming information. The diagonal blocks $\mathcal{D}^{+}$ and $\mathcal{D}^{-}$ contain the positive and negative diagonal entries, respectively. Note that the signs of the entries in the matrix $\mathcal{D}$ (4.12) depend entirely on the sign of the normal velocity. This clear separation of the positive and negative contributions to the energy at the boundary provides a suitable setting (not present in (4.13) or (4.14)) for discussing boundary conditions for the nonlinear problem. It will also serve to illustrate the difference between a linear and nonlinear analysis.

Our goal is to guarantee energy stability at open boundaries with a minimal number of boundary conditions. As a starting point, we follow the linear analysis path and assume that the number of boundary conditions required is equal to the size of $\mathcal{D}^{-}$ [19, 50]. However, in the course of the forthcoming analysis we will demonstrate that this assumption can be misleading. That is, it may be the case that more (or fewer) boundary conditions are required to guarantee that the inequality (5.2) holds. The subtlety of the nonlinear analysis is that the entries of the matrix $\mathcal{D}$ in general depend upon the solution quantities (for instance the normal velocity), whereas such a dependency does not exist in the linear case [19].

Remark 5.1*.*

As discussed in Section 3, an energy estimate can always be obtained by specifying too many boundary conditions [19, 50]. However, in general that means that existence cannot be obtained for linear problems. The precise existence requirement is not known for nonlinear problems, but it is reasonable to assume (unless proven otherwise) that similar requirements as for linear problems hold.

A general form of the boundary conditions is [20]

[TABLE]

where $\mathsf{R}$ is a coefficient matrix with the number of rows equal to the minimal number of required boundary conditions. The vector $\mathbf{g}$ contains known external data at the boundary. Essentially, the boundary condition (5.3) represents incoming boundary information as a linear combination of the outgoing information and external data. With the rewritten stability condition (5.1), its corresponding inequality (5.2) and the general boundary condition (5.3) we have

Theorem 5.2**.**

The inhomogeneous, general boundary conditions (5.3) for the nonlinear problem are energy stable provided

[TABLE]

and there exists a positive semi-definite matrix $\Gamma$ such that

[TABLE]

For homogeneous boundary conditions, (5.4) can be relaxed to

[TABLE]

and (5.5) becomes redundant.

Proof.

See [20].∎

Next, we will discuss the different types of boundaries that can exist.

5.1 Inflow boundaries

First, we consider inflow at open boundaries, i.e., $u_{n}<0$, $|u_{n}|>\delta$ where $\delta$ is a nonzero, possibly small, value. From the form of the boundary contribution (5.1) we require two boundary conditions since

[TABLE]

Further, from (5.3) we have for the inflow boundary case that

[TABLE]

From the first statement of Theorem 5.2 we find that

[TABLE]

becomes a scalar. Thus, the coefficients must satisfy the constraint $R_{1}^{2}+R_{2}^{2}<1$ to guarantee positivity. Additionally, we compute the matrix $\Gamma$ where

[TABLE]

where $C=R_{1}^{2}+R_{2}^{2}-1$. To determine if $\Gamma$ is positive semi-definite we compute its eigenvalues to find

[TABLE]

which are both positive and bounded because the coefficients for $\mathsf{R}$ are constrained to lie within the unit disc from (5.8) and $|u_{n}|>\delta$.

Remark 5.3*.*

The analysis above shows that two boundary conditions are sufficient to guarantee a bound on the energy for an inflow boundary. However, it is conceivable that uniqueness might require three. As mentioned above, these issues are presently not clear for general nonlinear hyperbolic problems.

Remark 5.4*.*

From the linear setting (4.19), three boundary conditions are required to bound supercritical flow, i.e. $|\bar{u}_{n}|>\sqrt{\overline{\phi}}$ where $\sqrt{\overline{\phi}}=\sqrt{g\bar{h}}$ is the gravity wave speed, as

[TABLE]

For subcritical flow, i.e. $|\bar{u}_{n}|<\sqrt{\overline{\phi}}$, only two are required [54] and the diagonal matrices are

[TABLE]

In the nonlinear case, the bound is not dependent on the magnitude of the flow and always requires two boundary conditions. From the details given above for the linear case we see that the minimum number of boundary conditions needed for a bound is at least two but possibly three depending on the flow regime. Consequently, the linear and nonlinear analysis are not consistent.

5.2 Outflow boundaries

Next, we examine outflow conditions at open boundaries, i.e., $u_{n}>\delta>0$. From the statement (5.1) we require one boundary condition since

[TABLE]

where all entries of the matrix $\mathcal{D}$ are multiplied by the positive quantity $\phi/u_{n}$. Thus, the general boundary condition (5.3) takes the form

[TABLE]

To apply the result of Theorem 5.2 we first compute the matrix

[TABLE]

and calculate its eigenvalues to find

[TABLE]

Here, the constraint $R_{1}^{2}+R_{2}^{2}<1$ is, again, needed to guarantee positivity. Further, we compute the scalar $\Gamma$ to be

[TABLE]

which is also positive and bounded whenscaled by $\phi/u_{n}$ since the coefficients given in $\mathsf{R}$ are contained within the unit disc and $u_{n}>\delta$.

Alternatively, if we do not apply any boundary condition then the boundary term becomes

[TABLE]

which still satisfies the inequality (5.2) because of the positive normal velocity. That is, one could view an outflow boundary as a so-called “free” boundary [41] because no boundary conditions are required.

Remark 5.5*.*

Hence, to obtain an energy bound at outflow boundaries no boundary condition is required. However, uniqueness might still require the imposition of one condition (as was indicated by the eigenvalue analysis above).

Remark 5.6*.*

In the linear case (4.19) no boundary condition is needed for supercritical outflow because the diagonal matrices at the boundary are

[TABLE]

One boundary condition is needed for subcritical outflow in the linear case [42] with diagonal matrices

[TABLE]

For outflow boundaries, we find, as for inflow boundaries, that the linear and nonlinear analysis are inconsistent. In the nonlinear case, the magnitude of the normal velocity is irrelevant and a bound is obtained with no boundary condition.

5.3 Glancing boundaries

For glancing flow, i.e. $u_{n}=0$, there is a singularity in $\mathcal{D}$, however this causes no problem. First, if the flow transitions from outflow to glancing flow, i.e. $u_{n}\to 0^{+}$, then the formulation (5.9) shows that a (neutral) energy bound exists for $u_{n}=0$. Secondly, if the flow transfers from inflow into a glancing flow, i.e. $u_{n}\to 0^{-}$, then the formulation (5.7) implies that $R_{1}=1$ and $R_{2}=0$ and the statement (5.8) becomes

[TABLE]

which is consistent with (5.6). Hence, Theorem 5.2 implies that glancing boundaries are (neutrally) stable, and no singularities exist.

Remark 5.7*.*

The nonlinear and linear analysis are consistent for glancing boundaries, i.e. no boundary condition is required provided we know that $u_{n}=0$. The situation is similar if we instead specify that $u_{n}=0$ as a boundary condition, the consistency remains. One can intuitively understand this by noting that the additional scalar term in (4.8) vanishes for $u_{n}=0$.

5.4 Discussion

We demonstrated that the linear analysis will, in general, yield a different number of boundary conditions compared to the nonlinear problem. In the linear analysis, the number of boundary conditions directly corresponds to the number of negative eigenvalues [41, 42, 50, 44]. For the linear case this leads to maximal semi-boundedness and well-posedness, i.e. a bound on the solution as well as uniqueness and existence. However, in the nonlinear analysis, the minimal number of boundary conditions required for an energy estimate could be different than that determined by the linear case. In fact, it was possible to obtain an energy estimate by prescribing only two boundary conditions for the inflow case and zero for the outflow case.

We also found that the magnitude of the normal velocity played no role in determining the number of boundary conditions in the nonlinear case. This is quite surprising, again not consistent and quite different from the results of the linear analysis.

Furthermore, we found that the subtleties of the nonlinear energy analysis are not obvious in the traditional entropy analysis. This is because the contraction of the hyperbolic system of conservation laws into entropy space occurs “too quickly” and hides information necessary in the boundary condition analysis. We see this directly in the boundary statement from the entropy analysis where the sign of the boundary contribution is dictated entirely by the normal velocity. This provides no clue as to how many conditions one should apply at inflow, in contrast to the number of boundary conditions indicated by the nonlinear energy analysis.

From the linear and nonlinear energy analysis in this work we demonstrated that the minimal number of boundary conditions required to guarantee a semi-bounded operator differ. In the linear analysis, maximal semi-boundedness together with certain smoothness assumptions on the initial data is enough to imply well-posedness. For the nonlinear problem maximal semi-boundedness leads to an energy estimate but, unfortunately, it is not known whether it leads to well-posedness since uniqueness and/or existence might be lacking. This is because the underlying well-posedness theory for hyperbolic systems of conservation laws, like the SWEs, is incomplete. In the linear case, the minimal number of boundary conditions together with an energy estimate is sufficient to guarantee well-posedness [55]. For the nonlinear system, maximal semi-boundedness is a necessary requirement for well-posedness but it may not be sufficient.

Specifically regarding the SWEs, the analysis in this work considered the model without bottom topography terms. For practical simulations in ocean modeling a nonzero bottom topography must be properly discretized and accounted for to develop an energy estimate, see, e.g., [24]. This is because the bottom topography contributes to the potential energy in the specific total energy (2.3). The analysis in this work is amenable to the case of nonzero bottom topography for the SWEs; however, it requires one to introduce the total geopotential $\Phi=g(h+b)$ as a primitive variable in (2.2) and invoke the assumption that the bottom topography is steady, i.e., $b_{t}=0$. We plan to extend the nonlinear energy analysis for the SWEs with nonzero bottom topography in future work.

6 Concluding remarks

We derived energy stable boundary conditions for the nonlinear shallow water equations. The classical energy method was applied to determine the stability conditions. We demonstrated that there are several possible ways to rewrite terms from the nonlinear analysis into a matrix form. Further, we demonstrated that one of these matrix forms was preferable for the boundary condition analysis because it could be rotated into a particularly simple diagonal form. It was shown that the traditional energy analysis is consistent with the entropy stability analysis where the specific total energy plays the role of a generalized entropy function. We also found that the energy analysis provides additional details but is more cumbersome, than the entropy analysis. The details brought forward by the nonlinear energy analysis also made it possible to pinpoint the difference with the linear analysis.

Specifically, this work investigated and showed:

The differences and similarities between the nonlinear energy and entropy methods. 2. 2.

Why the linear and nonlinear energy analysis differ and how this affects the boundary conditions. 3. 3.

That a linear analysis cannot be trusted to determine the boundary conditions for a nonlinear problem. 4. 4.

That the flow magnitude does not influence the number of boundary conditions in the nonlinear case. 5. 5.

That the eigenvalues and their corresponding signs are not the end all, be all to indicate the number of boundary conditions necessary for semi-boundedness in the nonlinear context.

The use of a nonlinear energy analysis offers two interesting future avenues of study. Firstly, for the nonlinear shallow water equations there are interesting theoretical questions concerning open boundary conditions for flow over non-constant bottom topography terms as well as possible friction on said bathymetries. This is because the bottom topography contributes to the potential energy of the shallow water equations and the energy analysis requires further investigation. Also, the addition of friction introduces the gradient of the geopotential (or water height) into the model. Therefore, this may expand the possible types of boundary conditions beyond the Dirichlet-type considered herein, to Neumann and Robin type.

Secondly, we want to expand this investigation to other nonlinear hyperbolic systems, like the compressible Euler equations, where the mathematical entropy function is not equal to the specific total energy. The marriage of traditional energy analysis and nonlinear entropy analysis requires careful attention, particularly when the mathematical entropy function is not the total energy of the system. However, the close relationship between the linear/nonlinear energy analysis and the entropy analysis could possibly pave the way for truly energy and/or entropy stable approximations.

Acknowledgments

Jan Nordström was supported by Vetenskapsrådet, Sweden [award no. 2018-05084 VR] and Andrew Winters by Vetenskapsrådet, Sweden [award no. 2020-03642 VR].

Appendix A Proof of Lemma 4.1

Proof.

To rewrite the scalar term as a quadratic form we assume the existence of a symmetric matrix of the form

[TABLE]

The tangential component of the velocity $u_{s}$ has no role in the remaining scalar term (and play no part in the boundary contribution to the energy [19, 20]). Thus, $n_{3}=n_{5}=n_{6}=0$. Therefore, the matrix ansatz for $\mathcal{N}$ becomes

[TABLE]

This gives three choices on which constant to solve for in the above relationship. First, we solve for $n_{1}$ to find

[TABLE]

This solution yields the matrix

[TABLE]

where the values of $n_{2}$ and $n_{4}$ are arbitrary. However, when multiplying from the left and right with $\widetilde{\mathbf{q}}$ this remaining matrix term vanishes, i.e.,

[TABLE]

Therefore, the values of $n_{2}$ and $n_{4}$ are irrelevant to the boundary contribution and we arrive at the first matrix in (4.9).

Alternatively, we can elect to solve for $n_{2}$ in (A.1) to obtain

[TABLE]

This produces the matrix

[TABLE]

where, now, the values of $n_{1}$ and $n_{4}$ remain arbitrary. But when substituting the matrix $\mathcal{N}_{2}$ into the quadratic form ansatz the remaining matrix vanishes

[TABLE]

and the values of $n_{1}$ and $n_{4}$ are irrelevant, yielding the second matrix in (4.9).

Lastly, we can solve for the value

[TABLE]

to find the third matrix form

[TABLE]

Substituting the matrix $\mathcal{N}_{3}$ into the quadratic form, the second matrix term vanishes, i.e.,

[TABLE]

and we find that the arbitrary values of $n_{1}$ and $n_{2}$ are irrelevant. ∎