# Energy versus entropy estimates for nonlinear hyperbolic systems of equations

We compare and contrast information provided by the energy analysis of Kreiss and the entropy theory of Tadmor for systems of nonlinear hyperbolic conservation laws. The two-dimensional nonlinear shallow water equations are used to highlight the similarities and differences since the total energy of the system is a mathematical entropy function. We demonstrate that the classical energy method is consistent with the entropy analysis, but significantly more fundamental as it guides proper boundary treatments. In particular, the energy analysis provides information on what type of and how many boundary conditions are required, which is lacking in the entropy analysis. For the shallow water system we determine the number and the type of boundary conditions needed for subcritical and supercritical flows on a general domain. As eigenvalues are augmented in the nonlinear analysis, we find that a flow may be classified as subcritical, but the treatment of the boundary resembles that of a supercritical flow. Because of this, we show that the nonlinear energy analysis leads to a different number of boundary conditions compared with the linear energy analysis. We also demonstrate that the entropy estimate leads to erroneous boundary treatments by over specifying and/or under specifying boundary data causing the loss of existence and/or energy bound, respectively. Our analysis reveals that the nonlinear energy analysis is the only one that provides an estimate for open boundaries. Both the entropy and linear energy analysis fail.

## Authors

• 9 publications
• 10 publications
09/21/2021

### Nonlinear and Linearised Primal and Dual Initial Boundary Value Problems: When are they Bounded? How are they Connected?

Linearisation is often used as a first step in the analysis of nonlinear...
04/16/2020

### A Crank-Nicolson type minimization scheme for a hyperbolic free boundary problem

We consider a hyperbolic free boundary problem by means of minimizing ti...
06/21/2021

### The implementation of a broad class of boundary conditions for non-linear hyperbolic systems

We propose methods that augment existing numerical schemes for the simul...
05/10/2021

### On the Theoretical Foundation of Overset Grid Methods for Hyperbolic Problems: Well-Posedness and Conservation

We use the energy method to study the well-posedness of initial-boundary...
04/26/2020

### Singular Diffusion with Neumann boundary conditions

In this paper we develop an existence theory for the nonlinear initial-b...
10/09/2020

### Deep Autoencoder based Energy Method for the Bending, Vibration, and Buckling Analysis of Kirchhoff Plates

In this paper, we present a deep autoencoder based energy method (DAEM) ...
04/06/2020

### Stable Boundary Conditions and Discretization for PN Equations

A solution to the linear Boltzmann equation satisfies an energy bound, w...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

Two general avenues are available in order to obtain estimates for hyperbolic systems of conservation laws. They are the energy method of Kreiss [16, 17] and the entropy stability theory of Tadmor [32, 33]. Traditionally, the method of Kreiss has been applied to linearized versions of systems of hyperbolic equations in order to develop boundary treatments that lead to an energy estimate. In practice, these boundary conditions are needed to develop energy stable numerical approximations that weakly impose boundary information, e.g., through simultaneous approximation terms [4, 23, 24, 25] or numerical interface flux functions [14, 19, 40, 41]. In contrast, the method of Tadmor has been applied to nonlinear hyperbolic systems on domains with periodic boundary conditions (or infinite domains) in order to obtain entropy conservation. This makes the investigation of entropy conservation similar to the classical von Neumann stability analysis in the sense that boundaries are ignored [12, 29]. By adding dissipation, entropy stability is obtained for periodic or infinite domains [33]. Hence, the use of entropy stability theory to develop provably stable boundary conditions has been limited [14, 27, 31]. The main strengths and weaknesses of these two approaches are: The energy analysis of Kreiss provides boundary conditions, but it is difficult to apply in the nonlinear case. It is straightforward to apply the entropy analysis of Tadmor for general systems of nonlinear equations, but it does not provide boundary condition information.

The focus of this work is to examine the connection between the classical energy method of Kreiss and Tadmor’s entropy analysis. Particular focus is given to comparing and contrasting these two strategies for deriving stable boundary treatments. In doing so, we demonstrate that the nonlinear energy method of Kreiss provides fundamental information about the hyperbolic system, that aid in choosing a minimal number of suitable boundary conditions, which is required for existence [11, 12, 25], and energy stability. The entropy stability theory of Tadmor is often portrayed as a nonlinear generalization of energy stability analysis [21]. However, we show that it gives no information about the characteristics of the hyperbolic problem and offers little or even erroneous information as to what type of and how many boundary conditions are needed on a general bounded domain [30, 31]. We also include the linear energy analysis in our comparison and discuss its weaknesses.

To perform this comparison we consider the two-dimensional nonlinear shallow water equations (SWEs). Shallow water models are of particular interest for flow configurations where the vertical scales of motion are much smaller than the horizontal scales, such as in rivers or lakes [37, 38]

. The SWEs are a system of nonlinear hyperbolic partial differential equations that represent the conservation (or balance) of mass and momentum, depending on the forces, e.g. bottom friction,

[6, 19, 37, 41]. An auxiliary conserved quantity, not explicitly built into the SWEs, is the total energy of the system. This additional conservation law can be used to create a stability estimate for the total energy [9, 34, 40] or build numerical approximations that respect the evolution of the total energy [8, 9, 22]. For the shallow water system the total energy also acts as a mathematical entropy function and fits into the entropy analysis framework of Tadmor [33].

Thus, the total energy and analysis of it for the SWEs act as a bridge between the classical energy method of Kreiss and Tadmor’s entropy stability theory. We apply the energy method [12, 16, 17] and derive a bound of the total energy, which for the SWEs is a particular scaled version of the norm of the solution. In the following, we will demonstrate that the energy method is consistent with the mathematical entropy analysis, but also that it provides additional information and guidance with respect to boundary treatments for the nonlinear problem. Investigations into energy consistent boundary conditions for the linearized SWEs are many [2, 3, 10, 20, 26, 37]. The linear analysis leads to a well-posed linear initial boundary value problem [10, 26, 37]. These linear boundary conditions can then be applied to the nonlinear case [10, 36]. However, as we will show, there are situations where the linear analysis cannot be applied to the nonlinear case. Similar to the entropy analysis, linear boundary treatments do not necessarily provide an energy estimate for the nonlinear problem.

The paper is organized as follows: The SWEs are given in Sect. 2. An estimate of the total energy for the shallow water system is provided in Sect. 3 using a minimal number of boundary conditions. In Sect. 4, we provide details, analysis, and discussion of the general open boundary conditions for the two-dimensional nonlinear SWEs in subcritical and supercritical flow regimes. In particular, we discuss the differences between the results from the linear energy analysis, the nonlinear energy analysis, and the entropy analysis. Concluding remarks are drawn in the final section.

## 2. Shallow water equations

We begin with the two-dimensional SWEs over a flat bottom topography written in conservative form [38]

 (2.1) ht+(hu)x+(hv)y =0, (hu)t+(hu2+g2h2)x+(huv)y−fhv =0, (hv)t+(huv)x+(hv2+g2h2)y+fhu =0,

which includes the continuity and momentum equations. Here is the water height, and are the fluid velocities in the - and -directions, and is the gravitational constant. The system of equations (2.1) are derived under the physical requirement that the water height [37, 38]. Additionally, we include the influence of Coriolis forces with the parameter which, for convenience, is assumed to be a constant. In practical applications is typically a function of latitude [37], which would not affect the subsequent energy analysis in this work.

In order to apply the (classical) energy method it is convenient to work with the equivalent non-conservative form of the governing equations [18, 23, 24, 26] which is

 (2.2) ϕt+ϕxu+ϕyv+ϕux+ϕvy =0, ut+uux+vuy+ϕx+fv =0, vt+uvx+vvy+ϕy−fv =0.

In (2.2), we formulate the non-conservative equations in terms of the geopotential to simplify the analysis [26]. Note, that 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

and

 qt+Aqx+Bqy+Cq=0,

where

 A=⎡⎢⎣uϕ01u000u⎤⎥⎦,B=⎡⎢⎣v0ϕ0v010v⎤⎥⎦,C=⎡⎢⎣00000f0−f0⎤⎥⎦.

The total energy (or entropy) of the SWEs is the sum of the kinetic and potential energy [9] where

 (2.3) ϵ=ϕ2g(u2+v2)+ϕ22g,

is an auxiliary conserved quantity of (2.1) or (2.2). The total energy (2.3) has associated total energy (or entropy) fluxes [8]

 (2.4)

that yield the total energy (or entropy) conservation law

 (2.5) ϵt+fϵx+gϵy=0.
###### Remark 2.1.

The Coriolis force is not present in the total energy equation (2.5). This agrees with the underlying physics of the problem because the Coriolis terms does not perform work on the fluid [37] and, thus, do not appear.

###### Remark 2.2.

In the entropy analysis of Tadmor, the total energy acts as a mathematical entropy function for the SWEs [8]. As such, it is possible to define a new set of entropy variables where

 s1=∂ϵ∂h,s2=∂ϵ∂(hu),s3=∂ϵ∂(hv).

Multiplying the conservative form of the SWEs (2.1) from the left with yields the auxiliary conservation law of the entropy function (in this case the total energy) (2.5

). This contraction of the conservative form of the SWEs into entropy space involves the chain rule and relies on certain compatibility conditions between the conservative fluxes from (

2.1) and the entropy fluxes (2.4) [8, 33].

## 3. Energy stability analysis

Next, we will apply the energy method to the non-conservative system (2.2), which require a suitable symmetrization matrix [23, 26]. Following [26], we select

 S=κ⎡⎢ ⎢⎣1000√ϕ000√ϕ⎤⎥ ⎥⎦,

where is a constant independent of the solution . The matrix simultaneously symmetrizes the flux matrices and

 AS\coloneqqSAS−1=⎡⎢ ⎢⎣u√ϕ0√ϕu000u⎤⎥ ⎥⎦,BS\coloneqqSBS−1=⎡⎢ ⎢⎣v0√ϕ0v0√ϕ0v⎤⎥ ⎥⎦.

To determine the scaling constant we examine the solution energy that will arise in the later analysis:

 (Sq)TSq=qTPq=κ2(ϕu2+ϕv2+ϕ2),

where . We want the solution energy to match the total energy (2.3)

 κ2(ϕu2+ϕv2+ϕ2)!=ϕu2+ϕv2+ϕ22g.

Therefore, we take and the final symmetrization matrix reads

 (3.1) S=1√2g⎡⎢ ⎢⎣1000√ϕ000√ϕ⎤⎥ ⎥⎦.

Now we are equipped to apply the classical energy method [11, 12, 16, 24, 25]. We pre-multiply (2.2) by to obtain

 (3.2) qTPqt+qTPAqx+qTPBqy+qTPCq=0.

From the skew-symmetry of the Coriolis matrix we immediately see that

. The flux matrices are now symmetrized and take the form

 PA=12g⎡⎢⎣uϕ0ϕϕu000ϕu⎤⎥⎦,PB=12g⎡⎢⎣v0ϕ0ϕv0ϕ0ϕv⎤⎥⎦.

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

 qTPqt =12(qTPq)t−12qTPtq, qTPAqx =12(qTPAq)x−12qT(PA)xq, qTPBqy =12(qTPBq)y−12qT(PB)yq.

The expression (3.2) becomes

 (3.3)

We compute the derivatives of the matrices to be

 Pt=12g⎡⎢⎣0000ϕt000ϕt⎤⎥⎦,(PA)x =12g⎡⎢⎣uxϕx0ϕxϕxu+ϕux000ϕxu+ϕux⎤⎥⎦, (PB)y =12g⎡⎢⎣vy0ϕy0ϕyv+ϕvy0ϕy0ϕyv+ϕvy⎤⎥⎦,

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

 Pt+(PA)x+(PB)y=12g⎡⎢⎣ux+vyϕxϕyϕx00ϕy00⎤⎥⎦.

It is straightforward to compute

 qT(Pt+(PA)x+(PB)y)q =12g[ϕ2(ux+vy)+2ϕϕxu+2ϕϕyv] =12g[(ϕ2u)x+(ϕ2v)y].

Therefore, (3.3) becomes

 (3.4) (qTPq)t+(qTPAq)x+(qTPBq)y−12g[(ϕ2u)x+(ϕ2v)y]=0.
###### Remark 3.1 (Connection to mathematical entropy analysis).

It is interesting to examine how the statement (3.4) compares to the energy conservation law created with the analysis tools in [8]. Due to the construction of the symmetrization matrix (3.1), we have the time evolution of the total energy

 (qTPq)t=(ϕu2+ϕv2+ϕ22g)t=ϵt.

Next, we look at the energy flux contribution in the -direction

 (qTPAq)x−12g(ϕ2u)x =(ϕu32g+ϕuv22g+ϕ2ug+ϕ2u2g)x−12g(ϕ2u)x =(ϕu2g(u2+v2)+ϕ2ug)x =fϵx.

Similarly, in the -direction we find

 (qTPBq)y−12g(ϕ2v)y=(ϕv2g(u2+v2)+ϕ2vg)y=gϵy,

So, we find that (3.4) is equivalent to the entropy conservation law (2.5).

###### Remark 3.2.

The main difference between (3.4) and (2.5) is the multitude of details and structure in (3.4) and the lack of it in (2.5).

Next, we examine the additional terms in (3.4) and try to incorporate them into the energy rate by rewriting the scalar terms into matrix-vector forms. We denote the required additional matrices by and , respectively. There are many possible ways to do this. We choose to follow a strategy where we require:

1. and to be simultaneously diagonalizable for any normal vector .

2. The scalar terms from (3.4) to be written in the normal direction as

 qTSˆNSq=qTP(N1nx+N2ny)q=ϕ22g(unx+vny).

The first requirement ensures that the matrices and

have the same eigenvectors. Combined with the matrix structure of the second requirement, this simplifies the derivation of boundary conditions, as will be demonstrated in Sect.

4. To create the matrices and we need

###### Lemma 3.3.

Consider two real matrices . If either the matrix or has distinct eigenvalues and the matrices commute

 WV=VW,

then the matrices are simultaneously diagonalizable

 W=RΛWR−1,V=RΛVR−1.
###### Proof.

See [15]. ∎

In the current analysis we know that the eigenvalues of are all distinct. Therefore, due to Lemma 3.3, it is sufficient to guarantee simultaneous diagonalizability if we can determine a matrix that commutes with . This leads to the following result.

###### Theorem 3.4.

The matrix

 ˆN=12⎡⎢ ⎢ ⎢⎣0nx√ϕny√ϕnx√ϕ00ny√ϕ00⎤⎥ ⎥ ⎥⎦

commutes with and is simultaneously diagonalizable with the same right eigenvector matrix , i.e.,

 ˆA=RΛˆARTandˆN=RΛˆNRT.

 N1=12⎡⎢⎣0ϕ0100000⎤⎥⎦,N2=12⎡⎢⎣00ϕ000100⎤⎥⎦,

commute with and , respectively. Furthermore, we can reformulate the scalar terms in (3.4) to be

 −(ϕ2u2g)x=−(qTPN1q)xand−(ϕ2v2g)y=−(qTPN2q)y.
###### Proof.

See Appendix A. ∎

From the result in Theorem 3.4, the energy equation (3.4) becomes

 (3.5)

Integrating (3.5) over , gives

We compactly write the time dependent term by introducing the norm

 ∥q∥2P=∫ΩqTPqdxdy,

for the symmetric positive definite matrix [23]. By applying Gauss’ theorem, we find

 ddt∥q∥2P+∮∂ΩqTP((A−N1)nx+(B−N2)ny)qdS=0,

with the outward pointing normal vector on the surface . We rewrite the line integral contribution in the normal direction to be

 (3.6) ddt∥q∥2P+∮∂ΩqTS(ˆA−ˆN)SqdS=0.

with the matrix

 ˆA=⎡⎢ ⎢ ⎢⎣unnx√ϕny√ϕnx√ϕun0ny√ϕ0un⎤⎥ ⎥ ⎥⎦,un=unx+vny,

and given in Theorem 3.4.

The eigenvalues of the matrices and are

 (λˆA)0 =un,(λˆA)±=un±√ϕ\coloneqqun±c, (λˆN)0 =0,(λˆN)±=±√ϕ2\coloneqq±c2,

where we define the wave celerity . From the construction of the matrix in Theorem 3.4 we know it has the same eigenvectors as , which are

 R\coloneqq[r+|r0|r−]=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣1√201√2nx√2−ny−nx√2ny√2nx−ny√2⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦.

This gives the eigendecomposition

 ˆA−ˆN=R(ΛˆA−ΛˆN)RT,

with

 ΛˆA=diag(un+c,un,un−c),ΛˆN=diag(c2,0,−c2).

From this information the expression (3.6) becomes

 (3.7) ddt∥q∥2P+ ∮∂ΩqTS(ˆA−ˆN)SqdS =ddt∥q∥2P+∮∂ΩqTSR(ΛˆA−ΛˆN)RTSqdS =ddt∥q∥2P+∮∂ΩwT(ΛˆA−ΛˆN)wdS=0,

now written in terms of the scaled characteristic variables

 (3.8) w\coloneqqRTSq=12√g⎡⎢ ⎢ ⎢⎣ϕ+√ϕ(unx+vny)√2ϕ(vnx−uny)ϕ−√ϕ(unx+vny)⎤⎥ ⎥ ⎥⎦=12√g⎡⎢⎣ϕ+cunc√2usϕ−cun⎤⎥⎦,

where is the tangential velocity.

###### Remark 3.5 (Relation to linear analysis).

We can recover the characteristic variables from the work in [10] if we take

 wLinear=√2gcw,

and then perform a linearization of the solution around a constant mean state.

The relation (3.7) implies that the nonlinear two-dimensional SWEs will be energy stable provided that the surface integral is made positive with a minimal number of energy stable boundary conditions.

###### Remark 3.6 (Entropy flux at the boundary).

A straightforward computation yields an alternative form of the boundary integral to be

 ∮∂ΩqTS(ˆA−ˆN)SqdS=∮∂Ωfϵnx+gϵnydS.

This shows that the nonlinear energy analysis is consistent with a weak form of the entropy conservation law (2.5) [40]. This form has been used to create entropy stable boundary conditions at non-penetrating walls [14, 27, 31]. However, as will be shown in Sect. 4.3, it offers no guidance as to what type of or how many boundary conditions should be imposed to guarantee energy or entropy stability, on general types of boundaries.

###### Remark 3.7 (Inconsistency with linear analysis).

If we use the linear analysis and neglect the terms contained in the matrices and , the contribution to the energy will contain an additional term at the boundary. In Sect. 4.3, we clarify what will happen if the linear results are applied in the nonlinear case.

## 4. Energy stable boundary conditions

The sign of the normal velocity determines whether there are inflow or outflow conditions at the domain boundary. That is, corresponds to inflow conditions and corresponds to outflow conditions. Additionally, we have that the argument of the surface integral from (3.7) is

 (4.1) wT(ΛˆA−ΛˆN)w =w21(un+c)+w22un+w23(un−c)−cw212+cw232 =w21(un+c2)+w22un+w23(un−c2) \coloneqqw21λ1+w22λ2+w23λ3,

where we introduce to be the new augmented eigenvalues.

###### Remark 4.1 (Solid wall boundary).

At a solid wall boundary (formally an outflow boundary) the normal velocity is zero and the statement (4.1) becomes

 w21λ1+w22λ2+w23λ3=(w21−w23)c2=0,

given the structure of the characteristic variables (3.8). Therefore, the solid wall boundary condition, , is energy (and entropy) stable for the nonlinear SWEs as pointed out in [26].

To determine energy stable open boundary conditions for a general domain

 (4.2) wT(ΛˆA−ΛˆN)w≥0,

must hold. The boundary conditions must be of Dirichlet type. Neumann and Robin type boundary conditions are not admissible due to the lack of gradients. We rewrite condition (4.2) into the form

 (4.3) wT(ΛˆA−ΛˆN)w=[w+w−]T[Λ+00Λ−][w+w−]≥0,

where contains the outgoing boundary information, the incoming information. The diagonal blocks and contain the positive and negative eigenvalues, respectively. Thi<s clear separation of the positive and negative eigenvalue contributions as well as the characteristic variables into incoming and outgoing boundary information provides a suitable setting to discuss energy stable boundary conditions for the nonlinear problem.

To start, we note that in order to bound (4.3), the minimal number of boundary conditions required is equal to the size of [23, 24]. The most general form of the boundary conditions is written [23, 24]

 (4.4) w−=Rw++gext,

where is a coefficient matrix with the number of rows equal to the minimal number of required boundary conditions (equal to the size of ). The vector contains known external data from the boundary. Essentially, the boundary condition (4.4) represents the incoming information as a linear combination of the outgoing information and boundary data. With the rewritten stability condition (4.3) and the general boundary condition (4.4) we have

###### Theorem 4.2.

The general boundary conditions (4.4) for the nonlinear problem lead to energy stability and a bound provided

 Λ++RTΛ−R≥0.

See [24]. ∎

###### Remark 4.3.

The proof in [24] for inhomogeneous boundary conditions involves some technical boundedness issues that we avoid in the current discussion. This is because it is straightforward to homogenize the boundary conditions [12, 29] and prove energy stability for the system augmented with a forcing function. Therefore, in the forthcoming analysis we take .

To present and discuss energy stable boundary conditions we introduce the Froude number for the normal flow component

 (4.5) Fr=|un|c,

which is a translation of the Mach number into the shallow water flow context [37, 41]. It serves to classify flow regimes as supercritical (or torrential) when and subcritical (or fluvial) when . Most shallow water flows are subcritical where, typically, [37]. However, under special circumstances, like near a dam failure or over a non-constant bottom topography, the flow can become supercritical, e.g. [1, 5, 41].

The sign and magnitude of the normal flow velocity determine the number of boundary conditions that are needed. To construct energy stable boundary conditions we must satisfy (4.2) with a minimal number of boundary conditions, i.e. a minimal number of rows in [23, 24]. We collect the different cases:

Supercritical inflow where :

The eigenvalues are all negative. This corresponds to

 (4.6) w+=0,w−=w,Λ+=0,Λ−=diag(λ1,λ2,λ3),

and we need three boundary conditions.

Supercritical outflow where :

The three eigenvalues are positive such that

 (4.7) w+=w,w−=0,Λ+=diag(λ1,λ2,λ3),Λ−=0.

Therefore, zero boundary conditions are required.

Subcritical inflow where :

We have and . Hence,

 (4.8) w+=w1,w−=(w2,w3)T,Λ+=λ1,Λ−=diag(λ2,λ3),

and we need two boundary conditions.

Subcritical inflow where :

Here the eigenvalues are all negative, i.e., . The sign of changed due to the strength of the normal velocity component. We again have situation (4.6) and require three boundary conditions, as in the supercritical case.

Subcritical outflow where :

We have and . Hence,

 (4.9) w+=(w1,w2)T,w−=w3,Λ+=diag(λ1,λ2),Λ−=λ3,

and we need one boundary condition.

Subcritical outflow where :

All the eigenvalues are positive, i.e., . The sign of flipped due to the strength of the normal velocity component such that the characteristics match (4.7) and now zero boundary conditions are needed, as in the supercritical case.

Compared to the linear analysis, e.g. [3, 10, 37], we see that the relative magnitude of the normal velocity, , with regards to the wave celerity, , affect the number of boundary conditions differently. This is due to the augmentation of the eigenvalues in the stability condition (4.2) by the additional matrix terms given in Theorem 3.4. Interestingly, the number of boundary condition needed only changes in the subcritical flow regime. As an example, the boundary treatment of a subcritical flow match that of a supercritical flow for subcritical outflow where . As previously mentioned, most shallow water flows are subcritical, which highlights the importance of this new finding.

For comparison with the current nonlinear energy analysis, we return to the stability requirement provided by the standard entropy analysis tools briefly described in Remarks 2.2 and 3.6:

 (4.10) ∫ΩϵtdΩ+∮∂Ωfϵnx+gϵnydS=0.

Stability relies on the sign of the integrand . If we substitute the form of the entropy fluxes (2.4) in (4.10) we find

 (4.11) fϵnx+gϵny=(ϕ2g(u2+v2)+ϕ2g)un =12g(ϕu2+ϕv2+2ϕ2)un =12g⎡⎢⎣ϕuv⎤⎥⎦T⎡⎢⎣2un000un000un⎤⎥⎦⎡⎢⎣ϕuv⎤⎥⎦.

For non-penetrating walls, where , (4.10) is energy stable as discussed in Remark 4.1. This is formally a subcritical outflow boundary, so the number of boundary conditions, i.e. one, is correct. However, in the general case, the statement (4.11) gives, erroneously, that the number of boundary conditions is entirely defined by the sign of the normal velocity. It indicates that three boundary conditions are required for and zero boundary conditions for . We see that the wave celerity or the flow being classified as subcritical or supercritical do not even appear in the boundary treatment provided by the entropy analysis.

### 4.1. Supercritical inflow and outflow boundaries

The boundary treatment for supercritical flows for the nonlinear equations matches the linear analysis as well as the entropy analysis. For completeness, we restate the boundary conditions: At a supercritical inflow boundary we need three boundary conditions to satisfy (4.3). Because , the term from (4.4) with the coefficient matrix vanishes and the boundary condition has the form

 (4.12) w−=0,

or any such linear combination.

At supercritical outflow, and the conditions of Theorem 4.2 are automatically satisfied, i.e no boundary conditions are required. It is a so-called “free” boundary [3].

### 4.2. Subcritical inflow and outflow boundaries

First, we consider subcritical inflow where . The number of boundary conditions to specify depends on the magnitude of the Froude number.

If we are in the regime where , then we must prescribe two boundary conditions. Following the ansatz (4.4) the coefficient matrix is

 R=[γInθIn],

such that the two incoming characteristic components, , are written in terms of the outgoing information, ,

 (4.13) [w2w3]=[γInθIn]w1,

with unknown coefficients and .

For a subcritical inflow with the sufficient condition for nonlinear energy stability from Theorem 4.2 becomes

 (4.14) Λ++RTΛ−R=λ1+γ2Inλ2+θ2Inλ3≥0.

We divide (4.14) by the wave celerity (a positive quantity) and substitute the propagation speeds (4.1) to find

 (4.15) unc+12+γ2Inunc+θ2In(unc−12)≥0.

At an inflow boundary, we know that the normal velocity is negative or . Thus, we can rewrite (4.15) into the form

 (4.16) γ2InFr+θ2In(12+Fr)≤12−Fr.

The expression (4.16) defines an ellipse in the plane as shown in Fig. 1. In order to guarantee energy stability the values of and must lie within this ellipse for a given value of the Froude number.

We now know how to select energy stable values of and for inflow boundaries due to (4.16). Three special cases are:

1. , : The boundary conditions in this case satisfy the condition from Theorem 4.2 where . The tangential variable is a scaled value of the internal information at the boundary whereas the left traveling characteristic component, , is taken entirely from external data, similar to Mcdonald [20].

2. : The boundary conditions in this case satisfy the condition from Theorem 4.2 where . Here the tangential and left characteristic variables are taken entirely from boundary data. This is a similar treatment as proposed by Blayo and Debreu [2].

3. , : This boundary condition sets the tangential velocity to zero and either prescribes the normal velocity component () or the geopotential (). These boundary conditions fail to satisfy the condition (4.14) since

 λ1+γ2Inλ2+θ2Inλ3=λ1+λ3=2un<0.

However, the energy stability theory described herein only gives a sufficient condition for stability. Thus, these boundary conditions are not excluded out-right and, in fact, they are used in practice, e.g. [3, 35, 39]. Moreover, Oliger and Sundström [26], analyzed these boundary conditions using normal-mode analysis [13] and found that (specification of normal velocity) is stable, whereas (specification of water height) is not.

If we are in the regime where , then we need three boundary conditions also for the subcritical inflow. This corresponds to the supercritical inflow boundary conditions (4.12).

Next, we examine subcritical outflow where . First we consider flows with and must prescribe one boundary condition. Now there are two outgoing characteristic components, , and a single incoming one, . We take the coefficient matrix from (4.4) to be

 R=[γOutθOut]

and find the boundary term to have the form

 (4.17) w3=[γOutθOut][w1w2].

For subcritical outflow with the sufficient condition for nonlinear energy stability from Theorem 4.2 becomes

 (4.18) Λ++RTΛ−R=[λ1+γ2Outλ3γOutθOutλ3γOutθOutλ3λ2+θ2Outλ3]≥0.

Again, the constants and must be determined. The contribution (4.1) becomes

 (4.19) w21λ1+w22λ2+w23λ3 =w21λ1+w22λ2+(γOutw1+θOutw2)2λ3 =w21(λ1+γ2Outλ3)+w22(λ2+θ2Outλ3) +2γOutθOutλ3w1w2.

To guarantee energy stability, (4.19) must be positive for all values of and . This is true only if the expression is positive definite, i.e. the discriminant must be negative and the coefficients of the square terms positive. It turns out that the positivity of the square terms is included in the condition on the discriminant. Thus, it is sufficient to investigate

 (2γOutθOutλ3)2−4(λ1+γ2Outλ3)(λ2+θ2Outλ3)≤0.

We expand and rearrange terms to find

 −γ2Outλ2λ3−θ2Outλ1λ3≤λ1λ2.

Next, we divide by , substitute the form of the three propagation speeds, and rewrite the expression in terms of the Froude number (4.5)

 (4.20) γ2OutFr(12−Fr)+θ2Out(12+Fr)(12−Fr)≤Fr(12+Fr).

Shown in Fig. 2, (4.20) is an ellipse in the plane that defines the possible energy stable values to construct the subcritical outflow boundary condition (4.17).

Two important subsets of subcritical outflow boundary condition emerge:

1. : Here, the condition for energy stability becomes

 Λ++RTΛ−R=[λ1+λ300λ2]=[2un00un]>0.

This corresponds to a special choice of the boundary condition where either