CONVERGENCE OF ENTROPY STABLE SCHEMES FOR DEGENERATE PARABOLIC EQUATIONS WITH A DISCONTINUOUS CONVECTION TERM

. Three-point entropy stable schemes are extended for partial differential equations of the degenerate convection-diffusion type where a discontinuous space-dependent function is incorporated in the convective flux. Using the compensated compactness theory, convergence of the proposed entropy stable approximations to the entropy weak solution is proved. Assuming the so-called potential condition in the jump discontinuities, an estimate for entropy functions is demonstrated. Finally, using benchmark tests a validation of the efficiency of the entropy stable scheme is provided by comparison with an upwind-type solution.


Introduction
In this work we focus on scalar degenerate parabolic-hyperbolic equations with discontinuous coefficients of the form   +  ((), )  = ()  , (, ) ∈ Π  := R × [0,  ], (1.1) where  = (, ) represents a physical property of a system (mass, velocity, energy, etc.),  is the transport or convective flux of the physical property ,  is a function modeling velocity variations dependent on space and  serves as the integrated diffusion coefficient of .We consider a degenerate diffusion, that is,  ′ () vanishes at one point or over an entire interval, which is commonly named strong degeneracy.Equations of the form (1.1) are used to model several phenomena such as flow in porous media [9,10,12], consolidation-sedimentation process in clarifier-thickener units [2,[4][5][6] and traffic flow in roads where the flux depends on the position or the driver's reaction [1].Notice that if there are discontinuities in the equation coefficients or a degenerate diffusion, classical solutions of equation (1.1) might not exist.For such reason, it is necessary to seek solutions in a variational sense which are known as weak solutions.
In this work, we are interested in proving the convergence of three-point entropy stable approximations to entropy solutions for the degenerate parabolic problem (1.1) and (1.2).An entropy stable framework for conservation laws was introduced by Tadmor in [25,27], for nonconservative hyperbolic systems was derived in [7], extended in [17] for degenerate parabolic conservation laws and high-order entropy stable essentially nonoscillatory schemes were proposed in [15].Entropy stable schemes proposed here is a family of space finite-difference semi-discretizations which are second-order accurate for smooth zones whereas the upwind discretization is firstorder accurate.This feature differentiates how both methods deal with discontinuities and smooth regions of the solution, as we will show later.For the convergence proof we will use the compensated compactness theory allowing to deal with sign changes of the discontinuous function  and the non convexity of the convective flux without further costs.This methodology has been used before to prove existence of weak solutions for degenerate parabolic equations [18] and convergence of numerical approximations for the scalar and multidimensional conservation laws with continuous and discontinuous coefficients [3,8,14,20].
( 1.5) The equivalence between the previous definition and Kružkov-type entropy solutions will be discussed in Appendix A.
In particular, Kružkov entropy function provides a 3-tuple with the previous desired symmetry property for the scalar case [21].
Let us now denote by  =  ′ () the entropy variable.The mapping  →  is one to one and by inverting such function, (), we define the potential function as follows: Let  and  be two entropy solutions of the Cauchy problem (1.1) and (1.2) satisfying the symmetry condition (1.6).The following two extra conditions are also imposed: Let us exemplify the previous assumptions for the particular case of the Kružkov entropy 3-tuple.For Assumption 1.3, its condition is equivalent to   ((), ) ′ () ≤ 0, ∀.
Therefore, if there exists a value  such that  →  ((), ) is constant, then Assumption 1.3 is always verified for the Kružkov 3-tuple.In particular,  ((), 0) = 0, ∀, is a physically expected condition.Now, Assumption 1.4 is proven for Krûzkov entropy 3-tuples in [21] if for any states  and  the crossing condition is assumed for any jump in , where ( − ,  + ) are the associated left and right limits.
The paper is organized as follows: In Section 2, we propose the entropy stable discretization for equation (1.1) which satisfies a discrete entropy inequality consistent with the entropy weak solution given in (1.2).In Section 3, we provide one of our main results: the  ∞ convergence of the entropy stable scheme.Which is proved based on the works of Karlsen et al. [18] and Fjordholm [14] using the compensated compactness methodology.Section 4 is devoted to prove the entropy estimate considering Assumptions 1.3 and 1.4.In Section 5 we illustrate the performance of the entropy stable approximation which captures correctly different dynamics like shock waves, rarefaction waves and oscillating solutions, improving in some cases the upwind approximation.Finally, conclusions are drawn in Section 6.

Entropy stable approximation
In this section, we adapt the entropy stable schemes proposed in [17] to the nonlinear parabolic equation (1.1).There are some previous works in the literature devoted to adapt numerical methods to treat the type of discontinuities in the advection term, for example upwind schemes [19,22,30], relaxation schemes [20] and Engquist-Osher schemes [5,29].In our case we consider some ideas from Karlsen et al.'s work [19] for the approximation of the function ().
Let us first discretize uniformly the spatial domain where the mesh points {  } ∈N are obtained by the stepsize ∆ =  +1 −   .To approximate the function  into the convection flux we consider a piecewise function defined as follows: for all  ∈ [  ,  +1 ).Here it is assumed that the possible discontinuities of  are placed at the cell interfaces.Approximation (2.1) is simple and has given satisfactory results [19].
We now present the entropy conservative scheme for equation (1.1) proposed in [17] where [] + 1 2 = ( +1 ) − (  ) and  + 1 2 is any numerical flux satisfying is the entropy variable and  is the potential function defined in (1.7).Moreover, we consider numerical fluxes of the form where  is an operator satisfying the Lipschitz condition: with   a positive constant depending on the operator .For instance, where Notice that the approximation of the integrated diffusion coefficient is valid for degenerate and nondegenerate diffusion terms.When the convection term dominates over the viscous term or the diffusive term vanishes (degenerate case), then some extra viscosity is required in order to get entropy stability [17,25].Thus, to prevent spurious oscillations, some numerical viscosity is added to the scheme (2.2) and the entropy stable solution is then given by the formula In the next two sections, we will prove results to determine the convergence of the proposed scheme to the entropy solution in the sense of Definition 1.2.

Discrete entropy inequality
Here we will demonstrate that scheme (2.2) satisfies a discrete entropy inequality consistent with entropy weak solutions given by (1.4).We start giving the following definition: Definition 2.1.Given an entropy 3-tuple (, , ), the numerical method (2.6) is said to be entropy stable for equation ( for some consistent entropy numerical fluxes,  + 1 2 and  +  2 , and with a function  that approximates the product   •  ′ . The numerical flux  + 1  2 given in (2.3) is consistent with the physical flux  (, ).Moreover, as a consequence of the facts that function  is increasing and  is convex ( ′ is increasing as well) then holds.Where   and   denote the left and right limits of  respectively.Following Theorem 3.2 given in [17], we will obtain that scheme (2.2) is an entropy stable approximation since it satisfies inequality (2.7).Let us first multiply scheme (2.6) to the left by the value   , obtaining (2.9) , we get and  =  −  and after some algebraic manipulations we obtain )︁ , and Using the previous identities, scheme (2.9) can be written as where where   =   is a positive real number) the right side of scheme (2.10) is negative and scheme (2.6) satisfies an inequality of the form (2.7).
After some algebraic manipulations, we can rewrite the previous scheme as follows:

2
. Regrouping terms, we get From condition (2.4) and Lipschitz continuity of  , we obtain Notice that for  = 0 and   = 1 (particular case of operator (2.5)), inequality (3.1) establishes a standard CFL condition for viscous conservation laws.Lemma 3.1 provides a generalization of such condition when numerical viscosity is introduced.
In the following result, we will obtain a bound for the diffusion increment.Lemma 3.3.Assume hypothesis (H1)-(H6) and consider  0 () ∈  ∞ (R) with compact support.Let  Δ be the entropy stable solution defined in (2.6) taking  = ∆ for  > 0, then there exists a constant  independent of ∆ and ∆ such that In particular, this implies that { Proof.Let us start finding a bound for the spatial variation of the diffusion coefficient.Multiplying by ∆∆ the entropy stable scheme (2.6) and discretizing the time derivative, we get ∆∆ We may consider when ∆ is sufficiently small.Adding now over the discretized domain Π  and using in the first term the formula We now use summation by parts and assume that  Δ vanishes sufficiently fast when  → ∞ and Thus, the term of the right hand side of the above equation can be bounded as follows: (3.4) The last term of the right hand side in the previous inequality is bounded since  is a Lipschitz function from definition (2.3) and  Δ is bounded.Discretizing the identity (  (())) Let us now bound the  2 -norm of the time variation of the diffusion coefficient over the domain Π  −Δ where we have used the summation by parts formula, Hölder's inequality, the Lipschitz continuity of the functions  and , inequality (3.5) and  = ∆.Thus, condition (3.3) is proved with the spatial and time variation bound obtained.Furthermore, by the Kolmogorov compactness criterion [16], it is concluded that { To complete the study of the compactness of the diffusion function, we state the fundamental theorem of compensated compactness theory.Furthermore, for any continuous function Φ :  → R, we have along a subsequence The function (, ) is known as a Young measure and Prob(R  ) is the set of all probability measures in the Borel sets.As a consequence of this theorem, a uniformly bounded sequence { Δ } Δ>0 converges to  a.e. on Π  if and only if  (,) =  (,) where   is a Dirac measure located at (, ).Finally, we are ready to prove the compactness property of the diffusive part.
Theorem 3.5.Assume hypothesis (H1)-(H6) and consider  0 () ∈  ∞ (R) with compact support.Let  Δ be a sequence of entropy stable approximations for problem (1.1) and (1.2).Then, a subsequence of { where  is the weak−⋆ limit of entropy stable solution Proof.It is straightforward from Theorems 3.1, 3.3 and 3.4.Now we study the compactness of {  +   }.For that, we extend Lemma 3.2 in Fjordholm [14] for parabolic conservation laws.Let us define {  +   } and {  +   −   + } as measures and let  be a borelian in R, then where  Δ is the numerical solution considered as a piecewise constant function.Let Ω ⊂ Π  be an open and bounded set with supp (︀  Δ )︀ ⊂ Ω for all ∆ > 0. For a function  ∈  −1,∞ (Ω) with  =0, = 0 we can define {  +   } and {  +   −   + } as functionals of the form The following result is known as the Murat's lemma and it will be useful to prove the compactness in  −1 loc (Ω) of the sequence {  +   }.
Let Ω ⊂ Π  be a bounded open set, (, ) an entropy pair and  an entropy numerical flux that satisfies the Lipschitz condition Assume the sequence of functions  Δ with ∆ > 0 satisfies: ⊂ Ω, for all ∆x > 0 (3.9)Then, the sequence Proof.Decomposing the functional as , and using inequality (3.11), we get that  1 () is bounded in ℳ(Ω) since .
We prove now that the entropy stable scheme defined in (2.6) satisfies conditions (3.7)-(3.11)from Lemma 3.7: (i) Substituting  + 1  2 and considering  =  + 1 2 and  = ∆ for some fixed  > 0, we get . The previous inequalities were obtained since  ′ and  ′′ are bounded in a compact set and  and  are Lipschitz continuous.Thus, condition (3.7) is satisfied for scheme (2.6).
(ii) By Lemma 3.1 and taking into account that the initial condition  0 is a compact support function, the approximation sequence  Δ satisfies conditions (3.8), (3.9) and Thus, we also have the condition (3.10).(iii) Using the definition of the entropy stable scheme given in (2.10), we get We now demonstrate that the succession of entropy stable approximations has as limit an entropy weak solution which satisfies Definition 1.2.Proof.Let us multiply the entropy stable scheme (2.6) by ∆∆   as follows where    is a compact support function verifying   0 = 0 with  the time index and    = 0 for all |  |  ∆.Approximating the partial derivative of time of (3.12) and using the summation by parts formula, we get Since  0 () is a compact support function, consequently, the approximation  Δ is also a compact support function for any .Thus, the previous formula approaches the weak solution (1.3) as ∆, ∆ → 0 since  = ∆.
where    is a compact support function.Considering the addition by parts in the previous inequality, we get the following: The second integral in this previous inequality can be written as the sum of two integrals.One that considers the continuous part of  and the other that includes its discontinuities.Thus, obtaining the entropy inequality (1.4).
Finally, we are already to prove the  ∞ -convergence of the entropy stable scheme.For that, we state the following result which is a direct extension of Lemma 2.2 given in [18], where it is proven for the case where  is multiplicative, so we omit its proof here.The proof uses the main theorem of compensated compactness (Thm.3.4).
We can now state our main convergence theorem.Theorem 3.11.Let (, , ) be an entropy 3-tuple associated to equation (1.1).Assume that  0 () ∈  ∞ (R) with compact support.Assume also that hypothesis (H1)-(H6) and conditions (3.1) and (3.7)-(3.11)hold.Considering that the discrete equation (2.10) conforms an entropy stable scheme such that  = ∆ for some  > 0.Then, there exists a weak solution  of the parabolic problem (1.1) that can be constructed as the limit of the approximations  Δ given by the entropy stable scheme (2.6) and this structure guarantees that  is also the entropy solution.
Remark 3.12.The previous convergence proof has been made considering the fully discrete numerical method obtained by applying forward Euler for the time stepping and entropy stable discretization for the spatial approximation.Thus, the resulting method is the first-order of accuracy in time and second-order in space.
Let us consider the change of variables Then, we can rewrite and inequality (4.5) can be written as where  Using the previous theorem, we will prove the a priori estimate mentioned before.)︀ .Now our goal is to find a stability result that can be guaranteed when the test function does not vanish around the discontinuities of .

Numerical examples
Here the entropy stable approximation (2.6) is applied to several degenerate parabolic tests with discontinuous coefficients in order to illustrate its performance.Taking into account that the upwind scheme has been proved convergent for such problems, it is also interesting to compare both numerical solutions.To obtain the numerical solutions, we consider a semi-discretization in space of equation (1.1) as follows: which is solved using the TVD-RK2 (1) =   + ∆ Δ (  ), )︁ .
The operator  Δ in each ( −1 ,  −1 ) is defined as in the case of the entropy stable method, being   +1/2 the numerical entropy flux defined in (2.3) and in the case of the upwind scheme, where It is important to point out that the proposed TVD-RK2 scheme is stable under the same CFL condition that the explicit Euler time-discretization [24].Therefore, the stability condition for our numerical scheme is (3.1).

Viscous Burgers' equation with continuous coefficients
We start analyzing the error and convergence rate of the proposed entropy stable (ES) scheme.For that, we consider a viscous Burgers' equation with a continuous coefficient () and an integrated diffusion coefficient without degeneracy along with the initial condition  0 () = (−2cos())/(2 + sin()). (5.2) In Table 1 it is shown the results of the numerical convergence by considering the L ∞ -norm.The error is obtained using an upwind reference solution with a spatial mesh of 1280 nodes.Notice that the convergence order of the entropy stable (ES) method is initially greater than three and at the end gets closer to two.In Figure 1, a comparison between the entropy stable approximation and the reference solution for problem (5.1) and (5.2) is presented.

Linear equation with oscillating solution
Let us consider the linear advection equation with initial condition For such Cauchy problem, the exact solution is known and it exhibits a fast oscillatory behavior, see Figure 2.
We show the convergence rates and error of the upwind (UP) and entropy stable methods for problem (5.3) and (5.4) in Table 2.We observe that the ES scheme is of second order accuracy while the UP method has difficulties in reaching the first order, which is related with the oscillating dynamics of the solution.
Based on the linear equation ( 5.3), we now consider the following degenerate parabolic equation: with the same initial condition, where To formulate the entropy stable approximation, we take the entropy function () =  2 /2 and from definition (2.3), we obtain the numerical flux In Figure 2, we can observe the numerical solutions for the hyperbolic and degenerate parabolic problems.It is interesting to notice that the dispersion of the upwind method for the hyperbolic case is maintained in the degenerate parabolic case.In Figure 3, we show the performance of both methods for two different advective fluxes considering two piecewise functions,  1 () and  2 ().It is important to note that the UP method requires further refinement of the mesh to leave behind the inherent scatter error that still appears in Figure 3. On the other hand, the stability of the ES approximation is regulated by the viscosity parameter , which must be higher as the degenerate diffusion and discontinuous coefficients are incorporated.

Burgers' equation with a rarefaction solution
Let us consider a well-known nonlinear hyperbolic equation  with initial condition The entropy solution of this problem is a continuous solution named as a rarefaction wave.It is well known that a standard upwind method develops a non-existent discontinuous jump, see Figure 4. Notice that the entropy stable approximation does not develop such discontinuity.We also propose a degenerate parabolic case associated to the previous equation where the integrated diffusion coefficient is (5.9) and The entropy stable approximation is obtained by considering () =  2 /2 and

𝑖+1
)︀ . (5.11) In Figure 4, it is shown how the wrong shock produced by the UP method is maintained despite the inclusion of a diffusion term in the equation.Moreover, it is important to point out that the ES approximation converges with CFL≤ 0.6, but the upwind approximation introduces several shocks for CFL values higher than 0.1.In Figure 5, we can observe the numerical solutions when a discontinuous flux, (), is considered.Both methods (ES and UP) reproduce well the jump of the solution as a consequence of the discontinuous flux.However, there is a small gap in the upwind solution caused by the wrong shock which is not eliminated even with a greater refinement in the mesh.For this problem, a rarefaction wave and a jump discontinuity connect two states of the solution,  = 0.6 and  ≈ 0.14, see Figure 6.The entropy stable method is constructed considering () = (log() − 1), obtaining

𝑖+1
)︀ +   )︀ .(5.14) In this example, we focus on studying the behavior of the UP and ES approximations when a mesh refinement is applied.In Figure 6, numerical simulations for the hyperbolic problem are shown.It is observed that the ES method needs further refinement to achieve the correct spatial position of the jump discontinuity.
For the degenerate parabolic problem, the refinement of the mesh does not offer a significant difference in the numerical solutions, see Figure 7. Another interesting feature is that the ES method does not require extra viscosity for the degenerate parabolic problem and the  parameter is zero, see Remark 3.2.

Conclusions
In this work we proposed three-point entropy stable finite-difference schemes for degenerate parabolic equations with a discontinuous convection term.The proposed entropy stable discretization is based on the entropy stable formulation given initially by Tadmor for conservation laws.Using the arithmetical mean as an estimation in each interval we included in a simple way the approximation of the discontinuous coefficient.The entropy stable solution was proven to converge to entropy weak solutions via the compensated compactness methodology.Such methodology allowed us to consider discontinuities and sign changes of the discontinuous coefficient.The entropy stable sequence has a free parameter, for which we found a theoretical bound; in the past, this parameter was found by trial and error.
To validate and analyze the efficiency of the proposed three-point entropy stable methods, numerical simulations were provided for benchmark problems.We conclude that the spatial entropy stable discretization has the advantages of a second-order accuracy method with the added bonus of capturing shocks correctly.Thus, the entropy stable solution can be an alternative to the upwind solution as a reference solution to solve degenerate convection-diffusion problems and to validate similar numerical methods.
In the near future, we would like to construct high-order entropy stable solutions for degenerate parabolic equations and extend them to solve multidimensional problems.For applications, it is important to search for entropy functions that perform well in several situations.
-Assume that a weak solution  satisfies the entropy 3-tuples inequality (1.4).To prove that  also satisfies the Kružkov inequality (A.1), we suppose conditions (H1)-(H5) and use similar arguments as given in Burger et al.By the Kružkov type entropy inequality (A.1) and (1.3), (1.4) is valid for   ,   and   .Then, since   ,   and   converge uniformly to ,  and  respectively on [a, b] we can interchange the limit with the integral.Therefore, the entropy inequality (1.4) is valid for the given , which was consider an arbitrary entropy function.

Theorem 3 . 8 .
Assume  0 () ∈  ∞ (R) with compact support.If the entropy stable scheme  Δ proposed in (2.6) taking  = ∆ ( > 0) has a limit, then it is a weak solution to the problem (1.1) and (1.2) in the sense of Definition 1.1.

Figure 5 .
Figure 5. Numerical solution for the parabolic problem (5.5) and (5.7) with () given in (5.10) (left) at time  = 0.35 taking ∆ = 0.01 and CFL = 0.9.On the right hand side, the total entropy decay rate is shown.We used  = 0.52 for both graphs.

Table 1 .
Error and convergence rate for the entropy stable scheme for problem (5.1) and (5.2) at time  = 0.15.

Table 2 .
Error and convergence rate for the entropy stable and upwind schemes for problem (5.3) and (5.4) at time  = 2.