**2. Numerical Method**

Consider the scalar conservation law

$$
\frac{
\partial
}{
\partial t
}\mu + \frac{
\partial
}{
\partial x
}f(\mu) = 0.
\tag{1}
$$

The weak solutions for Equation (1) are not unique and we are interested in the so-called entropy solution, which satisfies the entropy condition

$$
\frac{
\partial
}{
\partial t
}
\eta(u) + \frac{
\partial
}{
\partial x
}q(u) \le 0.
\tag{2}
$$

Here, the entropy function *η*(*u*) is convex and the entropy flux function *q*(*u*) satisfies *<sup>∂</sup><sup>q</sup> <sup>∂</sup><sup>u</sup>* <sup>=</sup> *<sup>v</sup> <sup>∂</sup> <sup>f</sup> <sup>∂</sup><sup>u</sup>* with the entropy variable *v* = *∂η <sup>∂</sup><sup>u</sup>* . To solve Equation (1) with the conservative difference method, we have the semi-discrete scheme

$$\frac{d\mu\_i}{dt} = -\frac{f\_{i+1/2} - f\_{i-1/2}}{\Delta x},\tag{3}$$

where *ui* represents the point valve at the node *xi* on a uniform Cartesian mesh with the mesh size *xi*+<sup>1</sup> − *xi* = Δ*x*. The scheme (3) is said to be entropy stable if it satisfies a discrete version of the entropy condition (2), namely,

$$\frac{\partial}{\partial t}\eta(u\_i) + \frac{q\_{i+1/2} - q\_{i-1/2}}{\Delta x} \le 0;\tag{4}$$

the scheme is said to be entropy conservative if it satisfies the discrete entropy equality

$$\frac{\partial}{\partial t}\eta(u\_i) + \frac{q\_{i+1/2} - q\_{i-1/2}}{\Delta x} = 0. \tag{5}$$

In this work, we utilize the familiar Runge–Kutta scheme of order four to solve the ordinary Equation (3). Now, we focus on how to build the entropy stable numerical flux *fi*+1/2.

First, the numerical flux is split into an entropy conservative part and a numerical diffusion part

$$f\_{i+1/2} = f\_{i+1/2}^{\text{EC}} - \frac{1}{2} \delta\_{i+1/2} \left( v\_{i+1/2}^+ - v\_{i+1/2}^- \right). \tag{6}$$

Here, *f EC <sup>i</sup>*+1/2 is a fourth order entropy conservative flux [15] expressed as

$$f\_{i+1/2}^{\text{EC}} = \frac{4}{3}\overline{f}\left(u\_i, u\_{i+1}\right) - \frac{1}{6}\left(\overline{f}\left(u\_{i-1}, u\_{i+1}\right) - \overline{f}\left(u\_i, u\_{i+2}\right)\right) \tag{7}$$

with ˜ *f*(*α*, *β*) being the second order entropy conservative flux defined by Tadmor [3]. Generally, *δi*+1/2 in the numerical diffusion part is chosen to be | *f* (*ui*<sup>+</sup>1/2)|. If the reconstructed values *v*<sup>±</sup> *<sup>i</sup>*+1/2 at the interfaces from the entropy variable *vi* satisfy

$$\text{sign}(v\_{i+1/2}^+ - v\_{i+1/2}^-) = \text{sign}(v\_{i+1} - v\_i),\tag{8}$$

the numerical flux (6) achieves entropy stability. Property (8) in which the jump of the reconstructed point values at each cell interface has the same sign as the jump of the underlying point values across that interface is called the sign property. We present the fourth order reconstruction satisfying the sign property as follows to obtain *v*± *<sup>i</sup>*+1/2.

Consider a polynomial reconstruction of the form

$$p\_i(\mathbf{x}) = v\_i + d\_i \left(\frac{\mathbf{x} - \mathbf{x}\_i}{\Delta \mathbf{x}}\right) + \left(\frac{v\_{i-1} - 2v\_i + v\_{i+1}}{2}\right) \left(\frac{\mathbf{x} - \mathbf{x}\_i}{\Delta \mathbf{x}}\right)^2 + \left(\frac{-v\_{i-1} + v\_{i+1} - 2d\_i}{2}\right) \left(\frac{\mathbf{x} - \mathbf{x}\_i}{\Delta \mathbf{x}}\right)^3 \tag{9}$$

with *di* being a slope function. For convenience, we introduce the following notation

$$\begin{aligned} ds\_i &= \frac{2}{3}v\_{i+1} - \frac{2}{3}v\_{i-1} - \frac{1}{12}v\_{i+2} + \frac{1}{12}v\_{i-2}, \\ \text{WC}\_i &= v\_{i+1} - v\_{i-1}, \quad \text{WR}\_i = v\_{i+1} - v\_i, \quad \text{WC}2\_i = v\_{i+2} - v\_{i-2}, \\ ds\_1 &= 0, \quad ds2\_i = \frac{1}{2} \left( \text{WC}\_i - 8 \text{WR}\_i \right), \quad ds3\_i = \frac{1}{2} \left( 8 \text{WR}\_i - 7 \text{WC}\_i \right), \\ S\_i &= \text{sign}(\text{WC}\_i), \\ \text{C}\_1 &= \frac{\sqrt{3}}{6}, \quad \text{C}\_2 = \frac{6}{12 + \sqrt{3}}, \end{aligned} \tag{10}$$

and define the slope *di* in the following way:


$$\text{(a)}\qquad\text{If }v\_i = \frac{v\_{i+1} + v\_{i-1}}{2}\text{, we define}$$

$$d\_i = \begin{cases} \max(ds1\_{i\prime}ds\_i)\_{\prime} & \text{if} \quad S\_i > 0\_{\prime} \\ \min(ds1\_{i\prime}ds\_i)\_{\prime} & \text{if} \quad S\_i < 0\_{\prime} \end{cases}$$

$$\text{(b)}\qquad\text{If }v\_i \neq \frac{v\_{i+1} + v\_{i-1}}{2}\text{, then the following hold:}$$

$$\text{i.t.} \quad \text{If } \left| \mathcal{WR}\_i - \frac{1}{2} \mathcal{WC}\_i \right| \ge \frac{1}{8} \left| \mathcal{WC2}\_i - 2 \mathcal{WC}\_i \right|, \text{ then}$$

$$d\_{i} = \begin{cases} \max(ds2\_{i}, ds3\_{\prime}ds\_{i})\_{\prime} & \text{if} \quad S\_{i} > 0\_{\prime}, \\\min(ds2\_{i}, ds3\_{\prime}ds\_{i})\_{\prime} & \text{if} \quad S\_{i} < 0. \end{cases}$$

$$\text{iii.} \qquad \text{If } \left| \mathcal{W} \mathcal{R}\_i - \frac{1}{2} \mathcal{W} \mathcal{C}\_i \right| < \frac{1}{8} \left| \mathcal{W} \mathcal{C} \mathcal{Z}\_i - 2 \mathcal{W} \mathcal{C}\_i \right|, \text{then}$$

$$d\_{i} = \begin{cases} \frac{W\mathbf{C}\_{i}}{2} - S\_{i} \cdot \mathbf{C}1 \cdot |2WR\_{i} - WC\_{i}| \, & \text{ if } \quad \left| \frac{WR\_{i}}{WC\_{i}} - \frac{1}{2} \right| \le C2, \\\frac{WC\_{i}}{2} \, & \text{ if } \quad \left| \frac{WR\_{i}}{WC\_{i}} - \frac{1}{2} \right| > C2. \end{cases}$$

**Theorem 1.** *With this definition of the slope, the polynomial reconstruction* (9) *fulfills the shape preserving property, namely,*


The proof of this theorem can be found in [16]. According to this property, the polynomial *pi*(*x*) does not create new extrema in the interior of [*xi*−1/2, *xi*<sup>+</sup>1/2]. In other words, the polynomial *pi*(*x*) may create new extrema at the interfaces {*xi*±1/2}. To keep away from such spurious extrema, *Entropy* **2019**, *21*, 508

we can employ a similar strategy as in [17]. Let us modify the polynomial reconstruction to be the following form:

$$
\varphi\_i(\mathbf{x}) = (1 - \theta\_i)\upsilon\_i + \theta\_i p\_i(\mathbf{x}), \quad 0 < \theta\_i < 1. \tag{11}
$$

The limiter *θ<sup>i</sup>* is determined by

$$\theta\_i = \begin{cases} \min\left\{\frac{M\_{i+1/2} - v\_i}{M\_i - v\_i}, \frac{m\_{i-1/2} - v\_i}{m\_i - v\_i}, 1\right\}, & \text{if} \quad v\_{i-1} < v\_i < v\_{i+1/2} \\ \min\left\{\frac{M\_{i-1/2} - v\_i}{M\_i - v\_i}, \frac{m\_{i+1/2} - v\_i}{m\_i - v\_i}, 1\right\}, & \text{if} \quad v\_{i-1} > v\_i > v\_{i+1/2} \\ 1, & \text{otherwise}, \end{cases} \tag{12}$$

with

$$M\_{\boldsymbol{i}} = \max\_{\boldsymbol{x} \in \left[\underline{\boldsymbol{x}}\_{i-1/2}, \underline{\boldsymbol{x}}\_{i+1/2}\right]} p\_{\boldsymbol{i}}(\boldsymbol{x}), \quad m\_{\boldsymbol{i}} = \min\_{\boldsymbol{x} \in \left[\underline{\boldsymbol{x}}\_{i-1/2}, \underline{\boldsymbol{x}}\_{i+1/2}\right]} p\_{\boldsymbol{i}}(\boldsymbol{x}).$$

and

$$\begin{cases} M\_{i \pm 1/2} = \max\left\{ \frac{1}{2} (v\_i + v\_{i \pm 1}), p\_{i \pm 1}(\chi\_{i \pm 1/2}) \right\}, \\ m\_{i \pm 1/2} = \min\left\{ \frac{1}{2} (v\_i + v\_{i \pm 1}), p\_{i \pm 1}(\chi\_{i \pm 1/2}) \right\}. \end{cases}$$

**Theorem 2.** *The polynomial reconstruction* (11) *is fourth order accurate and fulfills the sign property* (8)*.*

**Theorem 3.** *With the definitions of the reconstructed values v*<sup>+</sup> *<sup>i</sup>*+1/2 = *ϕi*+1(*xi*<sup>+</sup>1/2) *and v*<sup>−</sup> *<sup>i</sup>*+1/2 = *ϕi*(*xi*<sup>+</sup>1/2)*, the numerical flux* (6) *is entropy stable and fourth order accurate.*

**Remark 1.** *The proofs for Theorems 2 and 3 can be carried out very similarly as in [12,17]. We omit the simple but trivial procedure here.*

**Remark 2.** *For hyperbolic systems, the reconstruction procedure should be implemented in the local characteristic directions for the purpose of achieving entropy stability. The detailed steps can be found in our previous paper [12].*
