Processing math: 42%
Research article Special Issues

Homogenization of composite materials reinforced with unidirectional fibres with complex curvilinear cross section: a virtual element approach

  • Received: 07 December 2023 Revised: 03 June 2024 Accepted: 13 July 2024 Published: 29 July 2024
  • The paper presents an augmented curvilinear virtual element method to determine homogenized in-plane shear material moduli of long-fibre reinforced composites in the framework of asymptotic homogenization method. The new virtual element combine an exact representation of the curvilinear computational geometry for complex fibre cross section shapes through an innovative two-dimensional cubature suite for NURBS-like polygonal domains. A selection of representative numerical tests supports the accuracy and efficiency of the proposed approach for both doubly periodic and random fibre arrangement with matrix domain.

    Citation: E. Artioli, G. Elefante, A. Sommariva, M. Vianello. Homogenization of composite materials reinforced with unidirectional fibres with complex curvilinear cross section: a virtual element approach[J]. Mathematics in Engineering, 2024, 6(4): 510-535. doi: 10.3934/mine.2024021

    Related Papers:

    [1] Raimund Bürger, Kenneth H. Karlsen, John D. Towers . On some difference schemes and entropy conditions for a class of multi-species kinematic flow models with discontinuous flux. Networks and Heterogeneous Media, 2010, 5(3): 461-485. doi: 10.3934/nhm.2010.5.461
    [2] Raimund Bürger, Christophe Chalons, Rafael Ordoñez, Luis Miguel Villada . A multiclass Lighthill-Whitham-Richards traffic model with a discontinuous velocity function. Networks and Heterogeneous Media, 2021, 16(2): 187-219. doi: 10.3934/nhm.2021004
    [3] Raimund Bürger, Stefan Diehl, María Carmen Martí . A conservation law with multiply discontinuous flux modelling a flotation column. Networks and Heterogeneous Media, 2018, 13(2): 339-371. doi: 10.3934/nhm.2018015
    [4] Mauro Garavello, Roberto Natalini, Benedetto Piccoli, Andrea Terracina . Conservation laws with discontinuous flux. Networks and Heterogeneous Media, 2007, 2(1): 159-179. doi: 10.3934/nhm.2007.2.159
    [5] Giuseppe Maria Coclite, Lorenzo di Ruvo, Jan Ernest, Siddhartha Mishra . Convergence of vanishing capillarity approximations for scalar conservation laws with discontinuous fluxes. Networks and Heterogeneous Media, 2013, 8(4): 969-984. doi: 10.3934/nhm.2013.8.969
    [6] Maya Briani, Emiliano Cristiani . An easy-to-use algorithm for simulating traffic flow on networks: Theoretical study. Networks and Heterogeneous Media, 2014, 9(3): 519-552. doi: 10.3934/nhm.2014.9.519
    [7] Christophe Chalons, Paola Goatin, Nicolas Seguin . General constrained conservation laws. Application to pedestrian flow modeling. Networks and Heterogeneous Media, 2013, 8(2): 433-463. doi: 10.3934/nhm.2013.8.433
    [8] Clément Cancès . On the effects of discontinuous capillarities for immiscible two-phase flows in porous media made of several rock-types. Networks and Heterogeneous Media, 2010, 5(3): 635-647. doi: 10.3934/nhm.2010.5.635
    [9] Wen Shen . Traveling wave profiles for a Follow-the-Leader model for traffic flow with rough road condition. Networks and Heterogeneous Media, 2018, 13(3): 449-478. doi: 10.3934/nhm.2018020
    [10] Adriano Festa, Simone Göttlich, Marion Pfirsching . A model for a network of conveyor belts with discontinuous speed and capacity. Networks and Heterogeneous Media, 2019, 14(2): 389-410. doi: 10.3934/nhm.2019016
  • The paper presents an augmented curvilinear virtual element method to determine homogenized in-plane shear material moduli of long-fibre reinforced composites in the framework of asymptotic homogenization method. The new virtual element combine an exact representation of the curvilinear computational geometry for complex fibre cross section shapes through an innovative two-dimensional cubature suite for NURBS-like polygonal domains. A selection of representative numerical tests supports the accuracy and efficiency of the proposed approach for both doubly periodic and random fibre arrangement with matrix domain.



    It is the purpose of this work to introduce, and in part analyze, a numerical scheme for a system of conservation laws with source terms of the type

    t(A(z)(ϕψ))+z(A(z)(J(ϕ,z,t)˜F(ϕ,ψ,z,t)))=Kk=1QF,k(t)(ϕF,k(t)ψF,k(t))δ(zzF,k), (1.1)

    where t is time, z is spatial position, and ϕ and ψ are the volume fractions of the primary and secondary disperse phases, respectively. Both disperse phases move within the continuous phase of the one-dimensional flow. We let A(z) denote a variable cross-sectional area. The flux functions J and ˜F are discontinuous across the positions z=zU<zF,1<<zF,K<zE, and due to constitutive assumptions of the model, are nonlinear functions of ϕ and ψ. The right-hand side of Eq (1.1) describes singular sources located at z=zF,k, k=1,,K, and is composed of given functions. It is assumed that QF,k(t) is the volumetric bulk flow of the mixture (of the continuous and two disperse phases) injected at z=zF,k, and that ϕF,k(t) and ψF,k(t) are the volume fractions of the primary and secondary disperse phases in the feed flow, respectively. The system Eq (1.1) is posed on ΠT:=R×(0,T) together with initial conditions

    ϕ(z,0)=ϕ0(z)for all zR, (1.2a)
    ψ(z,0)=ψ0(z)for all zR, (1.2b)

    where we assume that

    0ϕ0(z)1,0ψ0(z)1ϕ0(z)for allzR (1.3)

    along with

    TV(ϕ0)<,TV(ψ0)<. (1.4)

    Likewise, we assume that ϕF,k and ψF,k are piecewise continuous functions of bounded variation with a finite number of discontinuities and that they satisfy the bounds

    0ϕF,k(t)1,0ψF,k(t)1ϕF,k(t)for allk=1,,Kandt[0,T]. (1.5)

    (In later parts of the analysis we will assume that these functions and the bulk flows QF,k are constants.) If θ denotes the volume fraction of the continuous phase, then we assume that

    0ϕ,ψ,θ1;ϕ+ψ+θ=1, (1.6)

    which motivates assumptions Eq (1.3) and Eq (1.5). (Of course, satisfaction of Eq (1.6) by exact or numerical solutions of Eq (1.1), Eq (1.2) on ΠT needs to be proved.)

    A specific application that gives rise to the system Eq (1.1) is a model of a flotation column [8,9], where ϕ denotes the volume fraction of bubbles and ψ that of solid particles (Figure 1). The bottom of the column has the coordinate zU (the underflow) and the top zE (the effluent). The primary disperse phase of bubbles – specifically, aggregate bubbles, to which hydrophobic valuable particles (minerals) are attached – is assumed to flow through the suspension of solid particles and liquid independently of the volume fraction of solids. The secondary disperse phase consists of solid hydrophilic particles (ore) that move in the remaining space outside the bubbles. If the solid particles of the secondary disperse phase have a density larger than that of the fluid, the two disperse phases undergo counter-current, and otherwise, co-current flow. The distinction between primary and secondary disperse phase also becomes evident in the flux functions: the flux J of the primary disperse phase depends on ϕ only (besides z and t), while that of the secondary disperse phase, ˜F, depends both on ϕ and ψ. Thus, the system Eq (1.1) is triangular; however, it is generally non-strictly hyperbolic; see [9], where a counter-current model of the form Eq (1.1) is studied.

    Figure 1.  Schematic of a one-dimensional column with K=3 inlets and K+1=4 zones, where QU is the downwards volumetric outflow, QF,j is the volumetric flow at the inlet zF,j, for each j=1,,K, and QE is the upwards volumetric outflow. Note that the distances between the inlets/outlets are arbitrary and the cross-sectional area A=A(z) may vary piecewise continuously (although the figure shows a piecewise constant example).

    The main contribution of this work is an easily implemented explicit monotone numerical scheme for Eq (1.1). To properly address the theoretical support we are able to provide for this scheme, we refer to the complete model Eq (1.1), with all assumptions stated so far in effect, as "Model 1", and the corresponding scheme (that handles Model 1) as "Scheme 1". Additional properties of Scheme 1 can be established for two successively simplified versions of Model 1, named "Model 2" and "Model 3", for which the corresponding versions of Scheme 1 are addressed as "Scheme 2" and "Scheme 3", respectively.

    The scheme is supported by three partial theoretical arguments. Firstly, it is proved that Scheme 1 satisfies an invariant-region property, i.e., the approximate volume fractions satisfy a discrete analogue of Eq (1.6) at every point. Secondly, the assumption of a constant cross-sectional area, i.e.,

    Aconstant,A>0, (1.7)

    and time-independent feed and volume rates defines Model 2, and it is shown that the corresponding scheme for the primary disperse phase (the "ϕ-scheme" of Scheme 2; in short, "ϕ-Scheme 2") converges to a suitably defined entropy solution. Thirdly, we additionally assume that there are no flux discontinuities, so that Model 2 is reduced to the triangular system of conservation laws ("Model 3")

    tϕ+zJ(ϕ)=0, (1.8a)
    tψ+z˜F(ϕ,ψ)=0,(z,t)ΠT, (1.8b)

    where J and ˜F are z- and t-independent versions of the fluxes arising in Eq (1.1) and Eq (1.8) is equipped with the initial conditions Eq (1.2), where assumptions Eq (1.3) remain in effect. The corresponding reduced version of Scheme 2 that handles Model 3 is called "Scheme 3." Under these additional assumptions, we may invoke arguments of compensated compactness to prove that the scheme for the secondary disperse phase (the "ψ-Scheme 3") converges to a weak solution of the corresponding conservation law Eq (1.8b). Summarizing all arguments, we prove that Scheme 3 converges to a weak solution of the system Eq (1.8) in the sense of the following definition.

    Definition 1.1. The pair (ϕ,ψ) is called a weak solution of the initial value problem Eq (1.8), Eq (1.2) if

    (ⅰ) The functions ϕ and ψ belong to L(ΠT).

    (ⅱ) The functions ϕ and ψ satisfy Eq (1.8), Eq (1.2) in the sense of distributions on ΠT, that is, for each smooth test function ζ with compact support in ΠT, the following identities hold:

    ΠT(ϕtζ+J(ϕ)zζ)dzdt+Rϕ0(z)dz=0, (1.9)
    ΠT(ψtζ+˜F(ϕ,ψ)zζ)dzdt+Rψ0(z)dz=0. (1.10)

    (ⅲ) The function ϕ is an entropy solution of the single conservation law Eq (1.8a), that is, for each smooth and nonnegative test function ζ with compact support in ΠT, the following inequality holds for all cR:

    ΠT(|ϕc|tζ+sgn(ϕc)(J(ϕ)J(k))zζ)dzdt+R|ϕ0(z)c|dz0. (1.11)

    Numerical experiments illustrate that Scheme 1 for the full model Eq (1.1) (Model 1), approximates expected solution behaviour for counter-current and co-current flows and that approximate numerical errors tend to zero as the mesh is refined.

    The system Eq (1.1) models the evolution of the primary unknown ϕ independently of the secondary unknown ψ. One application of such triangular systems is the process of column flotation, which is a solid-liquid separation process used in mineral processing, environmental and chemical engineering [10,11,27,28,42,45]. The model Eq (1.1) restricted to three-phase counter-current flow in a flotation column was originally proposed in [9]. Its nonlinear constitutive assumptions come from the drift-flux theory (used to analyze the bubbly and froth regions [41,50,51]) and the solids-flux theory (for particles in a liquid [24,25,38]). In [9], the construction of steady-state solutions is detailed, where conservation laws with discontinuous flux are a key ingredient with a specific entropy condition [2,22,29]. The most interesting desired steady states are classified in [9] and visualized in graphical so-called operating charts that show how the control variables QU, QF:=QF,1 (feed mixture of gas, solids and water) and QW:=QF2 (feed washwater) should be chosen.

    The mathematical and numerical difficulties associated with Eq (1.1) are twofold; namely, one has to deal with discontinuities of the fluxes with respect to z, as well as with the definition of the governing model by a (triangular) system of conservation laws (in contrast to otherwise similar, known scalar two-phase models arising in flotation or sedimentation [8,13,22,23]). The well-known difficulty of conservation laws with discontinuous flux lies in the appropriate formulation of admissibility conditions of jumps of the solution across discontinuities of the flux such that the resulting concept of weak (discontinuous) solutions supplied with an entropy condition would admit a uniqueness result. There exist many criteria for selecting unique solutions (e.g., [1,22]), each of which corresponds to a particular physical reality and relies on specific assumptions on the fluxes adjacent to a discontinuity. A unified treatment of this problem is advanced in [2]. While scalar conservation laws with discontinuous flux have been studied widely, only a few analyses of systems with discontinuous flux are available (e.g., [16,47]). That said, its triangular nature makes Eq (1.1) potentially easier to treat than a full 2×2 system of conservation laws (where the flux of each component would depend on both unknowns).

    The triangular system with discontinuous flux studied in [9] was solved numerically with a staggered-grid scheme that utilizes the triangular structure of Eq (1.1). Such a semi-Godunov scheme for general triangular hyperbolic systems is one of the two suggested schemes by Karlsen et al. [32,33], who proved convergence of the numerical solutions under certain assumptions on the flux functions. We here propose a simpler numerical scheme (on a single grid) that is easier to implement and analyze. The analysis (of the scheme proposed under simplifying assumptions) relies on the aligned version of the scheme introduced in [33] and in particular on the convergence analysis of an Engquist-Osher scheme for multi-dimensional triangular system of conservation laws by Coclite et al. [18]. (The proof of convergence done in [18] is motivated by the more easily proven convergence of a vanishing viscosity approximation for the same model, see [17].) These analyses, and the present treatment for the reduced model Eq (1.8), rely on compactness techniques that use discrete entropy inequalities and the compensated compactness framework.

    Further applications and results on the analysis of triangular systems include two-component chromatography [3]. Furthermore, polymer flooding in oil recovery is modelled by a 2×2 system [31], which can be converted to a triangular system in Lagrange coordinates [49]. In [20,40,48], the authors study the delta shock wave formation in solutions of triangular system of conservation laws from the so-called generalized pressureless gas dynamics model. Bressan et al. [5] established the existence and uniqueness of vanishing viscosity solutions for scalar conservation laws for a Cauchy problem and their results can be applied to a triangular system under suitable assumptions. The results of Karlsen et al. [32,33] for general triangular systems can be applied to models of three-phase flows in porous media, for example, in oil recovery.

    The remainder of this paper is organized as follows. In Section 2, the model of [9] of gas-solid-liquid three-phase flow in a flotation column from is written in a slightly more general form. Starting from the balance equations of the three phases we outline the derivation of the algebraic forms of the fluxes J(ϕ,z,t) and ˜F(ϕ,ψ,z,t) arising in the governing PDE system Eq (1.1). In Section 3 the numerical method (Scheme 1) proposed for the approximation of solutions to the initial value problem Eq (1.1), Eq (1.2) (Model 1) is detailed, where computational effort is essentially reduced to the interior of the vessel (cf. Figure 1). After outlining the discretization in Section 3.1, we specify the numerical fluxes and update formulas for the primary and secondary disperse phases in Sections 3.2 and 3.3, respectively. Both formulas are adapted to the particular algebraic form of the fluxes J(ϕ,z,t) and ˜F(ϕ,ψ,z,t) and involve upwind discretizations, a particular monotone discretization for "concentration times velocity" fluxes from [6], and the Engquist-Osher numerical flux [26]. We then prove in Section 3.4 that Scheme 1 is monotone and that the numerical solutions satisfy a so-called invariant-region property (Theorems 3.1 and 3.2), that is, a discrete analogue of Eq (1.6), provided, of course, that the initial data satisfy Eq (1.3) and the time step and spatial meshwidth satisfy a CFL condition. Section 4 provides further partial results of the convergence analysis of the numerical scheme based on additional simplifying assumptions, namely those of a constant cross-sectional area A and constant bulk and feed flows QU, QF,k, ϕF,k and ψF,k (k=1,,K). Thus, the scheme under discussion is Scheme 2. We can then prove convergence of the ϕ-Scheme 2 (the one that discretizes the ϕ-component of the governing PDE; Section 4.1) and L1 Lipschitz continuity of the ψ-Scheme 2 (Section 4.2). If in addition all z-dependent flux discontinuities are removed, so that the governing PDE system is Eq (1.8) (Model 3) and the scheme reduces to Scheme 3, we may apply compensated compactness techniques to prove convergence of the ψ-scheme (Section 5). For the simplified problems, the initial conditions Eq (1.2) and assumptions Eq (1.3) and Eq (1.5) are imposed, so Theorems 3.1 and 3.2 remain in effect. While in that case the convergence of (the monotone) ϕ-Scheme 3 to an entropy solution of Eq (1.8a) follows by standard arguments (for monotone schemes), the principal result of Section 5 is convergence of ψ-Scheme 3 to a weak solution of Eq (1.8b) (Lemma 5.5 and Theorem 5.1). Estimations of errors and convergence order of the numerical method can be found in Section 6.1. Some numerical examples are presented in Section 6, starting with preliminaries (Section 6.1). First, in Section 6.3, we use a smooth solution to estimate the order of convergence. Later on, we present two numerical examples that illustrate the model predictions for counter-current (Section 6.4) and co-current flows (Sections 6.5 and 6.6). Finally, some conclusions are drawn in Section 7, and Appendix A contains the proofs of Lemma 5.5 and Theorem 5.1.

    The density of each phase is assumed constant, so the conservation of mass can be expressed by the balance equations (vϕ, vψ, and vθ are the phase velocities)

    t(A(z)ϕ)+z(A(z)ϕvϕ)=Kk=1QF,k(t)ϕF,k(t)δ(zzF,k), (2.1)
    t(A(z)ψ)+z(A(z)ψvψ)=Kk=1QF,k(t)ψF,k(t)δ(zzF,k), (2.2)
    t(A(z)θ)+z(A(z)θvθ)=Kk=1QF,k(t)(1ϕF,k(t)ψF,k(t))δ(zzF,k), (2.3)

    where the right-hand sides contain Dirac symbols, the feed volume fractions ϕF,k and ψF,k of the disperse phases, and the corresponding volume fraction 1ϕF,k(t)ψF,k(t) of the continuous phase, at the inlet located at z=zF,k, k=1,,K.

    We define the volume-average velocity, or bulk velocity, of the mixture by

    q:=ϕvϕ+ψvψ+θvθ,

    and replace Eq (2.3) by the sum of Eqs (2.1)–(2.3), which is

    z(A(z)q)=Kk=1QF,k(t)δ(zzF,k), (2.4)

    hence q varies with z due to the K inlets and the variable cross-sectional area. We define Q(z,t):=A(z)q(z,t) and integrate Eq (2.4) from any point z0<zU to obtain

    Q(z,t)=Q(z0,t)+Kk=1QF,k(t)H(zzF,k),

    where H() is the Heaviside function. If the volumetric underflow QU(t) is given, then Q(z,t)=QU(t) for z<zU, and

    Q(z,t)=QU(t)+Kk=1QF,k(t)H(zzF,k)=QU(t)+k:zF,kzQF,k(t).

    This continuity equation of the mixture replaces Eq (2.3). Next, Eq (2.1) and Eq (2.2) are rewritten in terms of q and two constitutive functions. We refer to the continuous phase and the secondary disperse phase as "secondary mixture", and define the volume fraction of the secondary disperse phase within the secondary mixture as

    φ:=ψψ+θ=ψ1ϕ(when ϕ<1),

    where 0φ1 by Eq (1.6). The volume-average velocity of the secondary mixture is

    qs:=ψvψ+θvθψ+θ=φvψ+1ϕψ1ϕvθ=φvψ+(1φ)vθ.

    It is then assumed that within [zU,zE), the relative velocity vϕs:=vϕqs of the primary disperse phase with respect to the secondary mixture is a given constitutive function ˜vϕs(ϕ), while outside that interval, both phases move at the same velocity, so their velocity difference is zero. Thus, in terms of the characteristic function

    γ(z):=χ[zU,zE)(z):={1for z[zU,zE),0for z[zU,zE),

    this assumption can be expressed as vϕs=γ(z)˜vϕs(ϕ). Within [zU,zE), the relative velocity of the secondary disperse phase with respect to the continuous phase vψθ:=vψvθ is supposed to be a given function ˜vψθ of φ, that is, vψθ=γ(z)˜vψθ(φ).

    The definitions of all velocities imply the identities

    ϕvϕ=ϕq+γ(z)ϕ(1ϕ)˜vϕs(ϕ),ψvψ=ψq+γ(z)ψ((1φ)˜vψθ(φ)ϕ˜vϕs(ϕ)) (2.5)

    for the (unweighted) fluxes ϕvϕ and ψvψ arising in Eq (2.1) and Eq (2.2), respectively. It is then useful to introduce the velocity and flux functions

    W(ϕ):=(1ϕ)˜vϕs(ϕ),V(φ):=σ(1φ)˜vψθ(φ),j(ϕ):=ϕW(ϕ),f(φ):=φV(φ), (2.6)

    where σ=±1 is chosen depending on the application such that V(φ),f(φ)0 (for standard convenience, e.g., when plotting their graphs); σ=1 for co-current flows (upwards) and σ=1 for counter-current flows. The velocity and flux of the secondary disperse phase with respect to z are therefore σV(φ) and σf(φ), respectively. We assume that W,V0 and V(1)=W(1)=0, as well as that

    fhas one local maximumωand one inflection point˜ω,0<ω<˜ω<1. (2.7)

    Combining Eq (2.5) and Eq (2.6) we obtain the expressions

    ϕvϕ=ϕq+γ(z)ϕW(ϕ)=:J(ϕ,z,t),ψvψ=(1ϕ)φq+γ(z)((1ϕ)φσV(φ)φϕW(ϕ))=:F(ϕ,φ,z,t) (2.8)

    for the total fluxes of Eq (2.1) and Eq (2.2). For ϕ<1, we define the final flux function

    ˜F(ϕ,ψ,z,t):=F(ϕ,ψ1ϕ,z,t)=ψq+γ(z)(ψσV(ψ1ϕ)ψϕW(ϕ)1ϕ), (2.9)

    whereas for ϕ=1, we set ˜F(1,ψ,z,t):=0 (since F(1,φ,z,t)=0 for all φ[0,1]). Substituting Eq (2.8) and Eq (2.9) into Eq (2.1) and Eq (2.2), respectively, we obtain the final governing PDE system Eq (1.1).

    Illustrations and numerical examples are based on the expressions

    W(ϕ)=vterm,p(1ϕ)npfor 0ϕ1,np>1, (2.10)
    V(φ)=vterm,s(1φ)nsfor 0φ1,ns>1 (2.11)

    (see [46]), where vterm,p and vterm,s are the terminal velocities of a single particle of the primary and secondary disperse phases, respectively, in an unbounded fluid. We set np=3.2, vterm,p=2.7cm/s, ns=2.5, and vterm,s=0.5cm/s along with σ=1. These values are used in Applications 1 and 2 in Section 6. The resulting nonlinearities of J(ϕ,z,t) and ˜F(ϕ,ψ,z,t) in the different zones of the column are illustrated in [52,Figure 3.2].

    The discretization of the model is based on the triangularity of the system of conservation laws Eq (1.1). The numerical fluxes for ϕ are based on the particular treatment of conservation laws with fluxes having an explicit "concentration times velocity" structure [6]. In each time step, an approximate solution ϕ of the first PDE of Eq (1.1) is obtained and used as a given piecewise constant function in space and time in the second PDE of Eq (1.1), which is updated accordingly.

    We define a computational domain [0,zend) (to be used for the error calculation; see Section 6.1) consisting of N cells by covering the vessel with N2 cells and placing one cell each below and above; see Figure 2 (left). This setup, with a finite spatial domain, is introduced for practical reasons and is the minimal spatial domain that captures the interior of the tank and the concentrations in the underflow and effluent zones. The formulation of the scheme (i.e., Scheme 1) and subsequent proof of invariant region property are referred to this computational domain, but for the convergence analysis the model is specified as the initial value problem Eq (1.1), Eq (1.2) with the initial data posed on the real line. This distinction is merely a formal one since on (,0) and (zend,) the model reduces to linear advection equations describing that matter is transported away from the unit at constant velocity (if no changes in A in these zones arise).

    Figure 2.  (Left) Discretization of ϕ and ψ in the application to flotation, where the height of the vessel is H=zEzU, there are K inlets, and the cross-sectional area A(z) has two values separated by a discontinuity at z=zF,2; cf. the examples in Sections 6.4 and 6.5, (right) enlarged view illustrating cell division for error computations when Δrz is the discretization of the reference solution (see Section 6.1).

    Given the column height H, we define Δz:=H/(N2), the cell boundaries zi:=iΔz, i=0,1,,N, and the cells (intervals) Ii1/2:=[zi1,zi) and Ii:=[zi1/2,zi+1/2). We place the column between zU:=Δz=z1 and zE:=zU+H=(N1)Δz=zN1. Then the length of the interval of error calculation is zend:=H+2Δz=NΔz. Each injection point zF,k is assumed to belong to one cell Ii1/2 and we define the dimensionless quantity

    δk,i1/2:=Ii1/2δzF,k(z)dz:={1if zF,kIi1/2,0otherwise. (3.1)

    The cross-sectional area A=A(z) is allowed to have a finite number of discontinuities and it is discretized by

    Ai:=1ΔzIiA(z)dz,Ai+1/2:=1ΔzIi+1/2A(z)dz.

    We simulate NT time steps up to the final time T:=NTΔt, with the fixed time step Δt satisfying the Courant-Friedrichs-Lewy (CFL) condition

    Δt(2Q,TAmin+M(max{V(0),V}+W+W))Δz, (3.2)

    where

    M:=maxi=1,2,,N{Ai1Ai1/2,AiAi1/2},Amin:=mink=0,1/2,1,3/2,,NAk,Q,T:=max0tTKk=1QF,k(t),W:=max0ϕ1|W(ϕ)|.

    Finally, we set tn:=nΔt for n=0,1,,NT.

    The time-dependent feed functions are discretized as

    QnF,k:=1Δttn+1tnQF,k(t)dt,ϕnF,k:=1Δttn+1tnϕF,k(t)dt,

    for k=1,,K, and the same is made for ψF,k.

    The first equation of Eq (1.1) is discretized by combining upwind discretizations of qϕ with the particular scheme proposed in [6] for models with a "concentration times velocity" flux, as is the case for the term ϕW(ϕ).

    The initial data are discretized by

    ϕ0i1/2:=1Ai1/2ΔzIi1/2ϕ(z,0)A(z)dz.

    To advance from tn to tn+1 from given values ϕni1/2, i=1,,N, we define the numerical flux at z=zi by

    Jni:={ϕn1/2qn0for i=0,ϕni1/2qn+i+ϕni+1/2qni+γiϕni1/2W(ϕni+1/2)for i=1,,N1,ϕnN1/2qn+Nfor i=N, (3.3)

    where the notation

    a+:=max{a,0},a:=min{a,0},γi:=γ(zi),andqn+i:=(q(zi,tn))+

    is used. Since the bulk fluxes above and below the tank are directed away from it,

    ϕn1/2qn+0=0andϕnN+1/2qnN=0for any values of ϕn1/2andϕnN+1/2.

    To simplify the presentation, we use the middle line of Eq (3.3) as the definition of Jni together with ϕn1/2:=0 and ϕnN+1/2:=0. With the notation λ:=Δt/Δz and Qn+i:=Aiqn+i etc., the conservation law on Ii1/2 implies the update formula

    ϕn+1i1/2=ϕni1/2+λAi1/2(Ai1Jni1AiJni+Kk=1QnF,kϕnF,kδk,i1/2)=:Hi1/2(ϕni3/2,ϕni1/2,ϕni+1/2),i=1,,N. (3.4)

    Then we define the piecewise constant approximate solution ϕΔz on R×[0,T) by

    ϕΔz(z,t):=i,nχIi1/2(z)χ[tn,tn+1)(t)ϕni1/2, (3.5)

    where χΩ denotes the characteristic function of the set Ω.

    We discretize the initial data by

    ψ0i1/2:=1Ai1/2ΔzIi1/2ψ(z,0)A(z)dz.

    The well-known Engquist-Osher numerical flux [26] for a given continuous, piecewise differentiable flux function g and real values a and b on the left/right is given by

    G(g;a,b):=g(0)+a0max{0,g(s)}ds+b0min{0,g(s)}ds. (3.6)

    Then a consistent numerical flux corresponding to Eq (2.9) is

    Fni:=ψni1/2qn+i+ψni+1/2qni+γi(Gni(ψni1/2,ψni+1/2)ϕni1/2ψni+1/21ϕni+1/2W(ϕni+1/2)),i=0,,N,

    where we set ψn1/2:=0 and ψnN+1/2:=0 with the same motivation as for ϕ above (these values are irrelevant). Here

    Gni(ψni1/2,ψni+1/2):=G(σfni;ψni1/2,ψni+1/2) (3.7)

    is the Engquist-Osher numerical flux associated with the function

    σfni(ψ):=σψ˜V(ψψnmax,i),˜V(u):={V(u)for u<1,0for u1, (3.8)

    where (ab:=min{a,b}, ab:=max{a,b})

    ψnmax,i:=(1ϕni1/2)(1ϕni+1/2)=1(ϕni1/2ϕni+1/2).

    The function ψσfni(ψ) is unimodal. Let ˆψni denote the maximum point of fni. For a given function ˜V the values ˆψni and ψnmax,i are related by the following lemma.

    Lemma 3.1. Assume that 0<ω<˜ω<1 are the unique local maximum and inflection point, respectively, of f(φ)=φV(φ) (cf. Eq (2.7)). Then ˆψni=ωψnmax,i for all i and n and all possible values 0ψnmax,i1. Moreover, the unique inflection point ψninfl,i(ˆψni,ψnmax,i) of fni satisfies ψninfl,i=˜ωψnmax,i for all i and n and all possible values 0ψnmax,i1. (See Figure 3.)

    Figure 3.  Illustration of Lemma 3.1.

    Proof. Assume that 0<ψnmax,i1. Since ˆψni is the unique solution ˆψni<ψnmax,i of

    ddψ(ψ˜V(ψψnmax,i))=0˜V(ψψnmax,i)+ψψnmax,i˜V(ψψnmax,i)=0,

    it follows that ω is the unique solution in (0,1) of ˜V(ω)+ω˜V(ω)=0 (cf. Eq (2.7)). By a similar argument, ˜ω is the unique solution of 2˜V(˜ω)+˜ω˜V(˜ω)=0.

    The Engquist-Osher numerical flux Eq (3.7), which in this form appears in Scheme 1 as well as in its reduced versions, Schemes 2 and 3, can now be computed as follows, where we recall that fni(0)=0. For σ=1 we get

    ψni1/20max{0,(fni)(s)}ds={fni(ψni1/2)if ψni1/2ˆψni,fni(ˆψni)if ψni1/2>ˆψni,ψni+1/20min{0,(fni)(s)}ds={0if ψni+1/2ˆψni,fni(ψni+1/2)fni(ˆψni)if ψni+1/2>ˆψni, (3.9)

    hence

    G(fni;ψni1/2,ψni+1/2)={fni(ψni1/2)if ψni1/2ˆψniandψni+1/2ˆψni,fni(ψni1/2)+fni(ψni+1/2)fni(ˆψni)if ψni1/2ˆψniandψni+1/2>ˆψni,fni(ˆψni)if ψni1/2>ˆψniandψni+1/2ˆψni,fni(ψni+1/2)if ψni1/2>ˆψniandψni+1/2>ˆψni. (3.10)

    By analogous reasoning we obtain for σ=1

    G(fni;ψni1/2,ψni+1/2)={fni(ψni+1/2)if ψni1/2ˆψniandψni+1/2ˆψni,fni(ˆψni)if ψni1/2ˆψniandψni+1/2>ˆψni,fni(ˆψni)fni(ψni1/2)fni(ψni+1/2)if ψni1/2>ˆψniandψni+1/2ˆψni,fni(ψni1/2)if ψni1/2>ˆψniandψni+1/2>ˆψni. (3.11)

    We define the difference operators Δai:=aiai1 and Δ+ai:=ai+1ai. Then the marching formula for ψ-Scheme 1 is

    ψn+1i1/2=ψni1/2+λAi1/2(Ai1Fni1AiFni+Kk=1QnF,kψnF,kδk,i1/2)=ψni1/2λAi1/2(Δ(ψni1/2Qn+i+ψni+1/2Qni+(Aγ)i(Gni(ψni1/2,ψni+1/2)ϕni1/2ψni+1/21ϕni+1/2W(ϕni+1/2)))Kk=1QnF,kψnF,kδk,i1/2),i=1,,N. (3.12)

    Then we define the piecewise constant approximate solution ψΔz on R×[0,T) by

    ψΔz(z,t):=i,nχIi1/2(z)χ[tn,tn+1)(t)ψni1/2. (3.13)

    We prove that Scheme 1, defined by the update formulas Eq (3.4) and Eq (3.12), is monotone, a property which then is used to prove the invariant-region property that the approximate solutions are positive and bounded.

    Theorem 3.1. If the CFL condition Eq (3.2) is satisfied, then the update formula for ϕ, Eq (3.4) (that is, ϕ-Scheme 1) is monotone and 0ϕni1/21 for i=1,,N and n=1,,NT.

    Proof. We recall the assumption Eq (1.3). We first prove monotonicity of the three-point scheme for ϕ Eq (3.4), i.e, that ϕn+1i1/2/ϕnk1/20 for all i=1,,N and k=i1,i,i+1. We have

    ϕn+1i1/2ϕni3/2=λAi1/2(Qn+i1+(Aγ)i1W(ϕni1/2))0,ϕn+1i1/2ϕni+1/2=λAi1/2(Qni(Aγ)iϕni1/2W(ϕni+1/2))0,ϕn+1i1/2ϕni1/2=1+λAi1/2(Qni1+(Aγ)i1ϕni3/2W(ϕni1/2)Qn+i(Aγ)iW(ϕni+1/2))1λ(2Q,TAmin+M(W+W))0,

    where we have used the CFL condition Eq (3.2).

    We now prove that if 0ϕni1/21 for all i, then 0ϕn+1i1/21 for all i. Clearly, Eq (1.3) implies that 0ϕ0i1/21 for all i. Since the scheme Eq (3.4) is monotone, Hi1/2 is non-decreasing in each argument. Since by assumption W(1)=0, we get the following estimation (where we use a++a=a):

    0λAi1/2Kk=1QnF,kϕnF,kδk,i1/2=Hi1/2(0,0,0)ϕn+1i1/2=Hi1/2(ϕni3/2,ϕni1/2,ϕni+1/2)Hi1/2(1,1,1)=1+λAi1/2((Qni1Qni)+Kk=1QnF,kϕnF,kδk,i1/2)1+λAi1/2(Kk=1(QnF,k)δk,i1/2+Kk=1QnF,kδk,i1/2)=1.

    Lemma 3.2. The function fni (cf. Eq (3.8)) satisfies (fni)max{V(0),V}.

    Proof. By Eq (2.7), the function f(φ)=φV(φ) has a single inflection point ˜ω(0,1) and by Lemma 3.1, fni has the inflection point ˜ωψnmax,i(0,ψnmax,i). We have (fni)(0)=V(0), (fni)(φ)=0 for ψnmax,iφ1 and the lowest (and negative) value of (fni) is obtained at its only critical point ˜ωψnmax,i, for which

    (fni)(˜ωψnmax,i)=˜V(˜ω)+˜ω˜V(˜ω)V.

    This concludes the proof.

    Lemma 3.3. There holds Gni(1ϕni1/2,1ϕni+1/2)=0 for all i and n.

    Proof. Assume that 0<ψnmax,i=(1ϕni1/2)(1ϕni+1/2)1. By Lemma 3.1, ˆψni<ψnmax,i, hence Eq (3.10), Eq (3.11), and

    ˜V((1ϕni1/2)/ψnmax,i)=˜V((1ϕni+1/2)/ψnmax,i)=0

    imply that

    Gni(1ϕni1/2,1ϕni+1/2)={fni(1ϕni+1/2)=0if σ=1,fni(1ϕni1/2)=0if σ=1.

    Theorem 3.2. Under the assumptions of Theorem 3.1, the update formula for ψ, Eq (3.12) (i.e., ψ-Scheme 1) is monotone and along with Eq (3.4) produces approximate solutions that satisfy 0ψni1/21ϕni1/2 for all i and n.

    Proof. Assumptions Eq (1.3) and Eq (1.5) imply that 0ψ0i1/21ϕ0i1/2 for all i and

    ψnF,k1ϕnF,kfor alln. (3.14)

    To prove that the scheme Eq (3.12) is monotone, we write it as

    ψn+1i1/2=Kni1/2(ψni3/2,ψni1/2,ψni+1/2) (3.15)

    and show that this expression is non-decreasing in each of its arguments.

    Since 0ϕni1/21 for a given n and all i, and appealing to Eq (3.6), we have

    ψn+1i1/2ψni3/2=λAi1/2(Qn+i1+(Aγ)i1Gni1ψni3/2)0,ψn+1i1/2ψni+1/2=λAi1/2(Qni(Aγ)iGniψni+1/2+(Aγ)iϕni1/21ϕni+1/2W(ϕni+1/2))0,ψn+1i1/2ψni1/2=1+λAi1/2(Qni1Qn+i+(Aγ)i1(Gni1ψni1/2ϕni3/2W(ϕni1/2)1ϕni1/2)(Aγ)iGniψni1/2)1λ(2Q,TAmin+M(Gniψni1/2Gni1ψni1/2+W(ϕni1/2)1ϕni1/2)).

    By Eq (3.6) and Lemma 3.2 we also obtain

    Gniψni1/2Gni1ψni1/2=(fni)(ψni1/2)+(fni)(ψni1/2)=|(fni)(ψni1/2)|(fni)max{V(0),V},

    and for the remaining term, we use that W(1)=0 to get

    W(ϕni1/2)1ϕni1/2=W(ϕni1/2)W(1)1ϕni1/2=W(ξ)Wfor some ξ(ϕni1/2,1).

    Hence, the CFL condition Eq (3.2) implies

    ψn+1i1/2ψni1/21λ(2Q,TAmin+M(max{V(0),V}+W))0.

    The inequalities proved imply that Kni1/2 is non-decreasing in each of its arguments. Now we use that 0ψni1/21ϕni1/2 for all i and Lemma 3.3 to obtain

    0λAi1/2Kk=1QnF,kψnF,kδk,i1/2=Hi1/2(0,0,0)ψn+1i1/2=Hi1/2(ψni3/2,ψni1/2,ψni+1/2)Hi1/2(1ϕni3/2,1ϕni1/2,1ϕni+1/2)=1ϕni1/2+λAi1/2(Ai1Fni1(1ϕni3/2,1ϕni1/2)AiFni(1ϕni1/2,1ϕni+1/2)+Kk=1QnF,kψnF,kδk,i1/2)=1ϕni1/2+λAi1/2((1ϕni3/2)Qn+i1+(1ϕni1/2)Qni1(Aγ)i1ϕni3/2W(ϕni1/2)(1ϕni1/2)Qn+i(1ϕni+1/2)Qni+(Aγ)iϕni1/2W(ϕni+1/2)+Kk=1QnF,kψnF,kδk,i1/2).

    Appealing to Eq (3.14) and the update formula for ϕ Eq (3.4), we get

    ψn+1i1/21ϕn+1i1/2+λAi1/2(Qn+i1+Qni1Qn+iQni+Kk=1QnF,kδk,i1/2)=1ϕn+1i1/2+λAi1/2{Qni1Qni+Kk=1QnF,kδk,i1/2}=1ϕn+1i1/2.

    The last equality holds since {}=0 irrespective of whether there is a source in the cell; Qni1Qni+QnF,k=0, or not; Qni1Qni=0.

    For ease of the argument, let us focus on the case of a constant interior cross-sectional area A, i.e., assume that Eq (1.7) is in effect. In addition, we assume that QnF,k, ϕnF,k, and ψnF,k (k=1,,K) are constant and therefore do not depend on n. The same is assumed for the underflow volumetric flow QU. (That is, we now study Scheme 2 suitable for Model 2.) Then Eq (3.4) and Eq (3.12) take the forms

    ϕn+1i1/2=ϕni1/2λΔ(ϕni1/2q+i+ϕni+1/2qi+γiϕni1/2W(ϕni+1/2))+λKk=1qF,kϕF,kδk,i1/2, (4.1)
    ψn+1i1/2=ψni1/2λΔ(ψni1/2q+i+ψni+1/2qi+γi(Gni(ψni1/2,ψni+1/2)ϕni1/2ψni+1/2W(ϕni+1/2)1ϕni+1/2))+λKk=1qF,kψF,kδk,i1/2, (4.2)

    where qF,k:=QF,k/A. To embed the treatment into available analyses of schemes for conservation laws with discontinuous flux, we absorb the feed terms into the numerical flux. That is, we define ik:=i if δk,i1/2=1 (see Eq (3.1)). Then

    qi={qUif ii11,qU+qF,1++qF,lif iliil+11,l=1,,K1,qU+qF,1++qF,Kfor iiK. (4.3)

    Furthermore, we define the feed flux

    hF,i:={0if ii11,qF,1ϕF,1++qF,lϕF,lif iliil+11,l=1,,K1,qF,1ϕF,1++qF,KϕF,Kfor iiK, (4.4)

    such that

    hF,ihF,i1=Kk=1qF,kϕF,kδk,i1/2.

    Consequently, we may write the scheme (i.e., ϕ-Scheme 2) as

    ϕn+1i1/2=ϕni1/2λΔ(ϕni+1/2qi+ϕni1/2q+i+γiϕni1/2W(ϕni+1/2)+hF,i). (4.5)

    For later use we define the piecewise constant functions

    q(z):=qkandhF(z):=hF,kfor zF,k<z<zF,k+1,k=0,,K,

    where zF,0:=, zF,K+1:=, and we define the function

    h(z,v,u):=q(z)v+q+(z)u+γ(z)uW(v)+hF(z) (4.6)

    that allows us to write Eq (4.5) as

    ϕn+1i1/2=ϕni1/2λΔh(zi,ϕni+1/2,ϕni1/2). (4.7)

    The PDE for ϕ within Model 2, that is when the simplification Eq (1.7) is applied to Model 1, is the conservation law

    tϕ+zJ(ϕ,z)=0,(z,t)ΠT (4.8)

    with discontinuous flux

    J(ϕ,z)={q(z)ϕKk=1qF,kϕF,kfor z>zE,q(z)ϕKk=1qF,kϕF,k+j(ϕ)for zF,K<z<zE,q(z)ϕlk=1qF,kϕF,k+j(ϕ)for zF,l<z<zF,l+1,l=1,,K1,qUϕ+j(ϕ)for zU<z<zF,1,qUϕfor z<zU. (4.9)

    posed along with the initial condition Eq (1.2a).

    The choice of the appropriate solution concept for weak solutions, and the ways we may relate the model to the available theory of conservation laws with discontinuous flux, requires verifying whether J(ϕ,z) as given by Eq (4.9) satisfies the so-called "crossing condition" across each discontinuity

    zZ:={zU,zF,1,,zF,K,zE}. (4.10)

    Certain early well-posedness (existence, stability, and uniqueness) results for conservation laws with discontinuous flux (and related equations) rely on satisfaction of this condition [35], although later developments advance solution concepts that do not rely on satisfaction of the crossing condition [4,36,39]. In the present context this condition is satisfied for a particular discontinuity at z if the adjacent fluxes to the right and the left, J(ϕ,z+) and J(ϕ,z), satisfy

    ϕ1,ϕ2[0,1]:J(ϕ1,z+)J(ϕ1,z)<0<J(ϕ2,z+)J(ϕ2,z)ϕ1<ϕ2, (4.11)

    which means either the graphs of J(,z) and J(,z+) do not intersect, or if they do, there is at most one flux crossing ϕχ and the graph of J(,z) lies above that of J(,z+) to the left of ϕχ. For J(ϕ,z) as given by Eq (4.9) this condition is clearly satisfied for z{zE,zU} (considering that j(ϕ)>0 for 0<ϕ<1 implies that J(,z) and J(,z+) do not intersect in this case), while

    J(ϕ,z+F,l)J(ϕ,zF,l)=qF,l(ϕϕF,l)for l=1,,K.

    Thus, the crossing condition is satisfied also for z=zF,l, l=1,,K, since either ϕF,l=0 and the adjacent fluxes do not intersect in (0,1), or the intersection takes place at ϕχ=ϕF,l and Eq (4.11) holds since qF,l>0 for all l. The preceding consideration is analogous to the one for the simpler clarifier-thickener model (equivalent to K=1 in the present notation) studied e.g. in [13,14]. With the present analysis it is clear that the crossing condition is satisfied at each flux discontinuity zZ.

    Some of the available analyses refer to initial-value problems of the type

    tu+xF(u,x)=0for (x,t)ΠT,u(x,0)=u0(x)for xR,whereF(u,x):=H(x)g(u)+H(x)f(u) (4.12)

    where f and g are Lipschitz continuous functions of u denoting the "right" and "left" flux adjacent to a flux discontinuity across x=0 and H denotes the Heavyside function. The model problem Eq (4.12) features, of course, only one flux discontinuity (sitting at x=0), while Eq (4.9), Eq (1.2a) includes several of them at separate spatial locations. The study of Eq (4.12) is, however, sufficient for the analysis of each single flux discontinuity.

    Here we start from the concept of entropy solutions of type V introduced by Karlsen and Towers [36]. This concept does not appeal to the existence of traces of the unknown with respect to the interfaces zZ across which J(ϕ,z) is discontinuous. To state its adaptation to the situation at hand, we define the sets

    Π(K+3/2)T:=(zE,)×(0,T),Π(K+1/2)T:=(zF,K,zE)×(0,T),Π(k1/2)T:=(zF,k1,zF,k)×(0,T),k=2,,K,Π(1/2)T:=(zU,zF,1)×(0,T),Π(1/2)T:=(,zU)×(0,T).

    Definition 4.1. A measurable function ϕ=ϕ(z,t)L1(ΠT) is an entropy solution of type V of the initial-value problem Eq (4.8), Eq (1.2a) if it satisfies the following conditions:

    (ⅰ) The function ϕ belongs to L(ΠT); for a.e. (z,t)ΠT there holds ϕ(z,t)[0,1].

    (ⅱ) The function ϕ is a weak solution of Eq (4.8), i.e., for all smooth test functions ζ with compact support in ΠT,

    ΠT(ϕtζ+J(ϕ,z)zζ)dzdt=0. (4.13)

    (ⅲ) For all l=0,,K+2, for any nonnegative smooth test function ζ(l) with compact support in Π(l)T and all c[0,1] there holds

    ΠT(|ϕc|tζ(l)+sgn(ϕc)(J(ϕ,z)J(c,z))zζ(l))dzdt+R|ϕ0c|ζ(l)(z,0)dt0. (4.14)

    (ⅳ) The following Kružkov-type [37] entropy inequality holds for all nonnnegative smooth test functions ζ with compact support in ΠT and all constants cR:

    ΠT(|ϕc|tζ+sgn(ϕc)(J(ϕ,z)J(c,z))zζ)dzdt+T0zZ|J(c,z+)J(c,z)|ζ(z,t)dt0. (4.15)

    Note that the entropy inequality Eq (4.15) does not imply the weak formulation Eq (4.13). The standard derivation of the weak formulation from the Kružkov entropy inequality (e.g., [30,Section 2.1]) does not apply here since some of the flux differences |J(c,z+)J(c,z)| are not compactly supported with respect to c, cf. [13,Rem. 1.1].

    Lemma 4.1. Consider ϕ-Scheme 2 applied to Eq (4.8), Eq (1.2a). There exists a constant C1, depending on TV(ϕ0), such that

    ΔziZ|ϕn+1i1/2ϕni1/2|ΔziZ|ϕ1i1/2ϕ0i1/2|C1Δt.

    Proof. Subtracting from Eq (4.1) its version from the previous time step, we get

    ϕn+1i1/2ϕni1/2=(ϕni3/2ϕn1i3/2)λBni1/2+(ϕni1/2ϕn1i1/2){1λBni+1/2+λCni1/2}+(ϕni+1/2ϕn1i+1/2){λCni+1/2},

    where we define

    Bni1/2:=q+i1+γi1W(ϕni1/2),Cni+1/2:={qi+γiϕn1i1/2W(ϕni+1/2)W(ϕn1i+1/2)ϕni+1/2ϕn1i+1/2if ϕni+1/2ϕn1i+1/2,0otherwise.

    Clearly Bni1/20, Cni+1/20, and due to the CFL condition,

    1λBni+1/2+λCni1/20,

    hence taking absolute values and summing over iZ we get, by appealing to standard arguments, that

    ΔziZ|ϕn+1i1/2ϕni1/2|ΔziZ|ϕni1/2ϕn1i1/2|ΔziZ|ϕ1i1/2ϕ0i1/2|.

    Furthermore, following the lines e.g. of the proof of [13,Lemma 3.2], we get that there exists a constant C2 that is independent of (Δt,Δz) such that

    iZ|ϕ1i1/2ϕ0i1/2|C2(TV(ϕ0)+TV(q)+TV(γ)),

    which completes the proof.

    A straightforward calculation yields that we can write the scheme in the form

    ϕn+1i1/2=ϕni1/2+CniΔ+ϕni1/2Dni1Δϕni1/2θni,

    where we define

    Cni:={λqiλγi1ϕni3/2ΔW(ϕni+1/2)Δ+ϕni1/2if Δ+ϕni1/20,λqiotherwise, Dni1:=λq+i+λγi1W(ϕni+1/2),θni:=λ(ϕni1/2Δqi+ϕni3/2Δq+i+ϕni1/2W(ϕni+1/2)ΔγiΔhF,i).

    The incremental coefficients satisfy Cni0 and Dni0; furthermore, the CFL condition ensures that Cni+Dni1 (in all cases for all i and n). Notice that θni=0 with the possible exception for those indices i at which Δqi0, Δq+i0, or Δγi0. According to the definition of γi and that of qi, see Eq (4.3), this may occur at most at a finite number of indices. Precisely, we may assert that (see Eq (4.10))

    θni=0if zi1,ziZ,

    hence for all indices i with the exception of finitely many indices i such that |zjζ|Δz for some ζZ, the scheme is given by the incremental form

    ϕn+1i1/2=ϕni1/2+CniΔ+ϕni1/2Dni1Δϕni1/2

    with incremental coefficients Cni0, Dni0, and Cni+Dni1. This property, in conjunction with Lemma 4.1, shows that we may apply [15,Lemma 5.3] (which is essentially Lemma 4.2 of [6], where a proof can be found) to the situation at hand. From [15,Lemma 5.3] we deduce the following lemma, where Vba(g) denotes the total variation of a function zg(z) over the interval (a,b).

    Lemma 4.2. Consider ϕ-Scheme 2 applied to Eq (4.8), Eq (1.2a). For any interval [a,b] such that [a,b]Z= and any t[0,T] there exists a total variation bound

    Vba(ϕΔz(,t))C(a,b),

    where C(a,b) is independent of (Δx,Δt) and t for t[0,T].

    Finally, we have shown in Theorem 3.1 that ϕ-Scheme 1 Eq (3.4) is monotone. This applies, in particular, to the reduced ϕ-Scheme 2 Eq (4.1) or equivalently, Eq (4.5) or Eq (4.7). Thus, ϕ-Scheme 2 satisfies a discrete entropy inequality. The proof of the following lemma is identical to that of [36,Lemma 5.2], and is therefore omitted.

    Lemma 4.3. The scheme Eq (4.7) (ϕ-Scheme 2) satisfies the following entropy inequality for any ci3/2,ci1/2,ci+1/2[0,1]:

    |ϕn+1i1/2ci1/2||ϕni1/2ci1/2|λΔHniλsgn(ϕn+1i1/2ci1/2)Δh(zi,ϕni+1/2,ϕni1/2),

    where h is defined in Eq (4.6) and the numerical entropy flux Hni is defined by

    Hni:=h(zi,ϕni+1/2ci+1/2,ϕni1/2ci1/2)h(zi,ϕni+1/2ci+1/2,ϕni1/2ci1/2).

    We now may appeal to the results of [36] and argue as follows. Theorem 3.1 and Lemmas 4.1–4.3 ensure convergence of the functions ϕΔz to a weak solution of Eq (4.8), Eq (1.2a) that satisfies items (ⅰ), (ⅱ) and (ⅲ) of Definition 4.1. It also satisfies the entropy inquality Eq (4.15) arising in part (ⅳ) of that definition by utilizing the discrete entropy inequality stated in Lemma 4.3. Thus, we have proved the following theorem.

    Theorem 4.1. Suppose that assumptions Eq (1.3) to Eq (1.5) are in effect and that ϕΔz is defined by Eq (3.5), where the values ϕni1/2 are defined by the scheme Eq (4.5) (that is, ϕ-Scheme 2). Let Δt,Δz0 with λ=Δt/Δz=const. such that the CFL condition Eq (3.2) is satisfied. Then ϕΔz converges in L1loc(ΠT) and a.e. in ΠT to an entropy solution of type V of the initial-value problem Eq (4.8), Eq (1.2a).

    Next, we deal with the marching formula Eq (4.2). To this end, we define a feed flux ˜hF,i exactly as in Eq (4.4) but with ϕF,i replaced by ψF,i for i=1,,K. Furthermore, we recall that ˜vϕs(ϕ)=W(ϕ)/(1ϕ). Thus, the scheme can be written as

    ψn+1i1/2=ψni1/2λΔ(˜hF,i+ψni1/2q+i+ψni+1/2qi+γi(Gni(ψni1/2,ψni+1/2)ϕni1/2˜vϕs(ϕni+1/2)ψni+1/2)). (4.16)

    Lemma 4.4 (Crandall and Tartar [19]). Assume that (Ω,μ) is some measure space and that D is a subset of L1(Ω) with the property that if u,vD, then (uv)=max{u,v}D. Assume that T is a map T:DuT(u)D such that

    ΩT(u)dμ=Ωudμfor all uD.

    Then the following statements, valid for all u,vD, are equivalent:

    (ⅰ) If uv, then T(u)T(v).

    (ⅱ) Ω((T(u)T(v))0)dμΩ((uv)0)dμ.

    (ⅲ) Ω|T(u)T(v)|dμΩ|uv|dμ.

    Following, for instance, [18], we utilize Lemma 4.4 for the following mapping. Assume that DL1(R) is the set of piecewise constant functions and that are constant on the intervals Ii1/2 for iZ, and that with the marching formula Eq (3.15) we associate an operator Kn:DD such that if ψΔz(,tn) is the piecewise constant function defined by Eq (3.13) for t=tn, we may write ψ-Scheme 2 as

    ψΔz(,tn+1)=Kn(ψΔz(,tn)).

    Clearly, the monotonicity of the scheme implies that if u, v \in D , then

    \begin{align*} u \leq v \Rightarrow \mathcal{K}^n (u ) \leq \mathcal{K}^n ( v). \end{align*}

    Thus, Lemma 4.4 (ⅰ) holds. For u = \psi^{\Delta z} (\cdot, t_n) and v = \psi^{\Delta z} (\cdot, t_{n-1}) , Lemma 4.4 (ⅲ) implies that

    \begin{align*} \Delta z \sum\limits_{i \in \mathbb{Z}} \bigl| \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^{n} \bigr| & = \int_{\mathbb{R}} \bigl| \psi^{\Delta z} ( z , t_{n+1}) - \psi^{\Delta z} ( z , t_n) \bigr| \, \mathrm{d} z \\ & \leq \int_{\mathbb{R}} \bigl| \psi^{\Delta z} ( z , t_{n}) - \psi^{\Delta z} ( z , t_{n-1}) \bigr| \, \mathrm{d} z = \Delta z \sum\limits_{i \in \mathbb{Z}} \bigl| \psi_{i-1/2}^{n} - \psi_{i-1/2}^{n-1} \bigr| \end{align*}

    and therefore

    \begin{align*} & \Delta z \sum\limits_{i \in \mathbb{Z}} \bigl| \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^{n} \bigr| \leq \Delta z \sum\limits_{i \in \mathbb{Z}} \bigl| \psi_{i-1/2}^{1} - \psi_{i-1/2}^{0} \bigr| . \end{align*}

    However, we may assert that there exists a constant C_3 , which is independent of (\Delta t, \Delta x) , such that

    \begin{align*} & \sum\limits_{i \in \mathbb{Z}} \bigl| \psi_{i-1/2}^{1} - \psi_{i-1/2}^{0} \bigr| \\ & = \sum\limits_{i \in \mathbb{Z}} \biggl| \Delta_- \bigg( \psi_{i-1/2}^{0}q_{i}^{+}+ \psi_{i+1/2}^0 q_{i}^{-} + \gamma_{i}\biggl( G_{i}^{0}\big(\psi_{i-1/2}^{0}, \psi_{i+1/2}^{0}\big) -\phi_{i-1/2}^n\frac{\psi_{i+1/2}^0 W(\phi_{i+1/2}^{0}) }{1-\phi_{i+1/2}^{0}} \biggr) \bigg)\\ & \quad - \lambda \sum\limits_{k = 1}^K \frac{Q_{\mathrm{F}, k}}{A} \psi_{\mathrm{F}, k}^0\delta_{k, i-1/2} \biggr| \leq C_3. \end{align*}

    Since Eq (1.4) is a sufficient condition for this bound on the initial discrete divergence to hold, we get

    \begin{align*} & \Delta z \sum\limits_{i \in \mathbb{Z}} \bigl| \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^{n} \bigr| \leq \Delta z C_3 = \frac{\Delta t}{\lambda} C_3 . \end{align*}

    Consequently, we have proved the following lemma.

    Lemma 4.5. There exists a constant C_4 that is independent of (\Delta t, \Delta z) such that the numerical approximations produced by \psi -Scheme 2 satisfy

    \begin{align*} \Delta z \sum\limits_{i \in \mathbb{Z}} \bigl| \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^{n} \bigr| \leq C_4 \Delta t. \end{align*}

    To write down Scheme 1 in the simplest setting possible, we consider the model and numerical scheme under the assumptions before, and additionally assume a constant bulk velocity q , that the feed terms (giving rise to the singular source) are zero, and set the parameter \gamma = 1 . Thus, the model reduces to the triangular system of conservation laws Eq (1.8) with the initial conditions Eq (1.2) (Model 3), where we recall that assumptions Eq (1.3) are in effect, and Scheme 2, which in turn is a reduced form of Scheme 1, now is further reduced to Scheme 3.

    Assume now that \eta = \eta (\psi) is a smooth convex entropy function and Q = Q(\phi, \psi) is the corresponding compatible entropy flux compatible with Eq (1.8b), i.e.,

    \begin{align} \partial_{\psi} {Q} ( \phi, \psi) = \eta' ( \psi) \partial_{\psi} \tilde{F} ( \phi, \psi). \end{align} (5.1)

    In what follows, we refer to (\eta, Q) as an entropy pair for Eq (1.8b). In particular we denote by (\eta_0, Q_0) the Kružkov entropy pair [37], that is

    \begin{align} \eta_0 (\psi) = |\psi- c|, \quad Q_0 ( \phi, \psi) = \operatorname*{sgn} ( \psi-c) \bigl( \tilde{F} ( \phi, \psi) - \tilde{F} ( \phi, c) \bigr), \end{align} (5.2)

    where c \in \mathbb{R} is a constant.

    The convergence proof is based on the following lemma, slightly adapted from [18,Lemma 2.2], which in turn is an adaptation of [43,Theorem 5] (see also [44]).

    Lemma 5.1. Let \phi be the unique entropy solution of the initial-value problem Eq (1.8a), Eq (1.2a), and assume that \smash{\{ \psi^{\nu} \}_{\nu > 0}} is a family of functions defined on \Pi_T . If \smash{\{ \psi^{\nu} \}} is bounded in L^{\infty} (\Pi_T) and \{ \partial_t \eta_0 (\psi^{\nu}) + \partial_z Q_0 (\phi, \psi^{\nu})\}_{\nu > 0} lies in a compact set of H_{\mathrm{loc}}^{-1} (\Pi_T) for all constants c , then there exists a sequence \{ \nu_n \}_{n \in \mathbb{N}} such that \nu_n \to 0 as n \to \infty and a function \psi \in L^{\infty} (\Pi_T) such that

    \begin{align*} \psi^{\nu_n} \to \psi \quad \text{a.e. and in } L_{\mathrm{loc}}^p ( \Pi_T) , 1 \leq p < \infty . \end{align*}

    Consistently with Eq (4.6), Eq (4.7) we assume that the scheme employed to approximate entropy solutions of Eq (1.8a) is \phi -Scheme 3, that is,

    \begin{align*} \phi_{i-1/2}^{n+1 } = \phi_{i-1/2}^{n} - \lambda \Delta_- h \bigl( \phi_{i+1/2}^n, \phi_{i-1/2}^n \bigr), \quad h(v, u) : = q^- v +q^+ u + u W(v). \end{align*}

    Clearly, under a suitable CFL condition, the \phi -Scheme 3 converges to the unique entropy solution of Eq (1.8a), Eq (1.2a). Our goal is to establish convergence of the corresponding scheme for \psi ( \psi -Scheme 3). We here write the scheme as

    \begin{align} \psi_{i-1/2}^{n+1}& = \psi_{i-1/2}^n - \lambda \Delta_- \mathcal{F} \bigl( \phi_{i-1/2}^n, \phi_{i+1/2}^n, \psi_{i-1/2}^n, \psi_{i+1/2}^n \bigr) \equiv \psi_{i-1/2}^n - \lambda \Delta_- \mathcal{F} ( \boldsymbol{\phi}_i^n, \boldsymbol{\psi}_i^n), \end{align} (5.3)

    where we define the four-argument numerical flux

    \begin{align} \mathcal{F} ( a, b, u, v ) : = q^+ u + q^- v + \bigl( G(a, b, u, v) - a \tilde{v}_{\phi \mathrm{s}} (b) v \bigr), \end{align} (5.4)

    denote pairs of neighboring \phi - and \psi -values by

    \begin{align*} \boldsymbol{\phi}_i^n : = \bigl( \phi_{i-1/2}^n, \phi_{i+1/2}^n \bigr) \quad \text{and} \quad \boldsymbol{\psi}_i^n : = \bigl( \psi_{i-1/2}^n, \psi_{i+1/2}^n \bigr), \end{align*}

    and replace the arguments " \smash{\phi_{i-1/2}^n, \phi_{i+1/2}^n} " by \smash{\boldsymbol{\phi}_i^n} (analogously for \psi ). In Eq (5.4) a and b play the roles of \smash{\phi_{i-1/2}^n} and \smash{\phi_{i+1/2}^n} , and u and v those of \smash{\psi_{i-1/2}^n} and \smash{\psi_{i+1/2}^n} , respectively, and we define G(a, b, u, v) as follows (cf. Eq (3.7), Eq (3.8)). Let

    \begin{align*} f( a, b , \psi) : = \psi \tilde{V} \left( \frac{\psi}{1- (a \vee b )} \right), \end{align*}

    then G(a, b, \cdot, \cdot) is the Engquist-Osher numerical flux [26] associated with f(a, b, \cdot) .

    The compensated compactness approach strongly depends on entropy inequalities satisfied by the scheme Eq (5.3). To prepare for the derivation of suitable uniform estimates, we multiply the scheme Eq (5.3) by \smash{\eta' (\psi_{i-1/2}^{n+1})} , where \eta is a smooth convex entropy function, and utilize the Taylor expansion

    \begin{align*} \eta' \bigl( \psi_{i-1/2}^{n+1} \bigr) \bigl( \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^{n} \bigr) & = \eta \bigl( \psi_{i-1/2}^{n+1} \bigr) - \eta \bigl( \psi_{i-1/2}^{n} \bigr) + \frac{1}{2} \eta'' \bigl( \xi_{i-1/2}^{n+1/2} \bigr) \bigl( \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^{n} \bigr)^2, \end{align*}

    where \smash{\xi_{i-1/2}^{n+1/2}} is an intermediate value between \smash{\psi_{i-1/2}^{n}} and \smash{\psi_{i-1/2}^{n+1}} . This yields

    \begin{align} \begin{split} & \eta \bigl( \psi_{i-1/2}^{n+1} \bigr) - \eta \bigl( \psi_{i-1/2}^{n} \bigr) + \frac{1}{2} \eta'' \bigl( \xi_{i-1/2}^{n+1/2} \bigr) \bigl( \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^{n} \bigr)^2 \\ & = - \lambda \eta' \bigl( \psi_{i-1/2}^{n+1} \bigr) \Delta_- \mathcal{F} ( \boldsymbol{\phi}_i^n, \boldsymbol{\psi}_i^n) \\ & = - \lambda \eta' \bigl( \psi_{i-1/2}^{n} \bigr) \Delta_- \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i) - \lambda \bigl( \eta' \bigl( \psi_{i-1/2}^{n+1} \bigr) - \eta' \bigl( \psi_{i-1/2}^{n} \bigr) \bigr) \Delta_- \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i). \end{split} \end{align} (5.5)

    We now define the functions \hat{f} and \check{f} as the partial derivatives

    \begin{align*} \hat{f} (a, b, u) & : = \partial_{u} \mathcal{F} (a, b, u, v) = q^+ + \partial_u G(a, b, u, v) \geq 0, \\ \check{f} (a, b, v) & : = \partial_{v} \mathcal{F} (a, b, u, v) = q^- + \big(\partial_v G(a, b, u, v) - a \tilde{v}_{\phi \mathrm{s}} (b) \leq 0. \end{align*}

    The dependence of \partial_{u} \mathcal{F} (a, b, u, v) and \partial_{v} \mathcal{F} (a, b, u, v) on u and v only, respectively, is crucial for the subsequent analysis. We define the functions

    \begin{align*} \hat{\mathcal{F}} (a, b, u): = \int_0^u \hat{f} (a, b, s) \, \mathrm{d} s, \quad \check{\mathcal{F}} (a, b, v) : = \int_0^v \check{f} (a, b, s) \, \mathrm{d} s \end{align*}

    and note that

    \begin{align} \mathcal{F} (a, b, u, v) = \hat{\mathcal{F}} (a, b, u) + \check{\mathcal{F}} (a, b, v). \end{align} (5.6)

    Next, we define

    \begin{align} \begin{split} \hat{\mathcal{Q}} (a, b, \psi)& : = \int_0^{\psi} \eta'(u) \hat{f} (a, b, u) \, \mathrm{d} u, \quad \check{\mathcal{Q}} (a, b, \psi ) : = \int_0^{\psi} \eta'(v) \check{f} (a, b, v ) \, \mathrm{d} v , \\ \mathcal{Q} (a, b, \psi_1, \psi_2 )& : = \hat{\mathcal{Q}} (a, b, \psi_1 ) + \check{\mathcal{Q}} (a, b, \psi_2 ). \end{split} \end{align} (5.7)

    The function \mathcal{Q} is a consistent numerical entropy flux for the scheme Eq (5.3) for the entropy function \eta since

    \begin{align*} \mathcal{Q} (a, a, \psi, \psi) & = \int_0^{\psi} \eta' (u) \bigl( \hat{f} (a, a, u) + \check{f} (a, a, u ) \bigr) \, \mathrm{d} u \\ & = \int_0^{\psi} \eta' (u) \partial_u \mathcal{F} (a, a, u, u) \, \mathrm{d} u = \int_0^{\psi} \eta' (u) \tilde{F} (a, u) \, \mathrm{d} u = Q(a, \psi). \end{align*}

    Furthermore, integration by parts yields

    \begin{array}{l} \hat{\mathcal{Q}} (a, b, \psi) - \hat{\mathcal{Q}} (a, b, \tilde{\psi} ) \\ = \eta' ( \psi) \bigl( \hat{\mathcal{F}} (a, b, \psi) - \hat{\mathcal{F}} (a, b, \tilde{\psi} ) \bigr) - \int_{\tilde{\psi}}^{\psi} \eta'' (u) \bigl( \hat{\mathcal{F}} (a, b, u ) - \hat{\mathcal{F}} (a, b, \tilde{\psi}) \bigr) \, \mathrm{d} u, \end{array} (5.8)
    \begin{array}{l} \check{\mathcal{Q}} (a, b, \psi) - \check{\mathcal{Q}} (a, b, \tilde{\psi} ) \\ = \eta' ( \psi ) \bigl( \check{\mathcal{F}} (a, b, \psi) - \check{\mathcal{F}} (a, b, \tilde{\psi} ) \bigr) - \int_{\tilde{\psi}}^{\psi} \eta'' (u) \bigl( \check{\mathcal{F}} (a, b, u ) - \check{\mathcal{F}} (a, b, \tilde{\psi} ) \bigr)\, \mathrm{d} u \end{array} (5.9)
    \begin{align} & = \eta' ( \tilde{\psi} ) \bigl( \check{\mathcal{F}} (a, b, \psi) - \check{\mathcal{F}} (a, b, \tilde{\psi} ) \bigr) - \int_{\tilde{\psi}}^{\psi} \eta'' (u) \bigl( \check{\mathcal{F}} (a, b, u ) - \check{\mathcal{F}} (a, b, \psi ) \bigr)\, \mathrm{d} u. \end{align} (5.10)

    Now, denoting by \Delta_-^{\phi} and \Delta_-^{\psi} difference operators that act on both \phi - and \psi -arguments only, respectively (leaving the two others unchanged), we observe that

    \begin{align} \Delta_- \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i) = \Delta_-^{\psi} \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i) + \Delta_+ ^{\phi} \mathcal{F} ( \boldsymbol{\phi}^n_{i-1}, \boldsymbol{\psi}^n_{i-1}). \end{align} (5.11)

    In light of Eq (5.8) and Eq (5.10),

    \begin{array}{l} \eta' \bigl( \psi_{i-1/2}^n \bigr) \Delta_-^{\psi} \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i \bigr) \\ = \hat{\mathcal{Q}} \bigl( \boldsymbol{\phi}_i^n, \psi_{i-1/2}^n \bigr) - \hat{\mathcal{Q}} \bigl( \boldsymbol{\phi}_i^n, \psi_{i-3/2}^n \bigr) + \check{\mathcal{Q}} \bigl(\boldsymbol{\phi}_i^n , \psi_{i+1/2}^n \bigr) - \check{\mathcal{Q}} \bigl(\boldsymbol{\phi}_i^n, \psi_{i-1/2}^n \bigr) \\ - \Biggl( \eta' ( \psi_{i-1/2}^n ) \bigl( \hat{\mathcal{F}} (\boldsymbol{\phi}_i^n, \psi_{i-1/2}^n ) - \hat{\mathcal{F}} (\boldsymbol{\phi}_i^n, \psi_{i-3/2}^n ) \bigr) \\ \qquad \qquad - \int_{\psi_{i-3/2}^n}^{\psi_{i-1/2}^n} \eta'' (u) \bigl( \hat{\mathcal{F}} ( \boldsymbol{\phi}_i^n, u ) - \hat{\mathcal{F}} (\boldsymbol{\phi}_i^n, \psi_{i-3/2}^n) \bigr) \, \mathrm{d} u \Biggr) \\ - \Biggl( \eta' ( \psi_{i-1/2}^n ) \bigl( \check{\mathcal{F}} (\boldsymbol{\phi}_i^n, \psi_{i+1/2}^n ) - \check{\mathcal{F}} (\boldsymbol{\phi}_i^n, \psi_{i-1/2}^n ) \bigr) \\ \qquad + \int_{\psi_{i+1/2}^n}^{\psi_{i-1/2}^n} \eta'' (u) \bigl( \check{\mathcal{F}} ( \boldsymbol{\phi}_i^n, u ) - \check{\mathcal{F}} (\boldsymbol{\phi}_i^n, \psi_{i+1/2}^n) \bigr) \, \mathrm{d} u \Biggr) \\ + \eta' \bigl( \psi_{i-1/2}^n \bigr) \bigl( \hat{\mathcal{F}} \bigl(\boldsymbol{\phi}_i^n, \psi_{i-1/2}^n \bigr) - \hat{\mathcal{F}} \bigl( \boldsymbol{\phi}_i^n, \psi_{i-3/2}^n \bigr) + \check{\mathcal{F}} \bigl( \boldsymbol{\phi}_i^n, \psi_{i+1/2}^n \bigr) - \check{\mathcal{F}} \bigl( \boldsymbol{\phi}_i^n, \psi_{i-1/2}^n \bigr) \bigr) \\ = \mathcal{Q} \bigl(\boldsymbol{\phi}_i^n, \psi_{i-1/2}^n, \psi_{i+1/2}^n \bigr) - \mathcal{Q} \bigl(\boldsymbol{\phi}_i^n, \psi_{i-3/2}^n, \psi_{i-1/2}^n \bigr) + \Theta_{i-1/2}^n \\ = \Delta_-^{\psi} \mathcal{Q} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i) + \Theta_{i-1/2}^n , \end{array} (5.12)

    where the notation for evaluations and differences for \mathcal{Q} is the same as for \mathcal{F} and \smash{\Theta_{i-1/2}^n : = \hat{\Theta}_{i-1}^n + \check{\Theta}_{i}^n} , where

    \begin{align*} \hat{\Theta}_{i-1}^n & : = \int_{\psi_{i-3/2}^n}^{\psi_{i-1/2}^n} \eta'' (u) \bigl( \hat{\mathcal{F}} \bigl(\boldsymbol{\phi}_i^n , u \bigr) - \hat{\mathcal{F}} \bigl( \boldsymbol{\phi}_i^n, \psi_{i-3/2}^n \bigr)\bigr) \, \mathrm{d} u, \\ \check{\Theta}_{i}^n & : = - \int_{\psi_{i+1/2}^n}^{\psi_{i-1/2}^n} \eta'' (u) \bigl( \check{\mathcal{F}} \bigl( \boldsymbol{\phi}_i^n, u \bigr) - \check{\mathcal{F}} \bigl( \boldsymbol{\phi}_i^n, \psi_{i+1/2}^n \bigr) \bigr) \, \mathrm{d} u. \end{align*}

    Since \hat{\mathcal{F}} is increasing and \check{\mathcal{F}} is decreasing in the respective third argument, there holds \smash{ \hat{\Theta}_{i-1}^n, \check{\Theta}_{i}^n \geq 0} and therefore \smash{\Theta_{i-1/2}^n \geq 0} . Furthermore, we notice that

    \begin{align} \eta' \bigl( \psi_{i-1/2}^n \bigr) \Delta_+^{\phi} \mathcal{F} ( \boldsymbol{\phi}^n_{i-1}, \boldsymbol{\psi}^n_{i-1}) = \Delta_+^{\phi} \bigl( \eta' \bigl( \psi_{i-1/2}^n \bigr) \mathcal{F} ( \boldsymbol{\phi}^n_{i-1}, \boldsymbol{\psi}^n_{i-1}) \bigr). \end{align} (5.13)

    From Eq (5.11) we obtain by taking into account Eq (5.12) and Eq (5.13)

    \begin{array}{l} \eta' \bigl( \psi_{i-1/2}^n \bigr) \Delta_- \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i) = \Delta_-^{\psi} \mathcal{Q} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i) + \eta' \bigl( \psi_{i-1/2}^n \bigr) \Delta_+^{\phi} \mathcal{F} ( \boldsymbol{\phi}^n_{i-1}, \boldsymbol{\psi}^n_{i-1}) + \Theta_{i-1/2}^n \nonumber \\ = \Delta_- \mathcal{Q} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i) - \Delta_+^{\phi} \mathcal{Q} ( \boldsymbol{\phi}^n_{i-1}, \boldsymbol{\psi}^n_{i-1}) + \eta' \bigl( \psi_{i-1/2}^n \bigr) \Delta_+^{\phi} \mathcal{F} ( \boldsymbol{\phi}^n_{i-1}, \boldsymbol{\psi}^n_{i-1}) + \Theta_{i-1/2}^n \nonumber \\ = \Delta_- \mathcal{Q} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i) + \Delta_+^{\phi} \bigl( \eta' \bigl( \psi_{i-1/2}^n \bigr) \mathcal{F} ( \boldsymbol{\phi}^n_{i-1}, \boldsymbol{\psi}^n_{i-1}) - \mathcal{Q} ( \boldsymbol{\phi}^n_{i-1}, \boldsymbol{\psi}^n_{i-1}) \bigr) + \Theta_{i-1/2}^n. \end{array}

    Consequently, Eq (5.5) can be written as

    \begin{align} \begin{split} & \eta \bigl( \psi_{i-1/2}^{n+1} \bigr) - \eta \bigl( \psi_{i-1/2}^{n} \bigr) + \frac{1}{2} \eta'' \bigl( \xi_{i-1/2}^{n+1/2} \bigr) \bigl( \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^{n} \bigr)^2 + \lambda \Theta_{i-1/2}^n \\ & = - \lambda \Delta_- \mathcal{Q} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i) - \lambda \bigl( \eta' \bigl( \psi_{i-1/2}^{n+1} \bigr) - \eta' \bigl( \psi_{i-1/2}^{n} \bigr) \bigr) \Delta_- \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i) \\ & \quad - \lambda \Delta_+^{\phi} \bigl( \eta' \bigl( \psi_{i-1/2}^n \bigr) \mathcal{F} ( \boldsymbol{\phi}^n_{i-1}, \boldsymbol{\psi}^n_{i-1}) - \mathcal{Q} ( \boldsymbol{\phi}^n_{i-1}, \boldsymbol{\psi}^n_{i-1}) \bigr) . \end{split} \end{align} (5.14)

    Multiplying Eq (5.14) by \Delta z and summing over (n, i) \in \mathcal{I}_1 , where

    \begin{align*} \mathcal{I}_k : = \{ (n, i) \mid n = 0 , \dots, N_T -k, \, i \in \mathbb{Z} \}, \end{align*}

    we get

    \begin{align*} & \Delta z \sum\limits_{i \in \mathbb{Z}} \eta \bigl( \psi_{i-1/2}^{N} \bigr) - \Delta z \sum\limits_{i \in \mathbb{Z}} \eta \bigl( \psi_{i-1/2}^{0} \bigr) + \frac{\Delta z}{2} \sum\limits_{ \mathcal{I}_1 } \eta'' \bigl( \xi_{i-1/2}^{n+1/2} \bigr) \bigl( \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^{n} \bigr)^2 + \lambda \Delta z \sum\limits_{\mathcal{I}_1} \Theta_{i-1/2}^n \\ & = - \lambda \Delta z \sum\limits_{\mathcal{I}_1} \Delta_- \mathcal{Q}( \boldsymbol{\phi}_i^n, \boldsymbol{\psi}_i^n) - \lambda \Delta z \sum\limits_{\mathcal{I}_1} \bigl( \eta' \bigl( \psi_{i-1/2}^{n+1} \bigr) - \eta' \bigl( \psi_{i-1/2}^{n} \bigr) \bigr) \Delta_- \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i) \\ & \quad - \lambda \Delta z \sum\limits_{\mathcal{I}_1} \Delta_+^{\phi} \bigl( \eta' \bigl( \psi_{i-1/2}^n \bigr) \mathcal{F} ( \boldsymbol{\phi}_{i-1}^n, \boldsymbol{\psi}_{i-1}^n) - \mathcal{Q} ( \boldsymbol{\phi}_{i-1}^n, \boldsymbol{\psi}_{i-1}^n) \bigr), \end{align*}

    which implies the inequality

    \begin{align*} & \Delta z \sum\limits_{i \in \mathbb{Z}} \eta \bigl( \psi_{i-1/2}^{N} \bigr) + \frac{\Delta z}{2}\sum\limits_{\mathcal{I}_1} \eta'' \bigl( \xi_{i-1/2}^{n+1/2} \bigr) \bigl( \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^{n} \bigr)^2 + \lambda \Delta z \sum\limits_{\mathcal{I}_1} \Theta_{i-1/2}^n \\ & \leq \Delta z \sum\limits_{i \in \mathbb{Z}} \eta \bigl( \psi_{i-1/2}^{0} \bigr) +2 \| \eta' \|_{L^{\infty}} \Delta z \Delta t \sum\limits_{\mathcal{I}_1} \frac{1}{\Delta z} \bigl| \Delta_- \mathcal{F} ( \boldsymbol{\phi}_i^n, \boldsymbol{\psi}_i^n) \bigr| \\ & \quad + C \Delta z \Delta t \sum\limits_{\mathcal{I}_1} \frac{1}{\Delta z} \bigl|\Delta_+^{\phi} \bigl( \eta' \bigl( \psi_{i-1/2}^n \bigr) \mathcal{F} ( \boldsymbol{\phi}_{i-1}^n, \boldsymbol{\psi}_{i-1}^n) - \mathcal{Q} ( \boldsymbol{\phi}_{i-1}^n, \boldsymbol{\psi}_{i-1}^n) \bigr) \bigr| . \end{align*}

    The last term on the right-hand side is uniformly bounded since \phi^{\Delta z} has bounded variation. Now let us choose \eta(v) = v^2 and take into account [34] that there exists a constant C_{\mathcal{F}} such that

    \begin{align*} \hat{\Theta}_{i-1}^n \geq \frac{1}{C_{\mathcal{F}}} \Bigl( \Delta_-^{\psi} \hat{\mathcal{F}} \bigl(\boldsymbol{\phi}_i^n , \psi_{i-1/2} \bigr) \Bigr)^2, \quad \check{\Theta}_{i}^n \geq \frac{1}{C_{\mathcal{F}}} \Bigl( \Delta_+^{\psi} \check{\mathcal{F}} \bigl( \boldsymbol{\phi}_i^n , \psi_{i-1/2} \bigr) \Bigr)^2. \end{align*}

    Noticing that Lemma 4.5, applied to the present scheme ( \psi -Scheme 3), implies the bound on the discrete divergence of the numerical flux

    \begin{align} \Delta z \sum\limits_{i \in \mathbb{Z}} \frac{1}{\Delta z} \bigl| \Delta_- \mathcal{F}_i ( \boldsymbol{\phi}^n, \boldsymbol{\psi}^n) \bigr| \leq C, \end{align} (5.15)

    we obtain from Eq (5.14)

    \begin{align} \begin{split} & \Delta z \sum\limits_{i \in \mathbb{Z}} \eta \bigl( \psi_{i-1/2}^{N} \bigr) + \Delta z \sum\limits_{\mathcal{I}_1} \bigl( \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^{n} \bigr)^2 + \frac{\lambda}{C_{\mathcal{F}}} \Delta z \sum\limits_{\mathcal{I}_1} \Bigl( \bigl( \Delta_-^{\psi} \hat{\mathcal{F}} \bigl( \boldsymbol{\phi}_i^n, \psi_{i-1/2} \bigr) \bigr)^2 + \bigl( \Delta_+^{\psi} \check{\mathcal{F}} \bigl( \boldsymbol{\phi}_i^n, \psi_{i-1/2}^n \bigr) \bigr)^2 \Bigr) \\ & \leq \Delta z \sum\limits_{i \in \mathbb{Z}} \bigl( \psi_{i-1/2}^{0} \bigr)^2 + C_{T}. \end{split} \end{align} (5.16)

    Inequality Eq (5.16) implies the following estimate.

    Lemma 5.2. Consider numerical approximations produced by Scheme 3. There exists a constant C_7 that is independent of (\Delta z, \Delta t) such that

    \begin{align} \Delta t \Delta z \sum\limits_{n = 0}^{N_T-1} \sum\limits_{i \in \mathbb{Z}} \bigl( \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^{n} \bigr)^2 \leq C_7 \Delta z. \end{align} (5.17)

    Proof. The estimate for the "time variation" of \psi^{\Delta z} , Eq (5.17), follows immediately from Eq (5.16) if we consider that its right-hand side is uniformly bounded.

    Before proceeding, we prove the following lemma that is crucial for the subsequent analysis. For ease of notation we define the difference operators \smash{\Delta_{\pm}^{(3)}} and \smash{\Delta_{\pm}^{(4)}} that only act on the third or fourth argument of a function, respectively.

    Lemma 5.3. Consider numerical approximations produced by Scheme 3. There exist constants C_8 and C_9 that are independent of (\Delta z, \Delta t) such that for all i ,

    \begin{align} \bigl| \Delta_- \mathcal{Q} ( \boldsymbol{\phi}_i^n, \boldsymbol{\psi}_i^n ) \bigr| & \leq C_8 \bigl| (\Delta_-^{(3)} + \Delta_-^{(4)}) \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i ) \bigr| + C_9 \bigl( \bigl|\Delta_- \phi_{i-1/2}^n \bigr| + \bigl| \Delta_+ \phi_{i-1/2}^n \bigr| \bigr). \end{align} (5.18)

    Proof. We note that

    \begin{align} \Delta_- \mathcal{Q} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i ) = \Delta_-^{\psi} \mathcal{Q} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i ) + \Delta_+^{\phi} \mathcal{Q} ( \boldsymbol{\phi}^n_{i-1}, \boldsymbol{\psi}^n_{i-1} ). \end{align} (5.19)

    We first discuss

    \begin{align*} \Delta_-^{\psi} \mathcal{Q} ( \boldsymbol{\phi}_i^n, \boldsymbol{\psi}_i^n ) = \Delta_-^{(3)} \hat{\mathcal{Q}} \bigl( \boldsymbol{\phi}_i^n , \psi_{i-1/2}^n \bigr) + \Delta_-^{(3)} \check{\mathcal{Q}} \bigl( \boldsymbol{\phi}_i^n , \psi_{i+1/2}^n \bigr) . \end{align*}

    From Eq (5.8) we get

    \begin{align*} \bigl| \Delta_-^{(3)} \hat{\mathcal{Q}} \bigl( \boldsymbol{\phi}_i^n , \psi_{i-1/2}^n \bigr) \bigr| & = \Biggl| \eta' \bigl( \psi_{i-1/2}^n \bigr) \Delta_-^{(3)} \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i ) - \int_{\psi_{i-3/2}^n}^{\psi_{i-1/2}^n} \eta'' (u) \Bigl( \hat{\mathcal{F}} \bigl( \boldsymbol{\phi}_i^n, u \bigr) - \hat{\mathcal{F}} \bigl( \boldsymbol{\phi}_i^n , \psi_{i-1/2}^n \bigr) \Bigr) \, \mathrm{d} u \Biggr| \\ & \leq \bigl| \eta' \bigl( \psi_{i-1/2}^n \bigr) \bigr| \bigl| \Delta_-^{(3)} \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i ) \bigr| + \Biggl| \int_{\psi_{i-3/2}^n}^{\psi_{i-1/2}^n} \eta'' (u) \, \mathrm{d} u \Biggr| \bigl| \Delta_-^{(3)} \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i) \bigr| \\ & \leq 3 \| \eta' \|_{L^{\infty} (0, 1)} \bigl| \Delta_-^{(3)} \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i ) \bigr| \end{align*}

    and analogously

    \begin{align*} \bigl| \Delta_-^{(3)} \check{\mathcal{Q}} \bigl( \boldsymbol{\phi}_i^n , \psi_{i+1/2}^n \bigr) \bigr| \leq 3 \| \eta' \|_{L^{\infty} (0, 1)} \bigl| \Delta_-^{(4)} \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}_i^n ) \bigr|, \end{align*}

    hence

    \begin{align} \bigl| \Delta_-^{\psi} \mathcal{Q} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i ) \bigr| \leq 3 \| \eta' \|_{L^{\infty} (0, 1)} \bigl| \bigl( \Delta_-^{(3)} + \Delta_-^{(4)} \bigr) \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i ) \bigr|. \end{align} (5.20)

    On the other hand, we take into account that

    \begin{align*} \Delta_+^{\phi} \mathcal{Q}( \boldsymbol{\phi}_{i-1} ^n, \boldsymbol{\psi}_{i-1} ^n ) = \Delta_+^{\phi} \hat{\mathcal{Q}} \bigl( \boldsymbol{\phi}_{i-1}^n , \psi_{i-3/2}^n \bigr) + \Delta_+^{\phi} \check{\mathcal{Q}} \bigl( \boldsymbol{\phi}_{i-1}^n, \psi_{i-1/2}^n \bigr). \end{align*}

    Now

    \begin{array}{l} \Delta_+^{\phi} \hat{\mathcal{Q}} \bigl( \boldsymbol{\phi}_{i-1}^n, \psi_{i-3/2}^n \bigr) = \int_0^{\psi_{i-3/2}^n} \eta'(u) \bigl( \hat{f} \bigl( \boldsymbol{\phi}_i^n , u \bigr) - \hat{f} \bigl( \boldsymbol{\phi}_{i-1}^n , u \bigr) \bigr) \, \mathrm{d} u \\ = \Bigl[ \eta'(u) \bigl( \hat{\mathcal{F}} \bigl( \boldsymbol{\phi}_i^n, u \bigr) - \hat{\mathcal{F}} \bigl( \boldsymbol{\phi}_{i-1}^n, u \bigr) \bigr) \Bigr]_{u = 0}^{u = \psi_{i-3/2}^n} - \int_0^{\psi_{i-3/2}^n} \eta''(u) \bigl( \hat{\mathcal{F}} ( \boldsymbol{\phi}_i^n, u ) - \hat{\mathcal{F}} \bigl( \boldsymbol{\phi}_{i-1}^n , u ) \bigr) \, \mathrm{d} u \end{array} (5.21)

    and analogously

    \begin{array}{l} \Delta_+^{\phi} \check{\mathcal{Q}} \bigl( \boldsymbol{\phi}_{i-1}^n, \psi_{i-1/2}^n \bigr) = \int_0^{\psi_{i-1/2}^n} \eta'(u) \bigl( \check{f} \bigl( \boldsymbol{\phi}_i^n , v \bigr) - \check{f} \bigl( \boldsymbol{\phi}_{i-1}^n , v \bigr) \bigr) \, \mathrm{d} v \\ = \Bigl[ \eta'(v) \bigl( \check{\mathcal{F}} \bigl( \boldsymbol{\phi}_i^n, v \bigr) - \check{\mathcal{F}} \bigl( \boldsymbol{\phi}_{i-1}^n, v \bigr) \bigr) \Bigr]_{v = 0}^{v = \psi_{i-1/2}^n} - \int_0^{\psi_{i-1/2}^n} \eta''(v) \bigl( \check{\mathcal{F}} ( \boldsymbol{\phi}_i^n, v ) - \check{\mathcal{F}} ( \boldsymbol{\phi}_{i-1}^n , v ) \bigr) \, \mathrm{d} v. \end{array}

    Consequently,

    \begin{align} \bigl| \Delta_+^{\phi} \hat{\mathcal{Q}} \bigl( \boldsymbol{\phi}_{i-1}^n, \psi_{i-3/2}^n \bigr) \bigr| \leq 3 \| \eta' \|_{L^{\infty} (0, 1)} \max\limits_{0 \leq u \leq \psi_{i-3/2}^n } \bigl| \hat{\mathcal{F}} \bigl(\boldsymbol{\phi}_i^n, u \bigr) - \hat{\mathcal{F}} \bigl(\boldsymbol{\phi}_{i-1}^n, u \bigr) \bigr|, \end{align} (5.22)

    and by analogous reasoning for \check{\mathcal{Q}} ,

    \begin{align} \bigl| \Delta_+^{\phi} \check{\mathcal{Q}} \bigl( \boldsymbol{\phi}_{i-1}^n, \psi_{i-1/2}^n \bigr) \bigr| \leq 3 \| \eta' \|_{L^{\infty} (0, 1)} \max\limits_{0 \leq v \leq \psi_{i-1/2}^n } \bigl| \check{\mathcal{F}} \bigl( \boldsymbol{\phi}_i^n, v \bigr) - \check{\mathcal{F}} \bigl(\boldsymbol{\phi}_{i-1}^n, v \bigr) \bigr|. \end{align} (5.23)

    To estimate the right-hand side of Eq (5.22), we recall that

    \begin{align*} \hat{\mathcal{F}} \bigl(\boldsymbol{\phi}_i^n, u \bigr) - \hat{\mathcal{F}} \bigl(\boldsymbol{\phi}_{i-1}^n, u \bigr) = \int_0^u \bigl( ((f_i^n)')^+ (s) - ((f_{i-1}^n)')^+ (s) \bigr) \, \mathrm{d} s = : \hat{\mathcal{D}}_{i-1/2}^n. \end{align*}

    We now assume that \sigma = 1 and use Eq (3.9). To estimate \smash{ \hat{\mathcal{D}}_{i-1/2}^n} , we will appeal to the trivial but useful inequality | (\alpha \vee x) - (\beta \vee y) | \leq |\alpha - \beta| + |x-y| . We proceed by discussing all possible cases of the location of u in relation to the extrema \smash{\hat{\psi}_i^n} and \smash{\hat{\psi}_{i-1}^n} of f_i^n and f_{i-1}^n , respectively, and assume that \sigma = 1 .

    1. Assume that \smash{u \leq \hat{\psi}_i^n \wedge \hat{\psi}_{i-1}^{n} } . Since \smash{ \hat{\psi}_i^n \leq \psi_{\max, i}^n } and \smash{ \hat{\psi}_{i-1}^{n} \leq \psi_{\max, i-1}^{n} } , in this case \smash{\hat{\mathcal{D}}}_{i-1/2}^n = 0 if \smash{\psi_{\max, i}^n = \psi_{\max, i-1}^{n} } and otherwise

    \begin{align*} \bigl| \hat{\mathcal{D}}_{i-1/2}^n \bigr| & = \bigl| f_i^n (u) - f_{i-1}^n (u) \bigr| = u \bigl| \tilde{V} \bigl( u/ \psi_{\max, i}^n \bigr) - \tilde{V} \bigl( u/ \psi_{\max, i-1}^{n} \bigr) \bigr| \\ & \leq \frac{u^2}{\psi_{\max, i}^n \psi_{\max, i-1}^{n}} \| \tilde{V}' \|_{\infty} | \psi_{\max, i}^n - \psi_{\max, i-1}^{n}| \leq \| \tilde{V}' \|_{\infty} \bigl| \psi_{\max, i}^n - \psi_{\max, i-1}^{n} \bigr|. \end{align*}

    Noticing that

    \begin{align} \bigl| \psi_{\max, i}^n - \psi_{\max, i-1}^{n} \bigr| & = \bigl| \phi_{i-1/2}^{n} \vee \phi_{i+1/2}^{n} - \phi_{i-3/2}^{n} \vee \phi_{i-1/2}^{n} \bigr| \leq \bigl| \Delta_- \phi_{i-1/2}^{n} \bigr| + \bigl| \Delta_+ \phi_{i-1/2}^n \bigr|, \end{align} (5.24)

    we conclude that

    \begin{align*} \bigl| \hat{\mathcal{D}}_{i-1/2}^n \bigr| \leq \| \tilde{V}' \|_{\infty} \bigl( \bigl| \Delta_- \phi_{i-1/2}^{n} \bigr| + \bigl| \Delta_+ \phi_{i-1/2}^n \bigr| \bigr). \end{align*}

    2. If \smash{\hat{\psi}_i^n < u < \hat{\psi}_{i-1}^{n}} then

    \begin{align*} \bigl| \hat{\mathcal{D}}_{i-1/2}^n \bigr| & = \bigl| f_i^n \bigl( \hat{\psi}_i^n ) - f_{i-1}^n (u) \bigr| \leq \bigl| f_i^n \bigl( \hat{\psi}_i^n \bigr) - f_{i-1}^n \bigl( \hat{\psi}_{i-1}^n \bigr) \bigr| + \bigl| f_{i-1}^n \bigl( \hat{\psi}_{i-1}^n \bigr) - f_{i-1}^n (u) \bigr|. \end{align*}

    Since \smash{f_i^n (\hat{\psi}_i^n) = \psi_{\max, i}^n \omega \tilde{V} (\omega)} for all i , we conclude that

    \begin{align} \bigl| f_i^n \bigl( \hat{\psi}_i^n \bigr) - f_{i-1}^n \bigl( \hat{\psi}_{i-1}^n \bigr) \bigr| & \leq \omega \| \tilde{V} \|_{\infty} \left| \psi_{\max, i}^n - \psi_{\max, i-1}^n \right| \leq \| \tilde{V} \|_{\infty} \bigl( \bigl| \Delta_- \phi_{i-1/2}^{n} \bigr| + \bigl| \Delta_+ \phi_{i-1/2}^n \bigr| \bigr). \end{align} (5.25)

    On the other hand, in the present case

    \begin{align*} \bigl| f_{i-1}^n \bigl( \hat{\psi}_{i-1}^n \bigr) - f_{i-1}^n (u) \bigr| = f_{i-1}^n \bigl( \hat{\psi}_{i-1}^n \bigr) - f_{i-1}^n (u) \leq f_{i-1}^n \bigl( \hat{\psi}_{i-1}^n \bigr) - f_{i-1}^n \bigl( \hat{\psi}_{i}^n \bigr). \end{align*}

    Since for \smash{s \in [0, \hat{\psi}_{i-1}^n]} there holds \smash{(f_{i-1}^n)'(s) \leq (f_{i-1}^n)' (0) = \tilde{V} (0)} , we get

    \begin{align*} \bigl| f_{i-1}^n \bigl( \hat{\psi}_{i-1}^n \bigr) - f_{i-1}^n (u) \bigr| & = \int_{u}^{\hat{\psi}_{i-1}^n} (f_{i-1}^n)'(s) \, \mathrm{d} s \leq \int_{\hat{\psi}_i^n}^{\hat{\psi}_{i-1}^n} (f_{i-1}^n)'(s) \, \mathrm{d} s \leq \tilde{V}(0) \bigl(\hat{\psi}_{i-1}^n - \hat{\psi}_i^n \bigr). \end{align*}

    Lemma 3.1 (a) implies that

    \begin{align} \bigl| \hat{\psi}_i^n - \hat{\psi}_{i-1}^{n} \bigr| = \omega \bigl| \psi_{\max, i}^n - \psi_{\max, i-1}^{n} \bigr|, \end{align} (5.26)

    hence

    \begin{align*} \bigl| \hat{\mathcal{D}}_{i-1/2}^n \bigr| \leq 2 \| \tilde{V} \|_{\infty} \bigl( \bigl| \Delta_- \phi_{i-1/2}^n \bigr| + \bigl| \Delta_+ \phi_{i-1/2}^{n} \bigr| \bigr). \end{align*}

    The same estimate holds for \smash{\hat{\psi}_{i-1}^{n} < u < \hat{\psi}_i^{n}} .

    3. If \smash{u \geq \hat{\psi}_i^n \vee \hat{\psi}_{i-1}^{n}} then utilizing Eq (5.25) we get

    \begin{align*} \bigl| \hat{\mathcal{D}}_{i-1/2}^n \bigr| = \bigl| f_i^n \bigl( \hat{\psi}_i^n \bigr) - f_{i-1}^n \bigl( \hat{\psi}_{i-1}^n \bigr) \bigr| \leq \| \tilde{V} \|_{\infty} \bigl( \bigl| \Delta_- \phi_{i-1/2}^{n} \bigr| + \bigl| \Delta_+ \phi_{i-1/2}^n \bigr| \bigr). \end{align*}

    Combining all possible cases we deduce that

    \begin{align} \bigl| \hat{\mathcal{D}}_{i-1/2}^n \bigr| \leq \bigl( \| \tilde{V}' \|_{\infty} + 2 \| \tilde{V} \|_{\infty} \bigr) \bigl( \bigl| \Delta_- \phi_{i-1/2}^{n} \bigr| + \bigl| \Delta_+ \phi_{i-1/2}^n \bigr| \bigr). \end{align} (5.27)

    Next, we deal with Eq (5.23), recalling that (see Eq (5.4))

    \begin{align*} & \check{\mathcal{F}} ( \boldsymbol{\phi}_i^n, v ) - \check{\mathcal{F}} (\boldsymbol{\phi}_{i-1}^n, v ) = \check{\mathcal{D}}_{i-1/2}^n + \mathcal{E}_{i-1/2}^n, \end{align*}

    where we define

    \begin{align*} \check{\mathcal{D}}_{i-1/2}^n & : = \int_0^v \bigl( ((f_i^n)')^- (s) - ((f_{i-1}^n)')^- (s) \bigr) \, \mathrm{d} s, \quad \mathcal{E}_{i-1/2}^n : = \bigl( \phi_{i-1/2}^n \tilde{v}_{\phi \mathrm{s}} ( \phi_{i+1/2}^n) - \phi_{i-3/2}^n \tilde{v}_{\phi \mathrm{s}} ( \phi_{i-1/2}^n ) \bigr) v. \end{align*}

    The discussion of all possible cases of the location of v in relation to \smash{\hat{\psi}_i^n} and \smash{\hat{\psi}_{i-1}^n} gives rise to the following cases for the estimation of \smash{\check{\mathcal{D}}_{i-1/2}^n} .

    1. If v \leq \hat{\psi}_{i}^n \wedge \hat{\psi}_{i-1}^{n} or v \geq \psi_{\max, i}^n \vee \psi_{\max, i-1}^{n} then \smash{\check{\mathcal{D}}_{i-1/2}^n = 0} .

    2. To handle the case \smash{ \hat{\psi}_i^n \wedge \hat{\psi}_{i-1}^{n} \leq v \leq \hat{\psi}_i^n \vee \hat{\psi}_{i-1}^{n} } we assume that \smash{ \hat{\psi}_{i}^{n} < \hat{\psi}_{ i-1}^{n} } and \smash{ \hat{\psi}_{i}^{n} \leq v \leq \hat{\psi}_{ i-1}^{n} } . Then

    \begin{align*} \bigl| \check{\mathcal{D}}_{i-1/2}^n \bigr| & = \bigl| f_i^n (v) - f_i^n \bigl( \hat{\psi}_i^n \bigr) \bigr| = \left| \int^v_{\hat{\psi_i^n}} (( f_i^n)')^- (s) \, \mathrm{d} s \right| \leq \max\limits_{\hat{\psi}_i^n \leq s \leq \hat{\psi}_{i-1}^n } \bigl| (( f_i^n)')^- (s) \bigr| \bigl| \hat{\psi}_i^n - \hat{\psi}_{i-1}^{n} \bigr| \\ & \leq \bigl| (f_i^n)' ( \psi_{\mathrm{infl}, i}^n ) \bigr| \bigl| \hat{\psi}_i^n - \hat{\psi}_{i-1}^{n} \bigr|. \end{align*}

    By Lemma 3.1 (ⅱ),

    \begin{align*} (f_i^n)' ( \psi_{\mathrm{infl}, i}^n ) = (f_i^n)' ( \tilde{\omega} \phi_{\max, i}^n ) = \tilde{V} ( \tilde{\omega}) + \tilde{\omega} \tilde{V}' ( \tilde{\omega}), \end{align*}

    so \smash{(f_i^n)' (\psi_{\mathrm{infl}, i}^n)} does not depend on \smash{\phi_{\max, i}^n} and we conclude that

    \begin{align*} \bigl| \check{\mathcal{D}}_{i-1/2}^n \bigr| \leq \bigl( \| \tilde{V} \|_{\infty} + \| \tilde{V}' \|_{\infty} \bigr) \bigl| \hat{\psi}_i^n - \hat{\psi}_{i-1}^{n} \bigr|. \end{align*}

    Applying the argument of Eq (5.24) and Eq (5.26) yields

    \begin{align*} \bigl| \check{\mathcal{D}}_{i-1/2}^n \bigr| \leq \bigl( \| \tilde{V} \|_{\infty} + \| \tilde{V}' \|_{\infty} \bigr) \bigl( \bigl| \Delta_- \phi_{i-1/2}^n \bigr| + \bigl|\Delta_+ \phi_{i-1/2}^{n} \bigr| \bigr). \end{align*}

    The same inequality is also valid if \smash{ \hat{\psi}_{i-1}^{n} < \hat{\psi}_{ i}^{n} } and \smash{ \hat{\psi}_{i-1}^{n} \leq v \leq \hat{\psi}_{ i}^{n} } .

    3. Finally, assume that \smash{v > \hat{\psi}_{i}^n \vee \hat{\psi}_{i-1}^{n} } . In this case

    \begin{align*} \bigl| \check{\mathcal{D}}_{i-1/2}^n \bigr| & = \bigl| f_i^n (v) - f_{i-1}^n (v) - f_i^n \bigl( \hat{\psi}_{i}^n \bigr) + f_{i-1}^n \bigl( \hat{\psi}_{i-1}^n \bigr) \bigr| \leq \bigl| f_i^n (v) - f_{i-1}^n (v) \bigr| + \bigl| f_i^n \bigl( \hat{\psi}_{i}^n \bigr) - f_{i-1}^n \bigl( \hat{\psi}_{i-1}^n \bigr) \bigr|. \end{align*}

    Taking into account that f_i^n (\hat{\psi}_{i}^n) = \psi_{\max, i}^n \omega \tilde{V} (\omega) , we get

    \begin{align} \bigl| f_i^n \bigl( \hat{\psi}_{i}^n \bigr) - f_{i-1}^n \bigl( \hat{\psi}_{i-1}^n \bigr) \bigr|& = \omega \tilde{V} ( \omega) \left| \psi_{\max, i}^n - \psi_{\max, i-1}^n \right| \leq \| \tilde{V} \|_{\infty} \bigl( \bigl| \Delta_- \phi_{i-1/2}^n \bigr| + \bigl|\Delta_+ \phi_{i-1/2}^{n} \bigr| \bigr). \end{align} (5.28)

    If v > \psi_{\max, i}^n \vee \psi_{\max, i-1}^{n} , then f_i^n(v) = f_{i-1}^n(v) = 0 , hence Eq (5.28) means that

    \begin{align*} \bigl| \check{\mathcal{D}}_{i-1/2}^n \bigr| \leq \| \tilde{V} \|_{\infty} \bigl( \bigl|\Delta_- \phi_{i-1/2}^n \bigr| + \bigl| \Delta_+ \phi_{i-1/2}^{n} \bigr| \bigr). \end{align*}

    Suppose now that

    \begin{align} \hat{\psi}_{i}^n \vee \hat{\psi}_{i-1}^{n} \leq \psi_{\max, i}^n \leq v \leq \psi_{\max, i-1}^{n}. \end{align} (5.29)

    Since we know that \smash{v = \psi_{i+1/2}^n \leq 1-\phi_{i+1/2}^n} , the inequality \smash{ \psi_{\max, i}^n \leq v} can only be satisfied if

    \begin{align*} \psi_{\max, i}^n = (1- \phi_{i-1/2}^n) \wedge (1- \phi_{i+1/2}^n) = 1 - \phi_{i-1/2}^n. \end{align*}

    On the other hand,

    \begin{align*} \psi_{\max, i-1}^n = (1 - \phi_{i-1/2}^n) \wedge (1-\phi_{i-3/2}^n) \leq 1- \phi_{i-1/2}^n, \end{align*}

    so Eq (5.29) is only possible when \smash{\psi_{\max, i}^n = v = \psi_{\max, i-1}^{n} = 1- \phi_{i-1/2}^n} , which means \smash{ \check{\mathcal{D}}_{i-1/2}^n = 0} .

    If instead of Eq (5.29) we have

    \begin{align} \hat{\psi}_{i}^n \vee \hat{\psi}_{i-1}^{n} \leq \psi_{\max, i-1}^n \leq v \leq \psi_{\max, i}^{n}, \end{align} (5.30)

    then \smash{1/\psi_{\max, i-1}^n \leq 1/(\omega \psi_{\max, i}^n)} implies

    \begin{align*} \bigl| f_i^n (v) - f_{i-1}^n (v) \bigr| & = v \bigl| \tilde{V} \bigl( v/ \psi_{\max, i}^n \bigr) - \tilde{V} \bigl( v/ \psi_{\max, i-1}^{n} \bigr) \bigr| \\ & \leq v^2 \| \tilde{V}' \|_{\infty} \frac{| \psi_{\max, i}^n - \psi_{\max, i-1}^{n}|}{\psi_{\max, i}^n \psi_{\max, i-1}^{n}} \leq \frac{ \| \tilde{V}' \|_{\infty} }{\omega} \bigl| \psi_{\max, i}^n - \psi_{\max, i-1}^{n} \bigr|. \end{align*}

    The remainder of the estimate is based on Eq (5.24). Since \omega < 1 , we conclude that if Eq (5.29) holds, then

    \begin{align*} \bigl| f_i^n (v) - f_{i-1}^n (v) \bigr| \leq \frac{ \| \tilde{V}' \|_{\infty}}{\omega} \bigl( \bigl| \Delta_- \phi_{i-1/2}^n \bigr| + \bigl| \Delta_+ \phi_{i-1/2}^{n} \bigr| \bigr). \end{align*}

    In combination with Eq (5.28) we obtain in this case

    \begin{align} \bigl| \check{\mathcal{D}}_{i-1/2}^n \bigr| \leq \biggl( \| \tilde{V} \|_{\infty} + \frac{ \| \tilde{V}' \|_{\infty}}{\omega} \biggr) \bigl( \bigl| \Delta_- \phi_{i-1/2}^n \bigr| + \bigl| \Delta_+ \phi_{i-1/2}^{n} \bigr| \bigr). \end{align} (5.31)

    Next, suppose that instead of Eq (5.29) or Eq (5.30) there holds

    \begin{align*} \hat{\psi}_{i}^n \leq \psi_{\max, i}^n \leq \hat{\psi}_{i-1}^{n} \leq v \leq \psi_{\max, i-1}^{n}, \end{align*}

    then the discussion of Eq (5.29) can be applied again and we get that this ordering is only feasible if all terms are equal and zero, and therefore \smash{ \check{\mathcal{D}}_{i-1/2}^n = 0} . On the other hand, let us assume that

    \begin{align*} \hat{\psi}_{i-1}^n \leq \psi_{\max, i-1}^n \leq \hat{\psi}_{i}^{n} \leq v \leq \psi_{\max, i}^{n}. \end{align*}

    In this case,

    \begin{align*} \bigl| \check{\mathcal{D}}_{i-1/2}^n \bigr| & = \left| \int_0^v ((f_i^n)')^- (s) \, \mathrm{d} s \right| = \left| \int_{\hat{\psi}_i^n}^v ((f_i^n)')^- (s) \, \mathrm{d} s \right| \leq \int_{\hat{\psi}_i^n}^v \bigl| ((f_i^n)')^- (s) \bigr| \, \mathrm{d} s \\ & \leq \int_{\psi_{\max, i-1}^n}^{\psi_{\max, i}^n} \bigl| (f_i^n)' \bigr| \, \mathrm{d} s \leq \bigl( \| \tilde{V} \|_{\infty} + \| \tilde{V}' \|_{\infty} \bigr) \bigl( \psi_{\max, i}^{n} - \psi_{\max, i-1}^{n}\bigr) \\ & \leq \bigl( \| \tilde{V} \|_{\infty} + \| \tilde{V}' \|_{\infty} \bigr) \bigl( \bigl|\Delta_- \phi_{i-1/2}^n \bigr| + \bigl|\Delta_+ \phi_{i-1/2}^{n} \bigr| \bigr). \end{align*}

    It remains to treat the case \hat{\psi}_i^n \vee \hat{\psi}_{i-1}^{n} \leq v \leq \psi_{\max, i}^n \wedge \psi_{\max, i-1}^{n} . We then have \smash{v^2/(\psi_{\max, i-1}^{n} \psi_{\max, i}^n) \leq 1} , and analogously to the derivation of Eq (5.31) we get

    \begin{align*} \bigl| \check{\mathcal{D}}_{i-1/2}^n \bigr| \leq \bigl( \| \tilde{V} \|_{\infty} +\| \tilde{V}' \|_{\infty} \bigr) \bigl( \bigl|\Delta_- \phi_{i-1/2}^n \bigr| + \bigl|\Delta_+ \phi_{i-1/2}^{n} \bigr| \bigr). \end{align*}

    Collecting all estimates for \smash{ \check{\mathcal{D}}_{i-1/2}^n } , we see that

    \begin{align} \bigl|\check{\mathcal{D}}_{i-1/2}^{n} \bigr| & \leq \bigl( 3 \| \tilde{V} \|_{\infty} + (1+1/\omega) \| \tilde{V}' \|_{\infty} \bigr) \bigl( \bigl|\Delta_- \phi_{i-1/2}^n \bigr| + \bigl| \Delta_+ \phi_{i-1/2}^{n} \bigr| \bigr). \end{align} (5.32)

    Furthermore, we obtain

    \begin{align} \bigl| \mathcal{E}_{i-1/2}^{n} \bigr| \leq \| \tilde{v}_{\phi \mathrm{s}} '\|_{\infty} \bigl| \Delta_+ \phi_{i-1/2}^{n}\bigr| + \| \tilde{v}_{\phi \mathrm{s}} \|_{\infty} \bigl| \Delta_- \phi_{i-1/2}^n \bigr|. \end{align} (5.33)

    Combining the estimates Eq (5.27), Eq (5.32) and Eq (5.33), we obtain from Eq (5.22) and Eq (5.23) the bounds

    \begin{align} \bigl| \hat{\mathcal{F}} \bigl(\boldsymbol{\phi}_i^n, u \bigr) - \hat{\mathcal{F}} \bigl(\boldsymbol{\phi}_{i-1}^n, u \bigr) \bigr|, \, \bigl| \check{\mathcal{F}} \bigl(\boldsymbol{\phi}_i^n, v \bigr) - \check{\mathcal{F}} \bigl(\boldsymbol{\phi}_{i-1}^n, v \bigr) \bigr| \leq C_{10} \bigl( \bigl| \Delta_- \phi_{i-1/2}^n \bigr| + \bigl| \Delta_+ \phi_{i-1/2}^{n} \bigr| \bigr) \end{align} (5.34)

    and therefore

    \begin{align*} \bigl| \Delta_-^{\phi} \mathcal{Q} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i ) \bigr| \leq C_{11} \bigl( \bigl| \Delta_- \phi_{i-1/2}^n \bigr| + \bigl| \Delta_+ \phi_{i-1/2}^{n} \bigr| \bigr) \end{align*}

    with constants C_{10} and C_{11} . Combining the last inequality with Eq (5.19) and Eq (5.20) we arrive at the desired estimate Eq (5.18).

    From Eq (5.18), and considering that 0 \leq \phi_{i-1/2}^n \leq 1 for all i and n , we obtain

    \begin{array}{l} \bigl( \Delta_- Q ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i ) \bigr)^2 \leq 2 C_8^2 \bigl( \bigl( \Delta_-^{(3)} + \Delta_-^{(4)} \bigr) \mathcal{F} ( \boldsymbol{\phi}_i^n, \boldsymbol{\psi}_i^n \bigr) \bigr)^2 + 2 C_9^2 \bigl( \bigl| \Delta_- \phi_{i-1/2}^n \bigr| + \bigl| \Delta_+ \phi_{i-1/2}^{n} \bigr| \bigr)^2 \\ \leq 4 C_8^2 \Bigl( \bigl( \Delta_-^{(3)} \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i ) \bigr)^2 + \bigl( \Delta_-^{(4)} \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i ) \bigr)^2 \Bigr) + 4 C_9^2 \bigl( \bigl| \Delta_- \phi_{i-1/2}^n \bigr|^2 + \bigl| \Delta_+ \phi_{i-1/2}^{n} \bigr|^2 \bigr) \\ \leq C_{12} \Bigl( \bigl( \Delta_-^{(3)} \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i \bigr) \bigr)^2 + \bigl( \Delta_-^{(4)} \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i \bigr) \bigr)^2 + \bigl| \Delta_- \phi_{i-1/2}^n \bigr| + \bigl| \Delta_+ \phi_{i-1/2}^{n} \bigr| \Bigr). \end{array}

    Summing over (i, n) \in \mathcal{I}_0 we get

    \begin{align*} \sum\limits_{\mathcal{I}_0} \bigl( \Delta_- Q( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i ) \bigr)^2 & \leq C_{12} \sum\limits_{\mathcal{I}_0} \Bigl( \bigl( \Delta_-^{(3)} \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i ) \bigr)^2 + \bigl( \Delta_-^{(4)} \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i ) \bigr) \Bigr)^2 + 2 C_{12} \sum\limits_{ \mathcal{I}_1 } \bigl| \Delta_+ \phi_{i-1/2}^n \bigr| \\ & \leq C_{12} \sum\limits_{\mathcal{I}_0} \Bigl( \bigl( \Delta_-^{(3)} \mathcal{F} ( \boldsymbol{\phi}_i^n, \boldsymbol{\psi}_i^n \bigr) \bigr)^2 + \bigl( \Delta_-^{(4)} \mathcal{F} ( \boldsymbol{\phi}_i^n, \boldsymbol{\psi}_i^n \bigr) \Bigr)^2 + 2 C_{12} \sum\limits_{n = 0}^{N} \operatorname*{TV} ( \boldsymbol{\phi}^n). \end{align*}

    Multiplying this inequality by \Delta t \Delta z and taking into account Eq (5.16) and the uniform bound on \operatorname*{TV} (\boldsymbol{\phi}^n) we have proved the following lemma.

    Lemma 5.4. Consider numerical approximations produced by Scheme 3. There exists a constant C = C(T) that is independent of \Delta t and \Delta z such that the following estimate holds:

    \begin{align*} \Delta t \Delta z \sum\limits_{n = 0}^{N_T} \sum\limits_{i \in \mathbb{Z}} \bigl( \Delta_- \mathcal{Q} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i ) \bigr)^2 \leq C(T) \Delta z. \end{align*}

    In part following the proofs of Lemmas 5.5 and 5.9 in [33] and in particular that of Lemma 3.4 in [18] we now prove the H_{\mathrm{loc}}^{-1} compactness result.

    Lemma 5.5. Assume that \psi^{\Delta z} is generated by the scheme Eq (5.3) ( \psi -Scheme 3), and that \phi is the unique entropy solution of Eq (1.8a), Eq (1.2a) on \Pi_T . Furthermore, we denote by (\eta_0, Q_0) the Kružkov entropy pair Eq (5.2), and the distribution

    \begin{align*} \mu^{\Delta z} : = \partial_t \eta_0 ( \psi^{\Delta z} ) + \partial_z Q_0 ( \phi^{\Delta z} , \psi^{\Delta z}). \end{align*}

    Then the sequence \{ \mu^{\Delta z} \}_{\Delta z > 0} belongs to a compact subset of H_{\mathrm{loc}}^{-1} (\Pi_T) .

    The proof of Lemma 5.5 essentially follows the proof of [18,Lemma 3.4], with some slight modifications due the particular definition of the numerical flux in the present case. The proof is presented in Appendix A.

    Since \phi^{\Delta z} \to \phi strongly in L^p , we obtain that there exists a constant C such that

    \begin{align*} \bigl| \big\langle \partial_z \bigl( Q ( \phi^{\Delta z}, \psi^{\Delta z} ) - Q ( \phi, \psi^{\Delta z} ) \bigr) , \zeta \big\rangle \bigr| & \leq C \| \phi^{\Delta z} - \phi \|_{L^2( \Pi_T)} \| \zeta \|_{H^1 ( \Pi_T) } \to 0 \quad \text{as }\Delta z \to 0 , \end{align*}

    hence the sequence \smash{ \{ \tilde{\mu}^{\Delta z} \}_{\Delta z > 0}} , where we define

    \begin{align*} \tilde{\mu}^{\Delta z} : = \partial_t \eta_0 ( \psi^{\Delta z} ) + \partial_z Q ( \phi, \psi^{\Delta z} ), \end{align*}

    is compact in \smash{H_{\mathrm{loc}}^{-1} (\Pi_T)} . Now, by Lemma 5.1 there exists a subsequence \{ \psi^{\Delta z} \} (which we do not relabel) and a function \psi \in L^{\infty} (\Pi_T) such that

    \begin{align} \psi^{\Delta z} \to \psi \quad \text{as}\; \Delta z \to 0 , \; \text{a.e. and in } L^p_{\mathrm{loc}} ( \Pi_T)\; \text{for any } p \in [1, \infty) . \end{align} (5.35)

    Theorem 5.1. Assume that the maps \phi and \psi are the limit functions of \phi^{\Delta z} and of \psi^{\Delta z} as \Delta z \to 0 (the latter one being defined by Eq (5.35), that is, we consider Scheme 3). Then (\phi, \psi) is a weak solution of the initial-value problem Eq (1.8), Eq (1.2) in the sense of Definition 1.1.

    The proof follows that of [18,Lemma 3.5], again with slight modifications. We refer to Appendix A for details.

    To simplify error estimations we utilize a grid with the property that the boundaries of the tank agree with the boundaries of a cell (see Figure 2 (left)). Since an exact solution is frequently difficult to obtain, we use an approximate reference solution obtained with a large number N_\mathrm{ref} cells against which the error of other simulated solutions with N < N_\mathrm{ref} is measured. The error is estimated on a fixed interval [0, z_\mathrm{end}) slightly larger than the column of height H so that the outflow volume fractions are included. We define the coarsest grid of N_0 cells with \Delta z_0: = H/(N_0-2) and place the column between z_\mathrm{U}: = \Delta z_0 and z_\mathrm{E}: = z_\mathrm{U}+H = (N_0-1)\Delta z_0 . This corresponds to Figure 2 (left) with N = N_0 . We define the length of the interval of error estimation as z_\mathrm{end}: = H+2\Delta z_0 = N_0\Delta z_0 .

    To estimate the convergence order, we simulate with N_k = 2^kN_0 cells, k = 0, \ldots, k_\mathrm{ref}-1 , where the integer k_\mathrm{ref} defines the number of cells N_\mathrm{ref}: = N_{k_\mathrm{ref}} : = 2^{k_\mathrm{ref}}N_0 of the reference solution. Then we define \Delta z_k: = z_\mathrm{end}/N_k , \Delta z^\mathrm{r}: = \Delta z_{k_\mathrm{ref}} : = z_\mathrm{end}/N_{k_\mathrm{ref}} = \Delta z_{0}/2^{k_\mathrm{ref}} and the factor of refinement from N_k cells to N_{\mathrm{ref}} as m_k: = \Delta z_k/\Delta z^\mathrm{r} = N_{k_\mathrm{ref}}/N_k = 2^{k_\mathrm{ref}-k} . We note that z_{N_k}: = N_k\Delta z_k = z_\mathrm{end} for all k .

    We will now measure the error between the piecewise constant numerical solution obtained by N = N_k cells (we skip the index k for a moment) and the reference solution obtained with N_\mathrm{ref} cells on the grid refined by a factor m = \Delta z/\Delta z^\mathrm{r} . The refined grid satisfies z_0^\mathrm{r}: = z_0 = 0 and we have z_i = i\Delta z = im\Delta z^\mathrm{r} = :z^\mathrm{r}_{im} . The corresponding numerical solutions on the refined grids are denoted by (skipping the time index n ) \smash{\phi^\mathrm{r}_{i+1/2}} , \smash{\psi^\mathrm{r}_{i+1/2}} , etc., where \smash{A^\mathrm{r}_{i+1/2}} are defined by \Delta z^\mathrm{r} . The refined cells are numbered such that the first cell for \phi above z^\mathrm{r}_0 = 0 contains the value \smash{\phi^\mathrm{r}_{1/2}} . Then z_\mathrm{end} = Nm\Delta z^\mathrm{r} . This means that the cells within [0, z_\mathrm{end}) contain the values \smash{\phi^\mathrm{r}_{1/2}}, \dots, \smash{\phi^\mathrm{r}_{Nm-1/2}} , and analogously for \psi ; see Figure 2 (right).

    Note that the location of the spatial discontinuities z_\mathrm{U} and z_\mathrm{E} will coincide with a cell boundary for any mesh considered in the refinement process while the locations of the inlets z_\mathrm{F, 1} , etc. will be chosen in such a way that each of them lies inside a cell for the finest grid; hence, they do this for all the coarser meshes. In this way, the numerical fluxes at cell boundaries are well defined.

    We compute the estimated error at a time point t = T and define

    \begin{align*} \bigl\|\phi^{\Delta z}(\cdot, T)\bigr\| : = \int_{0}^{z_\mathrm{end}} A(z)\bigl|\phi^{\Delta z}(z, T)\bigr|\, \mathrm{d}z. \end{align*}

    The L^1 -difference between two numerical solutions computed on grids with cell sizes \Delta z and \Delta z^\mathrm{r} is calculated as follows for \phi :

    \begin{align*} E^\phi(\Delta z, \Delta z^\mathrm{r}, T) : = \bigl\|\phi^{\Delta z}(\cdot, t)-\phi^{\Delta z^\mathrm{r}}(\cdot, T)\bigr\| = \sum\limits_{i = 0}^{N-1} I^\phi_{i+1/2}(T) \end{align*}

    with the summands defined by

    \begin{align*} I^\phi_{i+1/2}(T)&: = \int_{z_i}^{z_{i+1}}A(z) \bigl|\phi^{\Delta z}(z, T) -\phi^{\Delta z^\mathrm{r}}(z, T)\bigr|\, \mathrm{d}z = \sum\limits_{k = 0}^{m-1}\int_{z^\mathrm{r}_{im+k}}^{z^\mathrm{r}_{im+k+1}}A(z)\bigl|\phi_{i+1/2} -\phi^{\mathrm{r}}_{im+1/2+k}\bigr|\, \mathrm{d}z\\ & = \Delta z^{\mathrm{r}} \sum\limits_{k = 0}^{m-1}A^\mathrm{r}_{im+1/2+k}\bigl|\phi_{i+1/2} -\phi^{\mathrm{r}}_{im+1/2+k}\bigr|. \end{align*}

    The approximate relative error for \phi in the interval [0, z_\mathrm{end}) is then defined as

    \begin{align*} e_{N_k}^{\phi}(T): = \frac{E^\phi(\Delta z_k, \Delta z^\mathrm{r}, T) }{\|\phi^{\Delta z^\mathrm{r}}(\cdot, T)\|} = \frac{ \|\phi^{\Delta z_k}(\cdot, T)-\phi^{\Delta z^\mathrm{r}}(\cdot, T)\|} {\|\phi^{\Delta z^\mathrm{r}}(\cdot, T)\|}. \end{align*}

    We define e_N^{\psi}(t) analogously and hence, the total relative error can be defined as

    \begin{align*} e^\mathrm{tot}_{N_k}(T): = e_{N_k}^{\phi}(T)+e_{N_k}^{\psi}(T) \end{align*}

    and the observed convergence order between two discretizations N_{k-1} and N_{k} as

    \begin{equation*} \Upsilon_k (T) : = - \dfrac{\log(e^\mathrm{tot}_{N_{k-1}}(T)/e^\mathrm{tot}_{N_{k}}(T))}{\log(N_{k-1}/N_k)}, \quad k = 1, \dots, k_{\mathrm{ref}}-1. \end{equation*}

    For smooth solutions and a constant A (see Eq (1.7)), we also use an alternative way [7] of calculating approximate errors and convergence orders in which a reference solution is not needed. One can use cubic interpolation to compute the quantities \smash{\tilde{\phi}_{i+1/2}^{{\Delta z}_{k}}} from a grid with N_{k+1} = 2\, N_k cells, k = 0, \ldots, \hat{k} , with \hat{k} an integer, taking into consideration that \smash{z^k_{i+1/2} = (z^{k+1}_{2i+1/2}+z^{k+1}_{2i+3/2})/2} . Then, \smash{\tilde{\phi}_{i+1/2}^{{\Delta z}_k}} is given by

    \begin{align*} \tilde{\phi}_{i+1/2}^{{\Delta z}_k} = \frac{9}{16}\bigl(\phi_{2i+3/2}^{{\Delta z}_{k+1}}+ \phi_{2i+1/2}^{{\Delta z}_{k+1}} \bigr) - \frac{1}{16}\bigl( \phi_{2i+5/2}^{{\Delta z}_{k+1}}+\phi_{2i-1/2}^{{\Delta z}_{k+1}}\bigr), \quad i = 1, \ldots, N_k. \end{align*}

    The alternative approximate relative L^1 -error for \phi can then be calculated as

    \begin{align*} \hat{e}_{N_k}^{\phi}(T): = \frac{1}{N_k}\sum\limits_{i = 1}^{N_k}\bigl| \tilde{\phi}_{i+1/2}^{{\Delta z}_k}(\cdot, T) -\phi_{i+1/2}^{{\Delta z}_k}(\cdot, T) \bigr|. \end{align*}

    We can define \smash{\tilde{\psi}_{i+1/2}^{N_k}} and \smash{\hat{e}_{N_k}^{\psi}(T)} analogously along with the alternative total approximate L^1 -error and convergence order

    \begin{align*} \hat{e}^\mathrm{tot}_{N_k}(T) : = \hat{e}_{N_k}^{\phi}(T)+\hat{e}_{N_k}^{\psi}(T), \quad \hat{\Upsilon}_k (T) : = \log_2(\hat{e}^\mathrm{tot}_{N_{k}}(T)/\hat{e}^\mathrm{tot}_{N_{k+1}}(T)) \quad \text{for } k = 0, \dots, \hat{k} . \end{align*}

    For the first example, in Section 6.3, we use a smooth solution away from spatial discontinuities, to estimate the order of convergence of the numerical scheme. For this example, we use N_0 = 500 , N_k = 2^k N_0 for k = 0, 1, \ldots, 5 and k_{\mathrm{ref}} = 8 ; hence, N_5 = 16\, 000 and N_{\mathrm{ref}} = 128\, 000 .

    In Sections 6.4 and 6.5, we exemplify counter- and co-current flows of the primary and secondary disperse phases, respectively. For these two examples, we use N_0 = 100 , and k_{\mathrm{ref}} = 7 . We set three inlets z_{\mathrm{F}, 1} , z_{\mathrm{F}, 2} and z_{\mathrm{F}, 3} dividing the tank into four equal parts each with the height H/4 , where H = 1\, \mathrm{m} is used. These three inlets are defined so that they lie inside a cell for any mesh size considered. A fixed quantity of the is introduced through inlet z_{\mathrm{F}, 1} , a fixed quantity of the secondary disperse phase through inlet z_{\mathrm{F}, 2} and some wash water through inlet z_{\mathrm{F}, 3} .

    Tables 1 and 2 show the estimated errors and convergence orders for the three scenarios studied. In the calculations of the alternative approximate error \smash{\hat{e}^\mathrm{tot}_{N_k}(T)} and convergence order \smash{\hat{\Upsilon}_{k} (T)} in Table 1, we use \hat{k} = 6 .

    Table 1.  Simulation of a smooth solution (Section 6.3). Total estimated relative L^1 -error \smash{e^\mathrm{tot}_{N_k}(T)} , alternative relative L^1 -error \smash{\hat{e}^\mathrm{tot}_{N_k}(T)} , estimated convergence order \Upsilon_k (T) , and its alternative counterpart \smash{\hat{\Upsilon}_{k} (T)} , calculated with N_{\mathrm{ref}} = 128\, 000 and T = 9 \, \mathrm{s} .
    N_k e^\mathrm{tot}_{N_k}(T) \Upsilon_k (T) \hat{e}^\mathrm{tot}_{N_k}(T) \hat{\Upsilon}_{k} (T)
    500 3.7212\times 10^{-2} —— 1.3041\times 10^{-3} 0.9513
    1000 1.8985\times 10^{-2} 0.9709 6.7443\times 10^{-4} 0.9657
    2000 9.5710\times 10^{-3} 0.9881 3.4533\times 10^{-4} 0.9781
    4000 4.7582\times 10^{-3} 1.0083 1.7531\times 10^{-4} 0.9870
    8000 2.3174\times 10^{-3} 1.0379 8.8448\times 10^{-5} 0.9927
    16000 1.0867\times 10^{-3} 1.0926 4.4447\times 10^{-5} ——

     | Show Table
    DownLoad: CSV
    Table 2.  Approximate total relative L^1 -errors \smash{e^\mathrm{tot}_{N_k}(T)} and convergence orders \smash{\Upsilon_k(T)} calculated between consecutive values of N_k , with N_{\mathrm{ref}} = 12\, 800 (a) for Application 1 (counter-current flow) at simulated time T = 350 \, \mathrm{s} , (b) for Application 2 (co-current flow) at simulated time T = 500 \, \mathrm{s} .

     | Show Table
    DownLoad: CSV

    We consider a vessel with a constant cross-sectional area of A(z) = 83.65 \, \mathrm{cm}^2 , and we set all inlet and outlet volumetric flows to zero, i.e, Q_{\mathrm{F, 1}} = Q_{\mathrm{F, 2}} = Q_{\mathrm{F, 3}} = Q_{\mathrm{U}} = Q_{\mathrm{E}} = 0 \, \mathrm{cm^3/s} . (Under these assumptions, the scheme reduces to Scheme 3 for Model 3.) For the velocity functions W and V , given by Eq (2.10) and Eq (2.11), respectively, we use the parameters n_{\mathrm{p}} = 2.2 , v_{\mathrm{term, p}} = 1.5\, \mathrm{cm/s} , n_{\mathrm{s}} = 2.2 and v_{\mathrm{term, s}} = 1.5\, \mathrm{cm/s} , and consider \sigma = -1 (counter-current flow). The initial datum is a sinusoidal function for both phases with support in the interval (z_\mathrm{U}, z_\mathrm{E}) .

    We simulate a short time, until t = 9\, \mathrm{s} , before the first discontinuity appears; see the first row in Figure 4 where N = 1000 is used. In the second row, we compare two approximate solutions obtained with a coarse mesh with N = 500 and a finer one with N = 8000 . Table 1 shows the estimated errors and convergence orders. Both \Upsilon_k (T) and \hat{\Upsilon}_{N_k} (T) assume values close to one as N_k increases, as expected, confirming that the scheme is first-order accurate for smooth solution.

    Figure 4.  Simulation of a smooth solution (Section 6.3). First row: Time evolution of the volume fractions of the primary disperse phase \phi (left) and the secondary disperse phase \psi (right) from t = 0 \, \mathrm{s} to t = 9 \, \mathrm{s} . Second row: Approximate solutions at time t = 9 \, \mathrm{s} computed with N = 500 (left) and N = 8000 (right).

    The remainder of examples refer to Model 1 discretized by Scheme 1. We first illustrate that the "crossing condition" is satisfied as mentioned in Section 4.1. For this we use the constant A \equiv 83.65\, \mathrm{cm}^2 and simulate a tank that initially contains only water, i.e., \phi(z, 0) = \psi(z, 0) = 0 for all z . At t = 0 we start pumping aggregates, solids, fluid and wash water with \phi_{\mathrm{F}, 1} = 0.9 , \psi_{\mathrm{F}, 1} = 0 , \phi_{\mathrm{F}, 2} = 0.2 , \psi_{\mathrm{F}, 2} = 0.4 , \phi_{\mathrm{F}, 3} = 0.1 and \psi_{\mathrm{F}, 3} = 0 . We choose the volumetric flows (Q_{\mathrm{U}}, Q_{\mathrm{F, 1}}, Q_{\mathrm{F, 2}}, Q_{\mathrm{F, 3}}) = (15, 20, 25, 15)\, \mathrm{cm^3/s} , so that the volumetric flows in the tank are positive in all zones but not in zone 1. Three inlets z_{\mathrm{F}, 1} , z_{\mathrm{F}, 2} and z_{\mathrm{F}, 3} divide the tank into four zones of equal height.

    Figure 5 shows the graphs of the flux functions on both sides of each discontinuity in z , for this case three inlets (z_{\mathrm{F}, 1}, z_{\mathrm{F}, 2}, z_{\mathrm{F}, 3}) and two outlets (z_{\mathrm{U}}, z_{\mathrm{E}}) . We see that the fluxes J(\phi, z^{\pm}) (defined in Eq (4.9)) intersect when \phi_{\mathrm{F}, 1} = 0.9 , \phi_{\mathrm{F}, 2} = 0.2 and \phi_{\mathrm{F}, 3} = 0.1 , and do not intersect in (0, 1) at neither z_{\mathrm{U}} nor z_{\mathrm{E}} . Figure 6 shows the simulation results during 200 \, \mathrm{s} .

    Figure 5.  Illustration of the crossing condition (Section 6.4). The crossing condition is satisfied at each of the five spatial discontinuities.
    Figure 6.  Simulation of the example in Section 6.4 during T = 200 \, \mathrm{s} .

    We consider now a complete tank with \sigma = -1 ; hence, the primary disperse phase will move upwards and the secondary disperse phase downwards with respect to the volume average velocity q of the mixture. A straightforward interpretation of this scenario is the flotation process used in the mineral industry to recover valuable minerals from crushed ore; see the model in [8,9]. In that model, the primary disperse phase consists of aggregates, which are air bubbles fully loaded with hydrophobic minerals, and the secondary disperse phase is the tailings, consisting of hydrophilic particles suspended in the fluid that do not attach to air bubbles. We consider three inlets z_{\mathrm{F}, 1} , z_{\mathrm{F}, 2} and z_{\mathrm{F}, 3} , dividing the tank into four regions with equal height. At z_{\mathrm{F}, 1} , only gas is fed, at z_{\mathrm{F}, 3} only wash water, while at z_{\mathrm{F}, 2} a slurry of solids and water is fed into the column.

    The cross-sectional area is discontinuous (cf. Figure 2) due to a centered pipe from the top down to z_{\mathrm{F}, 2} that introduces material into the tank. It is given by

    \begin{align*} A(z) = \begin{cases} 72.25\, \mathrm{cm^2} & \text{for }z \geq z_{\mathrm{F, 2}} , \\ 83.65\, \mathrm{cm^2} & \text{for }z < z_{\mathrm{F, 2}} . \end{cases} \end{align*}

    These values correspond to the reflux flotation cell studied in [21].

    We consider that the column is initially filled only with fluid, hence \phi(z, 0) = \psi(z, 0) = 0 for all z , when we start pumping aggregates and solids with concentrations \phi_{\mathrm{F}, 1} = 1.0 , \psi_{\mathrm{F}, 1} = 0 , \phi_{\mathrm{F}, 2} = 0 , \psi_{\mathrm{F}, 2} = 0.4 , \phi_{\mathrm{F}, 3} = 0 and \psi_{\mathrm{F}, 3} = 0 , along with fluid and/or wash water. We choose (Q_{\mathrm{U}}, Q_{\mathrm{F, 1}}, Q_{\mathrm{F, 2}}, Q_{\mathrm{F, 3}}) = (5, 15, 25, 10)\, \mathrm{cm^3/s} , so that the mixture flows in zones 2 and 3 are positive, i.e., directed upwards: Q_{\mathrm{F, 1}}-Q_{\mathrm{U}} = 10\, \mathrm{cm^3/s} in zone 2 and Q_{\mathrm{F, 2}}+Q_{\mathrm{F, 1}}-Q_{\mathrm{U}} = 35\, \mathrm{cm^3/s} in zone 3.

    Figure 7 shows the time evolution of the volume fractions of \phi and \psi . It can be seen that the aggregates rise fast to the top, while the solids are travelling both up and down the vessel, leaving through the effluent and the underflow.

    Figure 7.  Application 1: Counter-current flows. Time evolution of the volume fraction profiles of the primary disperse phase \phi (left) and secondary disperse phase \psi (right) from time t = 0 \, \mathrm{s} to t = 1800 \, \mathrm{s} seen from two different angles (first and second rows).

    At time t = 350\, \mathrm{s} , we change the volumetric flow from Q_{\mathrm{F}, 2} = 25\, \mathrm{cm^3/s} to Q_{\mathrm{F}, 2} = 7\, \mathrm{cm^3/s} . After this change, the solids settle and we obtain a steady state. We mention that this is not a desired steady state in the mining industry (the capacity of the device is not fully used); see [9] for more examples. Table 2 (a) shows the estimated errors and convergence orders for this simulation. As in the smooth example in Section 6.3, the convergence orders tend to one as N_k increases.

    For the last example, we consider \sigma = 1 , i.e., both the primary and secondary disperse phases have a density smaller than that of the fluid and therefore move upwards relative to the mixture. This scenario could be a flotation process with two buoyant phases differing in density and possibly also in size. We consider here the same flotation column as in Application 1 and choose n_{\mathrm{p}} = 3.2 , v_{\mathrm{term, p}} = 2.5\, \mathrm{cm/s} , n_{\mathrm{s}} = 2.5 , and v_{\mathrm{term, s}} = 1.5\, \mathrm{cm/s} so that we have two buoyant phases with different (upwards-directed) velocities relative to the mixture. As in the previous example, only the primary disperse phase is fed into the tank at z_{\mathrm{F}, 1} and only the secondary at z_{\mathrm{F}, 2} . The column is initially filled with only fluid at time t = 0\, \mathrm{s} , hence \phi(z, 0) = \psi(z, 0) = 0 for all z , when we start pumping both phases with the following volume fractions: \phi_{\mathrm{F}, 1} = 1.0 , \psi_{\mathrm{F}, 1} = 0.0 , \phi_{\mathrm{F}, 2} = 0.0 , \psi_{\mathrm{F}, 2} = 0.6 , \phi_{\mathrm{F}, 3} = 0 and \psi_{\mathrm{F}, 3} = 0 . We choose the volumetric flows (Q_{\mathrm{U}}, Q_{\mathrm{F, 1}}, Q_{\mathrm{F, 2}}, Q_{\mathrm{F, 3}}) = (15, 30, 20, 10)\, \mathrm{cm^3/s} , so that the volumetric flows in the tank are positive in all zones with the exception of zone 1.

    Figure 8 shows the time evolution of the volume fractions of both phases. It can be seen that, for times t < 350\, \mathrm{s} , the primary disperse phase leaves the tank through both the underflow and effluent outlets, while the secondary disperse phase quickly rises to the top part of the tank and leaves it just through the effluent outlet.

    Figure 8.  Application 2: Co-current flow. Time evolution of the volume fraction profiles of primary disperse phase \phi (left) and secondary disperse phase \psi (right) from time t = 0 \, \mathrm{s} to t = 1500 \, \mathrm{s} seen from two different angles (first and second rows).

    At t = 350 \, \mathrm{s} , we change the volumetric flow of the inlet z_{\mathrm{F}, 1} from Q_{\mathrm{F}, 1} = 30 to Q_{\mathrm{F}, 1} = 20\, \mathrm{cm^3/s} , maintaining the other volumetric flows. As a consequence we can see that the primary disperse phase \phi rises and leaves zone 1, exiting the tank only through the effluent while the secondary disperse phase maintains the same behaviour as before and is present only above the feed level z_{\mathrm{F, 2}} . Table 2 (b) shows the estimated errors and convergence orders for this simulation, which have the same behaviour as the ones in the numerical examples in Sections 6.3 and 6.4.

    The present study outlines a numerical method for a triangular system of two PDEs, whose flux functions have several spatial discontinuities due to in- and outflows of a one-dimensional tank with possibly varying cross-sectional area. The triangular structure is utilized in the following way in the numerical scheme. The numerical update formula corresponding to the first scalar equation contains, for the nonlinear term, a numerical flux where the the volume fraction in the left cell is multiplied with the velocity computed in the right cell; see [6]. The update formula for the second equation uses the Engquist-Osher numerical flux for the term modelling the nonlinear relative flux of the secondary disperse phase, chosen in a particular way since this flux also depends on the primary disperse phase volume fraction. The other terms of the second update formula are also chosen in such a way that the entire scheme is proved to be monotone under the CFL condition Eq (3.2). We prove that the numerically obtained volume fractions satisfy the invariant-region property that they stay between zero and one, as is physically expected.

    The numerical scheme is applied to simulate the hydrodynamic movement of simultaneously rising aggregates (air bubbles with attached hydrophobic particles) and settling hydrophilic particles in the fluid under in- and outflows of a flotation column. As a demonstration of the capabilities of the numerical method, three different settings are simulated. The convergence order of the numerical method is estimated. As expected, in regions where the solution is smooth, the order is one. The first-order scheme proposed in this paper could be improved to achieve second-order accuracy, for instance, by techniques of variable extrapolation [6,14].

    In [9], the authors proposed a staggered scheme to compute numerical solutions for a flotation column, following the approach of Karlsen et al. [32,33]. Although the staggered scheme worked for a single inlet for a mixture of aggregates and solids, we have, in the case of several feed inlets, found it difficult to get consistent numerical solutions with respect to different mesh sizes.

    We are currently [12] extending the model and numerical scheme to the explicit description of drainage of liquid from the foam forming at the top of a flotation column. This phenomenon gives rise to a model similar to Eq (1.1) but with an additional degenerating diffusion term. The numerical solution of the resulting system of non-linear convection-diffusion equations calls for semi-implicit discretizations to alleviate the severe restrictions in the CFL condition due to the diffusion term.

    R.B. is supported by ANID (Chile) through Anillo project ANID/PIA/ACT210030; Centro de Modelamiento Matemático (BASAL projects ACE210010 and FB210005); CRHIAM, project ANID/Fondap/15130015; and Fondecyt project 1210610. S.D. is supported by the Swedish Research Council (Vetenskapsrådet, 2019-04601). M.C.M. is supported by grant PID2020-117211GB-I00 funded by MCIN/AEI/10.13039/501100011033 and project CIAICO/2021/227. Y.V. is supported by IFARHU-SENACYT (Panama) and project UCO1866.

    The authors declare there is no conflict of interest.

    Proof of Lemma 5.5. Following [18], we work with smooth entropies instead of \eta_0 , so we denote by \eta_{\Delta z} a smooth convex approximation to \eta_0 , so that \eta_{\Delta z} (0) = 0 and |\eta_{\Delta z} | \leq 1 , and \| \eta_{\Delta z} - \eta_0 \|_{L^{\infty}} \leq \Delta z . Moreover, if Q_{\Delta z} is the entropy flux associated with \eta_{\Delta z} , then there also holds \| Q_{\Delta z} - Q_0 \|_{L^{\infty}} \to 0 as \Delta z \to 0 . Then we split \mu^{\Delta z} as \smash{\mu^{\Delta z} = \mu_1^{\Delta z} + \mu_2^{\Delta z}} , where we define

    \begin{align*} \mu_1^{\Delta z} : = \partial_t \bigl( \eta_0 ( \psi^{\Delta z} ) - \eta_{\Delta z} ( \psi^{\Delta z} ) \bigr) , \quad \mu_2^{\Delta z} : = \partial_t \eta_{\Delta z} ( \psi^{\Delta z} ) + \partial_z Q_0 ( \phi^{\Delta z}, \psi^{\Delta z} ). \end{align*}

    If \zeta \in C_0^1 (\Pi_T) denotes a test function with compact support, then as in [18], one verifies that

    \begin{align*} \bigl| \langle \mu_1^{\Delta z}, \zeta \rangle \bigr| \leq \iint_{\Pi_T} \bigl| \eta_{\Delta z} ( \psi^{\Delta z} ) - \eta_0 ( \psi^{\Delta z} ) \bigr| |\zeta_t | \, \mathrm{d} z \, \mathrm{d} t \leq C_{21} \| \zeta_t \|_{L^2 ( \Pi_T)} \| \eta_{\Delta z} - \eta_0 \|_{L^{\infty}} \to 0 \quad \text{as } \Delta z \to 0 , \end{align*}

    hence \smash{ \{ \mu_1^{\Delta z} \}_{\Delta z > 0} } is compact in H_{\mathrm{loc}}^{-1} (\Pi_T) . By an integration by parts we get

    \begin{align*} \langle \mu_2^{\Delta z}, \zeta \rangle & = - \iint_{\Pi_T} \bigl( \eta_{\Delta z} ( \psi^{\Delta z} ) \zeta_t + Q_0 ( \phi^{\Delta z}, \psi^{\Delta z} ) \zeta _z \bigr) \, \mathrm{d} z \, \mathrm{d} t \\ & = - \sum\limits_{n = 1}^{N_T-1} \int_{\mathbb{R}} \int_{t_n}^{t_{n+1}} \eta_{\Delta z} ( \psi^{\Delta z} ) \zeta_t \, \mathrm{d} t \, \mathrm{d} z - \sum\limits_{n = 1}^{N_T-1} \int_{t_n}^{t_{n+1}} \int_{\mathbb{R}} Q_{0} ( \phi^{\Delta z}, \psi^{\Delta z} ) \zeta_z \, \mathrm{d} z \, \mathrm{d} t \\ & = - \sum\limits_{n = 0}^{N_T-1} \int_{\mathbb{R}} \eta_{\Delta z} \bigl( \psi^{\Delta z} (z, t_n) \bigr) \bigl( \zeta ( z, t_{n+1} ) - \zeta ( z, t_n) \bigr) \, \mathrm{d} z \\ & \quad - \sum\limits_{n = 1}^{N_T-1} \sum\limits_{ i \in \mathbb{Z}} \int_{t_n}^{t_{n+1}} Q_{0} \bigl( \phi_{i-1/2}^n, \psi_{i-1/2}^n ) \bigl( \zeta (z_i, t) - \zeta( z_{i-1}, t) \bigr) \, \mathrm{d} z \, \mathrm{d} t, \end{align*}

    so we may finally write

    \begin{align} \begin{split} \langle \mu_2^{\Delta z}, \zeta \rangle & = \sum\limits_{n = 0}^{N_T-2} \sum\limits_{i \in \mathbb{Z}} \bigl( \eta_{\Delta z} ( \psi_{i-1/2}^{n+1} ) - \eta_{\Delta z} (\psi_{i-1/2}^{n} ) \bigr) \int_{I_{i-1/2}} \zeta (z, t_{n+1}) \, \mathrm{d}z \\ & \quad + \sum\limits_{n = 1}^{N_T-1} \sum\limits_{ i \in \mathbb{Z}} \bigl( \Delta_+ Q_0 ( \phi_{i-1/2}^n , \psi_{i-1/2}^n) \bigr)\int_{t_n}^{t_{n+1}} \zeta ( z_i, t) \, \mathrm{d} t. \end{split} \end{align} (A.1)

    We define the cell average

    \begin{align} \zeta_{i-1/2}^n : = \frac{1}{\Delta z \Delta t} \iint_{I_{j-1/2}^n} \zeta (z, t) \, \mathrm{d} z \, \mathrm{d}{t}. \end{align} (A.2)

    Replacing the integral in the first term of the right-hand side of Eq (A.1) by \Delta z \zeta_{i-1/2}^n produces the following error, where we follow the derivation of Eq (3.27) in [18]:

    \begin{align*} & \left| \sum\limits_{\mathcal{I}_2} \bigl( \eta_{\Delta z} ( \psi_{i-1/2}^{n+1} ) - \eta_{\Delta z} ( \psi_{i-1/2}^{n} ) \bigr) \biggl( \int_{I_{i-1/2}} \zeta (z, t_{n+1}) \, \mathrm{d} z - \Delta z \zeta_{i-1/2}^n \biggr) \right| \\ & \leq \| \eta_{\Delta z}' \|_{L^{\infty}} \sum\limits_{\mathcal{I}_2} \bigl| \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^{n} \bigr| \frac{1}{\Delta t} \iint_{I_{j-1/2}^n} \bigl| \zeta ( z, t_{n+1} ) - \zeta ( z, t) \bigr| \, \mathrm{d} z \, \mathrm{d} t \\ & \leq \sum\limits_{n, i} \bigl| \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^{n} \bigr| \frac{1}{\Delta t} \iint_{I_{j-1/2}^n} \int_t^{t_{n+1}} \bigl| \zeta_t (z, s) \bigr| \, \mathrm{d} s \, \mathrm{d} z \, \mathrm{d} t \\ & \leq \sum\limits_{\mathcal{I}_2} \bigl| \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^{n} \bigr| \frac{1}{\Delta t} \iint_{I_{j-1/2}} \int_{t_n}^{t_{n+1}} (t_{n+1}-t)^{1/2} \biggl( \int_{t_n}^{t_{n+1}} \bigl| \zeta_t (z, s) \bigr|^2 \, \mathrm{d} s\biggr)^{1/2} \mathrm{d} z \, \mathrm{d} t \\ & \leq \frac{2}{3} \sum\limits_{\mathcal{I}_2} \bigl| \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^{n} \bigr| \Delta t^{1/2} \int_{I_{i-1/2}} \biggl( \int_{t_n}^{t_{n+1}} \bigl| \zeta_t (z, s) \bigr|^2 \, \mathrm{d} s\biggr)^{1/2} \mathrm{d} z \\ & \leq \frac{2}{3} \sum\limits_{\mathcal{I}_2} \bigl| \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^{n} \bigr| \Delta z^{1/2} \Delta t^{1/2} \biggl( \iint_{I_{i-1/2}^n} \bigl( \zeta_t (z, s) \bigr)^2 \, \mathrm{d} z \, \mathrm{d} t \biggr)^{1/2} \\ & \leq \frac{2}{3} \Biggl( \Delta z \Delta t \sum\limits_{n, i} \bigl( \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^{n} \bigr)^2 \Biggr)^{1/2} \Biggl( \sum\limits_{\mathcal{I}_2} \iint_{I_{i-1/2}^n} \bigl( \zeta_t (z, s) \bigr)^2 \, \mathrm{d} z \, \mathrm{d} t \Biggr)^{1/2} \\ & \leq \frac{2}{3} (C_7 \Delta z)^{1/2} \| \zeta \|_{H^1 (\Pi_T)} \end{align*}

    (see Eq (5.17)). By similar arguments we obtain the bound

    \begin{align*} \Biggl| \sum\limits_{\mathcal{I}_1} \bigl( \Delta_+ Q_0 ( \phi_{i-1/2}^n , \psi_{i-1/2}^n) \bigr) \left(\int_{t_n}^{t_{n+1}} \zeta(z_{i+1/2}, t ) \, \mathrm{d} t - \Delta t \zeta_{i-1/2}^n \right) \Biggr| \leq \tilde{C}_T \Delta z^{1/2} \| \zeta \|_{H^1 ( \Pi_T)}. \end{align*}

    Consequently, and further following [18], we have shown that

    \begin{align*} \langle \mu_2^{\Delta z}, \zeta \rangle & = \Delta z \Delta t \sum\limits_{\mathcal{I}_1} \biggl\{ \frac{ \eta_{\Delta z} ( \psi_{i-1/2}^{n+1} ) - \eta_{\Delta z} ( \psi_{i-1/2}^{n} )}{\Delta t} + \frac{\Delta_+ Q_0 ( \phi_{i-1/2}^n, \psi_{i-1/2}^n )}{\Delta z} \biggr\} \zeta_{i-1/2}^n \\ & \quad + \text{terms which are compact in } H_{\mathrm{loc}}^{-1} ( \Pi_T ) . \end{align*}

    We now utilize the "scheme for \eta ", Eq (5.14), to rewrite the term in curled brackets as \smash{\{ \dots \} = \mathcal{A}_{i-1/2}^n + \mathcal{B}_{i-1/2}^n + \mathcal{C}_{i-1/2}^n + \mathcal{D}_{i-1/2}^n} , where we define

    \begin{align} \mathcal{A}_{i-1/2}^n & : = - \frac{1}{2\Delta t} \eta_{\Delta z}'' \bigl( \xi_{i-1/2}^{n+1/2} \bigr) \bigl( \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^n \bigr)^2 - \frac{1}{\Delta z} \Theta_{i-1/2}^n , \\ \mathcal{B}_{i-1/2}^n & : = - \frac{1}{\Delta z} \bigl( \eta' \bigl( \psi_{i-1/2}^{n+1} \bigr) - \eta' \bigl( \psi_{i-1/2}^{n} \bigr) \bigr) \Delta_- \mathcal{F} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i), \\ \mathcal{C}_{i-1/2}^n & : = - \frac{1}{\Delta z} \Delta_+^{\phi} \bigl( \eta' \bigl( \psi_{i-1/2}^n \bigr) \mathcal{F} ( \boldsymbol{\phi}_{i-1}^n, \boldsymbol{\psi}_{i-1}^n) - \mathcal{Q} ( \boldsymbol{\phi}_{i-1}^n, \boldsymbol{\psi}_{i-1}^n) \bigr) , \\ \mathcal{D}_{i-1/2}^n & : = - \frac{1}{\Delta z} \bigl( \Delta_- \mathcal{Q} ( \boldsymbol{\phi}_i^n, \boldsymbol{\psi}_i^n ) - \Delta_+ Q_0 \bigl( \phi_{i-1/2}^n, \psi_{i-1/2}^n \bigr) \bigr). \end{align} (A.3)

    Thus, \langle \mu_2^{\Delta z}, \zeta \rangle = \langle \mathcal{A}, \zeta \rangle + \langle \mathcal{B}, \zeta \rangle + \langle \mathcal{C}, \zeta \rangle + \langle \mathcal{D}, \zeta \rangle + \text{compact terms, } where

    \begin{align*} \langle \mathcal{A}, \zeta \rangle = \Delta z \Delta t \sum\limits_{n, i} \mathcal{A}_{i-1/2}^n \zeta_{i-1/2}^n \end{align*}

    and \langle \mathcal{B}, \zeta \rangle , \langle \mathcal{C}, \zeta \rangle , and \langle \mathcal{D}, \zeta \rangle are defined analogously. In view of Lemma 5.2, we get

    \begin{align*} \bigl| \langle \mathcal{A}, \zeta \rangle \bigr| & \leq \| \zeta \|_{L^{\infty} (\Pi_T) } \biggl( \frac{\Delta z}{2} \sum\limits_{\mathcal{I}_1} \eta_{\Delta z}'' \bigl( \xi_{i-1/2}^{n+1/2} \bigr) \bigl( \psi_{i-1/2}^{n+1} - \psi_{i-1/2}^n \bigr)^2 + \Delta t \sum\limits_{\mathcal{I}_1} \Theta_{i-1/2}^n \biggr) \leq C_T \| \zeta \|_{L^{\infty} ( \Pi_T ) }, \end{align*}

    and therefore \mathcal{A} \in \mathcal{M}_{\mathrm{loc}} (\Pi_T) . Appealing to the divergence bound of the numerical flux Eq (5.15) and taking into account the BV bound on \phi^{\Delta z} it also follows that | \langle \mathcal{B} + \mathcal{C}, \zeta \rangle | \leq C_T \| \zeta \|_{L^{\infty} (\Pi_T)} , and therefore \mathcal{B} + \mathcal{C} \in \mathcal{M}_{\mathrm{loc}} (\Pi_T) .

    Finally, to deal with \langle \mathcal{D}, \zeta \rangle we consider first \varepsilon > 0 and let Q_{\varepsilon} , \smash{\mathcal{Q}_{\varepsilon}^{\pm}} and \mathcal{Q}_{\varepsilon} denote the entropy and numerical entropy fluxes calculated from Eq (5.1) and Eq (5.7), respectively, where \eta = \eta_{\varepsilon} . Since \mathcal{Q}_{\varepsilon} is consistent with Q_{\varepsilon} ,

    \begin{align*} & \mathcal{Q}_{\varepsilon} ( \phi, \phi, \psi_1 , \psi_2) - Q_{\varepsilon} (\phi, \psi_2) \\ & = \mathcal{Q}_{\varepsilon} ( \phi, \phi, \psi_1 , \psi_2) - \mathcal{Q}_{\varepsilon} ( \phi, \phi, \psi_2 , \psi_2) = \int_{\psi_2}^{\psi_1} \eta_{\varepsilon}' (s) \check{f} ( \phi, \phi, s) \, \mathrm{d} s \\ & = \eta_{\varepsilon}' (\psi_1) \bigl( \check{\mathcal{F}} ( \phi, \phi, \psi_1) - \check{\mathcal{F}} ( \phi, \phi, \psi_2) \bigr) - \int_{\psi_2}^{\psi_1} \eta_{\varepsilon}'' (s) \bigl( \check{\mathcal{F}} ( \phi, \phi, s) - \check{\mathcal{F}} ( \phi, \phi, \psi_1) \bigr) \, \mathrm{d} s \end{align*}

    (cf. Eq (5.9)). By using the monotonicity of \check{\mathcal{F}} with respect to its \psi -argument we get

    \begin{align*} \bigl| \mathcal{Q}_{\varepsilon} ( \phi, \phi, \psi_1 , \psi_2) - Q_{\varepsilon} (\phi, \psi_2) \bigr| \leq 3 \| \eta_{\varepsilon}' \|_{L^{\infty}} \bigl| \check{\mathcal{F}} ( \phi, \phi, \psi_2) - \check{\mathcal{F}} ( \phi, \phi, \psi_1) \bigr| \leq 3 \bigl| \check{\mathcal{F}} ( \phi, \phi, \psi_2) - \check{\mathcal{F}} ( \phi, \phi, \psi_1) \bigr|, \end{align*}

    so in the limit \varepsilon \to 0 ,

    \begin{align} \bigl| \mathcal{Q} ( \phi, \phi, \psi_1 , \psi_2) - Q_0 (\phi, \psi_2) \bigr| & \leq 3 \bigl| \check{\mathcal{F}} ( \phi, \phi, \psi_2) - \check{\mathcal{F}} ( \phi, \phi, \psi_1) \bigr|. \end{align} (A.4)

    Noticing that

    \begin{align*} \Delta_- \mathcal{Q} ( \boldsymbol{\phi}^n_i, \boldsymbol{\psi}^n_i ) - \Delta_+ Q_0 \bigl( \phi_{i-1/2}^n, \psi_{i-1/2}^n \bigr) & = \Delta_- \bigl( \mathcal{Q} \bigl( \boldsymbol{\phi}_{i}^n, \psi_{i-1/2}^n, \psi_{i+1/2}^n \bigr) - Q_0 ( \phi_{i+1/2}^n, \psi_{i+1/2}^n \bigr) \bigr) \\ & = \Delta_- \bigl( \mathcal{Q} \bigl( \phi_{i+1/2}^n, \phi_{i+1/2}^n, \boldsymbol{\psi}_i^n \bigr) - Q_0 ( \phi_{i+1/2}^n, \psi_{i+1/2}^n \bigr) \bigr) \\ & \quad + \Delta_- \bigl( \mathcal{Q} \bigl( \boldsymbol{\phi}_{i}^n, \boldsymbol{\psi}_{i}^n \bigr) - \mathcal{Q} \bigl( \phi_{i+1/2}^n, \phi_{i+1/2}^n, \boldsymbol{\psi}_{i}^n \bigr) \bigr) \end{align*}

    we get from Eq (A.3)

    \begin{array}{l} \Biggl| \Delta z \Delta t \sum\limits_{\mathcal{I}_1} \mathcal{D}_{i-1/2}^n \zeta_{i-1/2}^n \Biggr| \leq \Biggl| \Delta z \Delta t \sum\limits_{\mathcal{I}_1} \frac{1}{\Delta z} \Delta_- \bigl( \mathcal{Q} \bigl( \phi_{i+1/2}^n, \phi_{i+1/2}^n, \boldsymbol{\psi}_{i}^n \bigr) - Q_0 ( \phi_{i+1/2}^n, \psi_{i+1/2}^n \bigr)\bigr) \zeta_{i-1/2}^n \Biggr| \\ + \Biggl| \Delta z \Delta t \sum\limits_{\mathcal{I}_1} \frac{1}{\Delta z} \Delta_- \bigl( \mathcal{Q} \bigl( \boldsymbol{\phi}_{i}^n , \boldsymbol{\psi}_{i}^n \bigr) - \mathcal{Q} \bigl( \phi_{i+1/2}^n, \phi_{i+1/2}^n, \boldsymbol{\psi}_{i}^n \bigr) \bigr)\zeta_{i-1/2}^n \Biggr| = : |\mathcal{S}_1| + |\mathcal{S}_2|. \end{array}

    By a summation by parts and applying Eq (A.4) we get

    \begin{align} |\mathcal{S}_1| & = \Biggl|\Delta z \Delta t \sum\limits_{\mathcal{I}_1} \bigl( \mathcal{Q} \bigl( \phi_{i+1/2}^n, \phi_{i+1/2}^n, \boldsymbol{\psi}_{i}^n \bigr) - Q_0 ( \phi_{i+1/2}^n, \psi_{i+1/2}^n \bigr)\bigr) \frac{\Delta_+ \zeta_{i-1/2}^n}{\Delta z} \Biggr| \\ & \leq 3 \Delta z \Delta t \sum\limits_{\mathcal{I}_1} \bigl| \Delta_+^{(3)} \check{\mathcal{F}} \bigl( \phi_{i+1/2}^n, \phi_{i+1/2}^n , \psi_{i-1/2}^n \bigr)\bigr| \frac{|\Delta_+ \zeta_{i-1/2}^n|}{\Delta z}. \end{align} (A.5)

    We now write

    \begin{align*} \Delta_+^{(3)} \check{\mathcal{F}} \bigl( \phi_{i+1/2}^n, \phi_{i+1/2}^n, \psi_{i-1/2}^n \bigr) = \Delta_+^{(3)} \check{\mathcal{F}} \bigl( \boldsymbol{\phi}_{i}^n, \psi_{i-1/2}^n \bigr) + \mathcal{Y}_{i+1/2}^n - \mathcal{Y}_{i-1/2}^n, \end{align*}

    where

    \begin{align*} \mathcal{Y}_{i \pm 1/2}^n : = \check{\mathcal{F}} \bigl( \phi_{i+1/2}^n, \phi_{i+1/2}^n, \psi_{i \pm 1/2}^n \bigr) -\check{\mathcal{F}} \bigl( \boldsymbol{\phi}_{i}^n, \psi_{i \pm 1/2}^n \bigr). \end{align*}

    From Eq (5.34), and considering \smash{\phi_{i-3/2}^n = \phi_{i-1/2}^n} in that bound, we deduce there exists a constant C such that \smash{| \mathcal{Y}_{i\pm 1/2}^n | \leq C | \Delta_+ \phi_{i-1/2}^n |} , therefore there exists (another) constant C such that

    \begin{align} \bigl| \mathcal{Y}_{i+1/2}^n - \mathcal{Y}_{i-1/2}^n \bigr| \leq C \bigl| \Delta_+ \phi_{i-1/2}^n \bigr|. \end{align} (A.6)

    Consequently, from Eq (A.5) we deduce that

    \begin{align*} | \mathcal{S}_1 | & \leq 3 \Delta z \Delta t \left( \sum\limits_{\mathcal{I}_1} \bigl| \Delta_+^{(3)} \check{\mathcal{F}} \bigl( \boldsymbol{\phi}_{i}^n, \psi_{i-1/2}^n \bigr)\bigr| \frac{|\Delta_+ \zeta_{i-1/2}^n|}{\Delta z} + \sum\limits_{\mathcal{I}_1} \bigl| \Delta_+ \phi_{i-1/2}^n \bigr| \frac{|\Delta_+ \zeta_{i-1/2}^n|}{\Delta z} \right) \\ & \leq 3 \Biggl( \Delta z \Delta t \Biggl( \sum\limits_{\mathcal{I}_1} \bigl( \Delta_+^{(3)} \check{\mathcal{F}} \bigl( \phi_{i+1/2}^n, \phi_{i+1/2}^n , \psi_{i-1/2}^n \bigr)\bigr)^2 + C \sum\limits_n \operatorname*{TV} ( \boldsymbol{\phi}^n ) \Biggr) \Biggr)^{1/2} \\ & \Biggl( \Delta z \Delta t \sum\limits_{\mathcal{I}_1} \biggl( \frac{|\Delta_+ \zeta_{i-1/2}^n|}{\Delta z} \biggr)^2 \Biggr)^{1/2} . \end{align*}

    From Eq (5.16) we infer that there exists a constant C_{T} such that

    \begin{align*} \Delta t \sum\limits_{\mathcal{I}_1} \bigl( \Delta_+^{(3)} \check{\mathcal{F}} \bigl( \phi_{i+1/2}^n, \phi_{i+1/2}^n , \psi_{i-1/2}^n \bigr)\bigr)^2 \leq C_{T}. \end{align*}

    Noticing that also

    \begin{align*} \Delta t \sum\limits_n \operatorname*{TV} ( \boldsymbol{\phi}^n ) \leq C_{t_N}, \end{align*}

    we conclude that there exists a constant C_{t_N} such that

    \begin{align*} \bigl| \langle \mathcal{S}_1 , \zeta \rangle \bigr| \leq C_{T_N} \Delta z^{1/2} \| \partial_z \zeta \|_{L^2 ( \Pi_T )}. \end{align*}

    Next, we deal with \mathcal{S}_2 . Applying again a summation by parts, we get

    \begin{align*} |\mathcal{S}_2 | & = \Biggl| \Delta z \Delta t \sum\limits_{n, i} \bigl( \mathcal{Q} \bigl( \boldsymbol{\phi}_{i}^n, \boldsymbol{\psi}_i^n \bigr) - \mathcal{Q} \bigl( \phi_{i+1/2}^n, \phi_{i+1/2}^n, \boldsymbol{\psi}_i^n \bigr) \bigr) \frac{\Delta_+ \zeta_{i-1/2}^n}{\Delta z} \Biggr|. \end{align*}

    The definition of \mathcal{Q} (see Eq (5.7)) yields

    \begin{align*} \bigl| \mathcal{Q} \bigl( \boldsymbol{\phi}_{i}^n, \boldsymbol{\psi}_i^n \bigr) - \mathcal{Q} \bigl( \phi_{i+1/2}^n, \phi_{i+1/2}^n, \boldsymbol{\psi}_i^n \bigr) \bigr| & \leq \bigl| \hat{\mathcal{Q}} ( \boldsymbol{\phi}_{i}^n , \psi_{i-1/2}^n ) - \hat{\mathcal{Q}} ( \phi_{i+ 1/2}^n, \phi_{i+1/2}^n, \psi_{i-1/2}^n ) \bigr| \\ & \quad + \bigl| \check{\mathcal{Q}} ( \boldsymbol{\phi}_{i}^n, \psi_{i+1/2}^n ) - \check{\mathcal{Q}} ( \phi_{i+ 1/2}^n, \phi_{i+1/2}^n, \psi_{i-1/2}^n ) \bigr|. \end{align*}

    By a computation similar to Eq (5.21) we get

    \begin{align*} \bigl| \hat{\mathcal{Q}} \bigl( \boldsymbol{\phi}_{i}^n , \psi_{i-1/2}^n \bigr) - \hat{\mathcal{Q}} \bigl( \phi_{i+ 1/2}^n, \phi_{i+1/2}^n, \psi_{i-1/2}^n \bigr) \bigr| \leq 3 \| \eta' \|_{\infty} \bigl| \mathcal{X}_{i-1/2}^n \bigr|, \end{align*}

    where

    \begin{align*} \mathcal{X}_{i-1/2}^n : = \hat{\mathcal{F}} \bigl( \boldsymbol{\phi}_{i}^n, \psi_{i-1/2}^n \bigr) - \hat{\mathcal{F}} \bigl( \phi_{i+1/2}^n, \phi_{i+1/2}^n, \psi_{i-1/2}^n \bigr). \end{align*}

    The discussion of \smash{\mathcal{X}_{i-1/2}^n} is similar to that of \smash{\mathcal{Y}_{i-1/2}^n} above, and appealing to Eq (5.34) we see that

    \begin{align*} \bigl| \mathcal{X}_{i-1/2}^n \bigr| \leq C \bigl| \Delta_+ \phi_{i-1/2}^n \bigr|. \end{align*}

    On the other hand, Eq (A.6) implies that

    \begin{align*} & \bigl| \check{\mathcal{Q}}\bigl( \boldsymbol{\phi}_{i}^n , \psi_{i-1/2}^n \bigr) - \check{\mathcal{Q}} \bigl( \phi_{i+ 1/2}^n, \phi_{i+1/2}^n, \psi_{i-1/2}^n \bigr) \bigr| \\ & \leq 3 \| \eta' \|_{\infty} \bigl| \check{\mathcal{F}} \bigl( \boldsymbol{\phi}_{i}^n, \psi_{i-1/2}^n \bigr) - \check{\mathcal{F}} \bigl( \phi_{i+1/2}^n, \phi_{i+1/2}^n, \psi_{i-1/2}^n \bigr) \bigr| \leq C \bigl| \Delta_+ \phi_{i-1/2}^n \bigr|. \end{align*}

    Thus

    \begin{align*} \bigl| \mathcal{Q} \bigl( \boldsymbol{\phi}_{i}^n, \boldsymbol{\psi}_i^n \bigr) - \mathcal{Q} \bigl( \phi_{i+1/2}^n, \phi_{i+1/2}^n, \boldsymbol{\psi}_i^n \bigr) \bigr| \leq C \bigl| \Delta_+ \phi_{i-1/2}^n \bigr|, \end{align*}

    and we deduce that \mathcal{S}_2 can be bounded in a similar way as \mathcal{S}_1 . In particular,

    \begin{align*} | \mathcal{S}_2 | & \leq 3 \Biggl( \Delta z \Delta t C \sum\limits_n \operatorname*{TV} ( \boldsymbol{\phi}^n ) \Biggr) \Biggr)^{1/2} \Biggl( \Delta z \Delta t \sum\limits_{n, i} \biggl( \frac{|\Delta_+ \zeta_{i-1/2}^n|}{\Delta z} \biggr)^2 \Biggr)^{1/2}, \end{align*}

    and we conclude that also

    \begin{align*} \bigl| \langle \mathcal{S}_2 , \zeta \rangle \bigr| \leq C_{T} \Delta z^{1/2} \| \partial_z \zeta \|_{L^2 ( \mathbb{R}^2 \times \mathbb{R}^+ )}, \end{align*}

    so \mathcal{D} is compact in H_{\mathrm{loc}}^{-1} (\Pi_T) . Thus, the sequence \{ \mu_2^{\Delta z} \}_{\Delta z > 0} , and therefore also the sequence \{ \mu^{\Delta z} \}_{\Delta z > 0} belong to a compact subset of H_{\mathrm{loc}}^{-1} (\Pi_T) .

    Proof of Theorem 5.1. We only need to verify that \psi is a weak solution of Eq (1.8b), that is, that Eq (1.10) holds. To this end, we choose a test function \zeta \in C_0^{\infty} (\Pi_T) , recall the definition Eq (A.2) of cell averages \smash{\zeta_{i-1/2}^n} , multiply the \psi -scheme Eq (5.3) by \Delta z \zeta_{i-1/2}^n , sum over i and n , and apply a summation by parts to obtain an identity \mathcal{J}_0 + \mathcal{J}_1 + \mathcal{J}_2 = 0 , where

    \begin{align*} & \mathcal{J}_0 : = \Delta z \sum\limits_{i \in \mathbb{Z}} \psi_{i-1/2}^0 \zeta_{i-1/2}^0, \quad \mathcal{J}_1 : = \Delta z \sum\limits_{i \in \mathbb{Z}} \sum\limits_{n = 1}^{N_T} \psi_{i-1/2}^n \bigl( \zeta_{i-1/2}^n- \zeta_{i-1/2}^{n-1} \bigr) , \\ & \mathcal{J}_2 : = \Delta z \Delta t \sum\limits_{n = 0}^{N_T-1} \sum\limits_{i \in \mathbb{Z}} \mathcal{F} ( \boldsymbol{\phi}_i^n, \boldsymbol{\psi}_i^n) {\frac{{ {{\Delta_+ \zeta_{i-1/2}^n}}}}{{ {{\Delta z}}}}}. \end{align*}

    By exactly following the estimates of terms I_0 and I_1 in the proof of [18,Lemma 3.5] and appealing to the bounded convergence theorem we may prove that

    \begin{align} \lim\limits_{\Delta z \to 0} \mathcal{J}_0 = \int_{\mathbb{R}} \psi_0( z) \zeta (z, 0) \, \mathrm{d} z, \quad \lim\limits_{\Delta z \to 0} \mathcal{J}_1 = \iint_{\Pi_T} \psi \partial_t \zeta \, \mathrm{d} z \, \mathrm{d} t. \end{align} (A.7)

    The treatment of \mathcal{J}_2 differs from that of the term I_2 in [18,Lemma 3.5] since here the numerical flux depends on four arguments (not three, as in [18]). We here get

    \begin{align*} \mathcal{J}_2 & = \iint_{\Pi_T} \tilde{F} ( \phi^{\Delta z}, \psi^{\Delta z} ) \partial_z \zeta \, \mathrm{d} z \, \mathrm{d} t + \mathcal{J}_{2, 1} + \mathcal{J}_{2, 2} + \mathcal{J}_{2, 3}, \end{align*}

    where we define

    \begin{align*} \mathcal{J}_{2, 1} & : = - \sum\limits_{\mathcal{I}_1} \tilde{F} \bigl( \phi_{i-1/2}^n, \psi_{i-1/2}^n \bigr) \iint_{I_{i-1/2}^n} \int_0^{\Delta z} \frac{\partial_z \zeta (z, t) - \partial_z \zeta ( z+ \xi, t)}{\Delta z} \, \mathrm{d} \xi \, \mathrm{d} z \, \mathrm{d} t, \\ \mathcal{J}_{2, 2} & : = - \Delta z \Delta t \sum\limits_{\mathcal{I}_1} \bigl( \tilde{F} \bigl( \phi_{i-1/2}^n, \psi_{i-1/2}^n \bigr) - \mathcal{F} \bigl( \boldsymbol{\phi}_i^n, \psi_{i-1/2}^n, \psi_{i-1/2}^n \bigr) \bigr) \frac{\Delta_+ \zeta_{i-1/2}^n}{\Delta z}, \\ \mathcal{J}_{2, 3} & : = - \Delta z \Delta t \sum\limits_{\mathcal{I}_1} \bigl( \mathcal{F} \bigl( \boldsymbol{\phi}_i^n, \psi_{i-1/2}^n, \psi_{i-1/2}^n \bigr) - \mathcal{F} \bigl( \boldsymbol{\phi}_i^n, \boldsymbol{\psi}_{i}^n \bigr)\bigr) \frac{\Delta_+ \zeta_{i-1/2}^n}{\Delta z} \\ & = \Delta z \Delta t \sum\limits_{\mathcal{I}_1} \Delta_+^{(3)} \check{\mathcal{F}} \bigl( \boldsymbol{\phi}_i^n, \psi_{i-1/2}^n \bigr) \frac{\Delta_+ \zeta_{i-1/2}^n}{\Delta z}. \end{align*}

    The term \mathcal{J}_{2, 1} can be estimated by choosing a constant M such that \zeta (z, t) = 0 for |z| > M and noting that

    \begin{align} |\mathcal{J}_{2, 1}| \leq \Delta z \| \partial_z^2 \zeta \|_{L^{\infty}} \| \tilde{F} (\phi^{\Delta z}, \psi^{\Delta z} ) \|_{L^1 ([-M, M] \times [0, T])} { \to 0\; {\rm{as}} \;\Delta z \to 0 .} \end{align} (A.8)

    Furthermore, in light of Eq (5.6) the difference arising in \mathcal{J}_{2, 2} can be written as

    \begin{array}{l} \tilde{F} \bigl( \phi_{i-1/2}^n, \psi_{i-1/2}^n \bigr) - \mathcal{F} \bigl( \boldsymbol{\phi}_i^n, \psi_{i-1/2}^n, \psi_{i-1/2}^n \bigr) \\ = \mathcal{F} \bigl( \phi_{i-1/2}^n, \phi_{i-1/2}^n, \psi_{i-1/2}^n, \psi_{i-1/2}^n \bigr) - \mathcal{F} \bigl( \phi_{i-1/2}^n, \phi_{i+1/2}^n, \psi_{i-1/2}^n, \psi_{i-1/2}^n \bigr) \\ = \hat{\mathcal{F}} \bigl( \phi_{i-1/2}^n, \phi_{i-1/2}^n, \psi_{i-1/2}^n \bigr) - \hat{\mathcal{F}} \bigl( \phi_{i-1/2}^n, \phi_{i+1/2}^n, \psi_{i-1/2}^n \bigr) \\ \quad + \check{\mathcal{F}} \bigl( \phi_{i-1/2}^n, \phi_{i-1/2}^n, \psi_{i-1/2}^n \bigr) - \check{\mathcal{F}} \bigl( \phi_{i-1/2}^n, \phi_{i+1/2}^n, \psi_{i-1/2}^n \bigr). \end{array}

    Utilizing the estimate Eq (5.34) with \phi_{i-3/2}^n = \phi_{i-1/2}^n yields that there exists a constant C_{12} such that

    \begin{align*} \bigl| \tilde{F} \bigl( \phi_{i-1/2}^n, \psi_{i-1/2}^n \bigr) - \mathcal{F} \bigl( \boldsymbol{\phi}_i^n, \psi_{i-1/2}^n, \psi_{i-1/2}^n \bigr) \bigr| \leq C_{12} \bigl| \Delta_+ \phi_{i-1/2}^n \bigr|, \end{align*}

    hence

    \begin{align} \bigl| \mathcal{J}_{2, 2} \bigr| \leq C_{12} \Delta z \Delta t \sum\limits_{n = 0}^{N_T -1} \operatorname*{TV} ( \boldsymbol{\phi}^n ) \frac{\Delta_+ \zeta_{i-1/2}^n}{\Delta z} \leq C \Delta z \| \partial_z \zeta \|_{L^{\infty}} \quad { \to 0 \;{\rm{as}} \;\Delta z \to 0 .} \end{align} (A.9)

    To estimate \mathcal{J}_{2, 3} , we utilize Eq (5.16). Then

    \begin{align} \begin{split} | \mathcal{J}_{2, 3} | & \leq \left( \Delta z \Delta t \sum\limits_{\mathcal{I}_1} \bigl( \Delta_+^{(3)} \check{\mathcal{F}} \bigl( \boldsymbol{\phi}_i^n, \psi_{i-1/2}^n \bigr) \bigr)^2 \right)^{1/2} \left( \Delta z \Delta t \sum\limits_{\mathcal{I}_1} \biggl( \frac{\Delta_+ \zeta_{i-1/2}^n}{\Delta z} \biggr)^2 \right)^{1/2} \\ & \leq C_T^{1/2} \Delta z^{1/2} \| \partial_z \zeta \|_{L^2 ( \Pi_T)} \quad { \to 0 \;{\rm{as}} \;\Delta z \to 0 .} \end{split} \end{align} (A.10)

    From Eq (A.8), Eq (A.9) and Eq (A.10) and appealing to the strong convergence of \phi^{\Delta z} and \psi^{\Delta z} we deduce that

    \begin{align} \lim\limits_{\Delta z \to 0} \mathcal{J}_2 = \iint_{\Pi_T} \tilde{F} ( \phi, \psi) \partial_z \zeta \, \mathrm{d} z \, \mathrm{d} t. \end{align} (A.11)

    The limits Eq (A.7) and Eq (A.11) imply that the limit \psi is a weak solution.



    [1] T. M. Apostol, Calculus, 2 Eds., Vol. Ⅱ, Blaisdell, 1969.
    [2] E. Artioli, Asymptotic homogenization of fibre-reinforced composites: a virtual element method approach, Meccanica, 53 (2018), 1187–1201. https://doi.org/10.1007/s11012-018-0818-2 doi: 10.1007/s11012-018-0818-2
    [3] E. Artioli, L. Beirão da Veiga, M. Verani, An adaptive curved virtual element method for the statistical homogenization of random fibre-reinforced composites, Finite Elem. Anal. Des., 177 (2020), 103418. https://doi.org/10.1016/j.finel.2020.103418 doi: 10.1016/j.finel.2020.103418
    [4] E. Artioli, A. Sommariva, M. Vianello, Algebraic cubature on polygonal elements with a circular edge, Comput. Math. Appl., 79 (2020), 2057–2066. https://doi.org/10.1016/j.camwa.2019.10.022 doi: 10.1016/j.camwa.2019.10.022
    [5] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, A. Russo, Basic principles of virtual element methods, Math. Models Methods Appl. Sci., 23 (2013), 199–214. https://doi.org/10.1142/S0218202512500492 doi: 10.1142/S0218202512500492
    [6] L. Beirão da Veiga, F. Brezzi, L. D. Marini, A. Russo, The Hitchhiker's guide to the virtual element method, Math. Models Methods Appl. Sci., 24 (2014), 1541–1573. https://doi.org/10.1142/S021820251440003X doi: 10.1142/S021820251440003X
    [7] L. Beirão da Veiga, A. Russo, G. Vacca, The virtual element method with curved edges, ESAIM: M2AN, 53 (2019), 375–404. https://doi.org/10.1051/m2an/2018052 doi: 10.1051/m2an/2018052
    [8] L. Beirão da Veiga, F. Brezzi, L. D. Marini, A. Russo, Polynomial preserving virtual elements with curved edges, Math. Models Methods Appl. Sci., 30 (2020), 1555–1590. https://doi.org/10.1142/S0218202520500311 doi: 10.1142/S0218202520500311
    [9] A. Bensoussan, J. Lions, G. Papanicolau, Asymptotic analysis for periodic structures, Studies in Mathematics and Its Applications, Vol. 5, North-Holland, 1978.
    [10] D. Bigoni, S. K. Serkov, M. Valentini, A. B. Movchan, Asymptotic models of dilute composites with imperfectly bonded inclusions, Int. J. Solids Struct., 35 (1998), 3239–3258. https://doi.org/10.1016/S0020-7683(97)00366-1 doi: 10.1016/S0020-7683(97)00366-1
    [11] P. J. Davis, A construction of nonnegative approximate quadratures, Math. Compt., 21 (1967), 578–582. https://doi.org/10.2307/2005001 doi: 10.2307/2005001
    [12] C. de Boor, A practical guide to splines, Springer-Verlag, 1978.
    [13] K. Deckers, A. Mougaida, H. Belhadjsalah, Algorithm 973: extended rational Fejér quadrature rules based on Chebyshev orthogonal rational functions, ACM Trans. Math. Software, 43 (2017), 1–29. https://doi.org/10.1145/3054077 doi: 10.1145/3054077
    [14] M. Dessole, M. Dell'Orto, F. Marcuzzi, The Lawson‐Hanson algorithm with deviation maximization: finite convergence and sparse recovery, Numer. Linear Algebra Appl., 30 (2023), e2490. https://doi.org/10.1002/nla.2490 doi: 10.1002/nla.2490
    [15] M. Dessole, F. Marcuzzi, M. Vianello, dCATCH–A numerical package for d-variate near G-optimal Tchakaloff regression via fast NNLS, Mathematics, 8 (2020), 1122. https://doi.org/10.3390/math8071122 doi: 10.3390/math8071122
    [16] M. Dessole, F. Marcuzzi, M. Vianello, Accelerating the Lawson-Hanson NNLS solver for large-scale Tchakaloff regression designs, Dolomites Res. Notes Approx., 13 (2020), 20–29. https://doi.org/10.14658/PUPJ-DRNA-2020-1-3 doi: 10.14658/PUPJ-DRNA-2020-1-3
    [17] T. Kanit, S. Forest, I. Galliet, V. Mounoury, D. Jeulin, Determination of the size of the representative volume element for random composites: statistical and numerical approach, Int. J. Solids Struct., 40 (2003), 3467–3679. https://doi.org/10.1016/S0020-7683(03)00143-4 doi: 10.1016/S0020-7683(03)00143-4
    [18] N. Hale, A. Townsend, Fast and accurate computation of Gauss-Legendre and Gauss-Jacobi quadrature nodes and weights, SIAM J. Sci. Comput., 35 (2013), A652–A674. https://doi.org/10.1137/120889873 doi: 10.1137/120889873
    [19] Z. Hashin, The spherical inclusion with imperfect interface, J. Appl. Mech., 58 (1991), 444–449. https://doi.org/10.1115/1.2897205 doi: 10.1115/1.2897205
    [20] C. L. Lawson, R. J. Hanson, Solving least squares problems, Classics in Applied Mathematics, SIAM, 1995.
    [21] F. Lene, D. Leguillon, Homogenized constitutive law for a partially cohesive composite material, Int. J. Solids Struct., 18 (1982), 443–458. https://doi.org/10.1016/0020-7683(82)90082-8 doi: 10.1016/0020-7683(82)90082-8
    [22] L. Mascotto, The role of stabilization in the virtual element method: a survey, Comput. Math. Appl., 151 (2023), 244–251. https://doi.org/10.1016/j.camwa.2023.09.045 doi: 10.1016/j.camwa.2023.09.045
    [23] M. Ostoja-Starzewski, Material spatial randomness: from statistical to representative volume element, Probab. Eng. Mech., 21 (2006), 112–132. https://doi.org/10.1016/j.probengmech.2005.07.007 doi: 10.1016/j.probengmech.2005.07.007
    [24] L. Piegl, W. Tiller, The NURBS book, 2 Eds., Springer-Verlag, 1997.
    [25] E. Sanchez-Palencia, Non-homogeneous media and vibration theory, Lecture Notes in Physics, Springer, 1980. https://doi.org/10.1007/3-540-10000-8
    [26] R. Sevilla, S. Fernández-Méndez, Numerical integration over 2D NURBS-shaped domains with applications to NURBS-enhanced FEM, Finite Elem. Anal. Des., 47 (2011), 1209–1220. https://doi.org/10.1016/j.finel.2011.05.011 doi: 10.1016/j.finel.2011.05.011
    [27] M. Slawski, Non-negative least squares: comparison of algorithms. Available from: https://sites.google.com/site/slawskimartin/code.
    [28] A. Sommariva, Indomain routines for NURBS, composite Bezier or bivariate parametric splines domains, alvisesommariva/inRS. Available from: https://github.com/alvisesommariva/inRS.
    [29] A. Sommariva, Software for computing algebraic cubature rules of degree n on domains defined parametrically by rational splines, alvisesommariva/CUB_RS. Available from: https://github.com/alvisesommariva/CUB_RS.
    [30] A. Sommariva, M. Vianello, inRS: Implementing the indicator function of NURBS-shaped planar domains, Appl. Math. Lett., 130 (2022), 108026. https://doi.org/10.1016/j.aml.2022.108026 doi: 10.1016/j.aml.2022.108026
    [31] A. Sommariva, M. Vianello, Low cardinality Positive Interior cubature on NURBS-shaped domains, Bit Numer. Math., 63 (2023), 22. https://doi.org/10.1007/s10543-023-00958-y doi: 10.1007/s10543-023-00958-y
    [32] A. Sommariva, M. Vianello, Computing Tchakaloff-like cubature rules on spline curvilinear polygons, Dolomit. Res. Notes Approximation, 14 (2021), 1–11.
    [33] V. Tchakaloff, Formules de cubatures mécaniques à coefficients non négatifs, Bull. Sci. Math., 81 (1957), 123–134.
    [34] D. R. Wilhelmsen, A nearest point algorithm for convex polyhedral cones and applications to positive linear approximation, Math. Comp., 30 (1976), 48–57. https://doi.org/10.1090/S0025-5718-1976-0394439-5 doi: 10.1090/S0025-5718-1976-0394439-5
    [35] D. Gunderman, K. Weiss, J. A. Evans, Spectral mesh-free quadrature for planar regions bounded by rational parametric curves, Comput.-Aided Des., 130 (2021), 102944. https://doi.org/10.1016/j.cad.2020.102944 doi: 10.1016/j.cad.2020.102944
    [36] A. Sommariva, M. Vianello, Compression of multivariate discrete measures and applications, Numer. Funct. Anal. Optim., 36 (2015), 1198–1223. https://doi.org/10.1080/01630563.2015.1062394 doi: 10.1080/01630563.2015.1062394
  • This article has been cited by:

    1. Raimund Bürger, Stefan Diehl, M Carmen Martí, Yolanda Vásquez, A degenerating convection–diffusion system modelling froth flotation with drainage, 2022, 87, 0272-4960, 1151, 10.1093/imamat/hxac033
    2. Li Feng, Yunjuan Jin, Yinzheng Sun, A class of piecewise constant Radon measure solutions to Riemann problems of compressible Euler equations with discontinuous fluxes: pressureless flow versus Chaplygin gas, 2024, 75, 0044-2275, 10.1007/s00033-024-02353-1
    3. Stefan Diehl, Jaime Manríquez, Catherine J. Paul, Tage Rosenqvist, A convection-diffusion-reaction system with discontinuous flux modelling biofilm growth in slow sand filters, 2025, 137, 0307904X, 115675, 10.1016/j.apm.2024.115675
    4. Raimund Bürger, Julio Careaga, Stefan Diehl, Romel Pineda, Numerical schemes for a moving-boundary convection-diffusion-reaction model of sequencing batch reactors, 2023, 57, 2822-7840, 2931, 10.1051/m2an/2023068
  • Reader Comments
  • © 2024 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(1426) PDF downloads(124) Cited by(1)

Figures and Tables

Figures(16)  /  Tables(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog