A difference scheme for a triangular system of conservation laws with discontinuous flux modeling three-phase flows

1.   Introduction

1.1.   Scope

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

where $t$ is time, $z$ is spatial position, and $\phi$ and $\psi$ 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 $\smash{\tilde{F}}$ are discontinuous across the positions $z = z_{\mathrm{U}} < z_{\mathrm{F}, 1} < \dots < z_{\mathrm{F}, K} < z_{\mathrm{E}}$, and due to constitutive assumptions of the model, are nonlinear functions of $\phi$ and $\psi$. The right-hand side of Eq (1.1) describes singular sources located at $z = z_{\mathrm{F}, k}$, $k = 1, \dots, K$, and is composed of given functions. It is assumed that $Q_{\mathrm{F}, k} (t)$ is the volumetric bulk flow of the mixture (of the continuous and two disperse phases) injected at $z = z_{\mathrm{F}, k}$, and that $\phi_{\mathrm{F}, k}(t)$ and $\psi_{\mathrm{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 $\Pi_T : = \mathbb{R} \times (0, T)$ together with initial conditions

where we assume that

along with

Likewise, we assume that $\phi_{\mathrm{F}, k}$ and $\psi_{\mathrm{F}, k}$ are piecewise continuous functions of bounded variation with a finite number of discontinuities and that they satisfy the bounds

(In later parts of the analysis we will assume that these functions and the bulk flows $Q_{\mathrm{F}, k}$ are constants.) If $\theta$ denotes the volume fraction of the continuous phase, then we assume that

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 $\Pi_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 $\phi$ denotes the volume fraction of bubbles and $\psi$ that of solid particles (Figure 1). The bottom of the column has the coordinate $z_\mathrm{U}$ (the underflow) and the top $z_\mathrm{E}$ (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 $\phi$ only (besides $z$ and $t$), while that of the secondary disperse phase, $\tilde{F}$, depends both on $\phi$ and $\psi$. 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.

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.,

and time-independent feed and volume rates defines Model 2, and it is shown that the corresponding scheme for the primary disperse phase (the "$\phi$-scheme" of Scheme 2; in short, "$\phi$-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")

where $J$ and $\tilde{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 "$\psi$-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 $(\phi, \psi)$ is called a weak solution of the initial value problem Eq (1.8), Eq (1.2) if

(ⅰ) The functions $\phi$ and $\psi$ belong to $L^{\infty} (\Pi_T)$.

(ⅱ) The functions $\phi$ and $\psi$ satisfy Eq (1.8), Eq (1.2) in the sense of distributions on $\Pi_T$, that is, for each smooth test function $\zeta$ with compact support in $\Pi_T$, the following identities hold:

(ⅲ) The function $\phi$ is an entropy solution of the single conservation law Eq (1.8a), that is, for each smooth and nonnegative test function $\zeta$ with compact support in $\Pi_T$, the following inequality holds for all $c \in \mathbb{R}$:

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.

1.2.   Related work

The system Eq (1.1) models the evolution of the primary unknown $\phi$ independently of the secondary unknown $\psi$. 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 $Q_\mathrm{U}$, $Q_\mathrm{F}: = Q_\mathrm{F, 1}$ (feed mixture of gas, solids and water) and $Q_\mathrm{W}: = Q_\mathrm{F2}$ (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 \times 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\times 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.

1.3.   Outline of the paper

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(\phi, z, t)$ and $\tilde{F} (\phi, \psi, 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(\phi, z, t)$ and $\tilde{F} (\phi, \psi, 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 $Q_{\mathrm{U}}$, $Q_{\mathrm{F}, k}$, $\phi_{\mathrm{F}, k}$ and $\psi_{\mathrm{F}, k}$ ($k = 1, \dots, K$). Thus, the scheme under discussion is Scheme 2. We can then prove convergence of the $\phi$-Scheme 2 (the one that discretizes the $\phi$-component of the governing PDE; Section 4.1) and $L^1$ Lipschitz continuity of the $\psi$-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 $\psi$-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) $\phi$-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 $\psi$-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.

2.   Three-phase flow model

The density of each phase is assumed constant, so the conservation of mass can be expressed by the balance equations (${v}_{\phi}$, ${v}_{\psi}$, and ${v}_{\theta}$ are the phase velocities)

where the right-hand sides contain Dirac symbols, the feed volume fractions $\phi_{\mathrm{F}, k}$ and $\psi_{\mathrm{F}, k}$ of the disperse phases, and the corresponding volume fraction $1 - \phi_{\mathrm{F}, k} (t) - \psi_{\mathrm{F}, k} (t)$ of the continuous phase, at the inlet located at $z = z_{\mathrm{F}, k}$, $k = 1, \dots, K$.

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

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

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 $z_0 < z_\mathrm{U}$ to obtain

where $\mathcal{H}(\cdot)$ is the Heaviside function. If the volumetric underflow $Q_{\mathrm{U}}(t)$ is given, then $Q(z, t) = - Q_{\mathrm{U}}(t)$ for $z < z_\mathrm{U}$, and

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

where $0 \leq \varphi \leq 1$ by Eq (1.6). The volume-average velocity of the secondary mixture is

It is then assumed that within $[z_{\mathrm{U}}, z_{\mathrm{E}})$, the relative velocity $v_{\phi \mathrm{s}} : = v_{\mathrm{\phi}} - q_{\mathrm{s}}$ of the primary disperse phase with respect to the secondary mixture is a given constitutive function $\tilde{v}_{\phi \mathrm{s}}(\phi)$, while outside that interval, both phases move at the same velocity, so their velocity difference is zero. Thus, in terms of the characteristic function

this assumption can be expressed as $v_{\phi \mathrm{s}} = \gamma (z) \tilde{v}_{\phi \mathrm{s}} (\phi)$. Within $[z_{\mathrm{U}}, z_{\mathrm{E}})$, the relative velocity of the secondary disperse phase with respect to the continuous phase $v_{\psi \theta} : = v_{\psi} - v_{\theta}$ is supposed to be a given function $\tilde{v}_{\psi \theta}$ of $\varphi$, that is, $v_{\psi \theta} = \gamma (z) \tilde{v}_{\psi \theta} (\varphi)$.

The definitions of all velocities imply the identities

for the (unweighted) fluxes $\phi v_{\mathrm{\phi}}$ and $\psi v_{\mathrm{\psi}}$ arising in Eq (2.1) and Eq (2.2), respectively. It is then useful to introduce the velocity and flux functions

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

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

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

whereas for $\phi = 1$, we set $\tilde{F}(1, \psi, z, t): = 0$ (since $F(1, \varphi, z, t) = 0$ for all $\varphi\in[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

(see ^[46]), where $v_{\mathrm{term, p}}$ and $v_{\mathrm{term, s}}$ are the terminal velocities of a single particle of the primary and secondary disperse phases, respectively, in an unbounded fluid. We set $n_{\mathrm{p}} = 3.2$, $v_{\mathrm{term, p}} = 2.7\, \mathrm{cm/s}$, $n_{\mathrm{s}} = 2.5$, and $v_{\mathrm{term, s}} = 0.5\, \mathrm{cm/s}$ along with $\sigma = -1$. These values are used in Applications 1 and 2 in Section 6. The resulting nonlinearities of $J(\phi, z, t)$ and $\tilde{F} (\phi, \psi, z, t)$ in the different zones of the column are illustrated in [52,Figure 3.2].

3.   Numerical method

3.1.   Discretization and CFL condition

The discretization of the model is based on the triangularity of the system of conservation laws Eq (1.1). The numerical fluxes for $\phi$ 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 $\phi$ 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, z_\mathrm{end})$ (to be used for the error calculation; see Section 6.1) consisting of $N$ cells by covering the vessel with $N-2$ 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 $(- \infty, 0)$ and $(z_{\mathrm{end}}, \infty)$ 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).

Given the column height $H$, we define $\Delta z: = H/(N-2)$, the cell boundaries $z_{i}: = i\Delta z$, $i = 0, 1, \ldots, N$, and the cells (intervals) $I_{i-1/2}: = [z_{i-1}, z_{i})$ and $I_{i}: = [z_{i-1/2}, z_{i+1/2})$. We place the column between $z_\mathrm{U}: = \Delta z = z_{1}$ and $z_\mathrm{E}: = z_\mathrm{U}+H = (N-1)\Delta z = z_{N-1}$. Then the length of the interval of error calculation is $z_\mathrm{end}: = H+2\Delta z = N\Delta z$. Each injection point $z_{\mathrm{F}, k}$ is assumed to belong to one cell $I_{i-1/2}$ and we define the dimensionless quantity

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

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

where

Finally, we set $t^n: = n\Delta t$ for $n = 0, 1, \ldots, N_T$.

The time-dependent feed functions are discretized as

for $k = 1, \ldots, K$, and the same is made for $\psi_{\mathrm{F}, k}$.

3.2.   Update of $\phi$-Scheme 1

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

The initial data are discretized by

To advance from $t^n$ to $t^{n+1}$ from given values $\smash{\phi_{i-1/2}^{n}}$, $i = 1, \ldots, N$, we define the numerical flux at $z = z_{i}$ by

where the notation

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

To simplify the presentation, we use the middle line of Eq (3.3) as the definition of $\mathcal{J}_{i}^n$ together with $\smash{\phi_{-1/2}^n: = 0}$ and $\smash{\phi_{N+1/2}^n: = 0}$. With the notation $\smash{\lambda: = \Delta t/\Delta z}$ and $\smash{Q_{i}^{n+}: = A_iq_{i}^{n+}}$ etc., the conservation law on $I_{i-1/2}$ implies the update formula

Then we define the piecewise constant approximate solution $\phi^{\Delta z}$ on $\mathbb{R}\times[0, T)$ by

where $\chi_\Omega$ denotes the characteristic function of the set $\Omega$.

3.3.   Update of $\psi$-Scheme 1

We discretize the initial data by

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

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

where we set $\smash{\psi_{-1/2}^n: = 0}$ and $\smash{\psi_{N+1/2}^n: = 0}$ with the same motivation as for $\phi$ above (these values are irrelevant). Here

is the Engquist-Osher numerical flux associated with the function

where ($a\wedge b: = \min\{a, b\}$, $a\vee b: = \max\{a, b\}$)

The function $\psi \mapsto \sigma f_{i}^{n} (\psi)$ is unimodal. Let $\smash{\hat{\psi}_{i}^n}$ denote the maximum point of $f_{i}^{n}$. For a given function $\tilde{V}$ the values $\smash{\hat{\psi}_{i}^n}$ and $\smash{\psi_{\mathrm{max}, i}^n}$ are related by the following lemma.

Lemma 3.1. Assume that $0 < \omega < \tilde{\omega} < 1$ are the unique local maximum and inflection point, respectively, of $f(\varphi) = \varphi V(\varphi)$ (cf. Eq (2.7)). Then $\smash{\hat{\psi}_i^n = \omega \psi_{\max, i}^n}$ for all $i$ and $n$ and all possible values $0 \leq \psi_{\max, i}^n \leq 1$. Moreover, the unique inflection point $\smash{\psi_{\mathrm{infl}, i}^n \in (\hat{\psi}_i^n, \psi_{\max, i}^n)}$ of $f_i^n$ satisfies $\smash{\psi_{\mathrm{infl}, i}^n = \tilde{\omega} \psi_{\max, i}^n}$ for all $i$ and $n$ and all possible values $0 \leq \psi_{\max, i}^n \leq 1$. (See Figure 3.)

Proof. Assume that $0 < \psi_{\max, i}^n \leq 1$. Since $\smash{\hat{\psi}_i^n}$ is the unique solution $\smash{\hat{\psi}_i^n < \psi_{\max, i}^n}$ of

it follows that $\omega$ is the unique solution in $(0, 1)$ of $\smash{\tilde{V} (\omega) + \omega \tilde{V}'(\omega) = 0}$ (cf. Eq (2.7)). By a similar argument, $\tilde{\omega}$ is the unique solution of $2 \tilde{V}' (\tilde{\omega}) + \tilde{\omega} \tilde{V}'' (\tilde{\omega}) = 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 $f_i^n (0) = 0$. For $\sigma = 1$ we get

hence

By analogous reasoning we obtain for $\sigma = -1$

We define the difference operators $\Delta_-a_i: = a_{i}-a_{i-1}$ and $\Delta_+a_i: = a_{i+1}-a_{i}$. Then the marching formula for $\psi$-Scheme 1 is

Then we define the piecewise constant approximate solution $\psi^{\Delta z}$ on $\mathbb{R}\times[0, T)$ by

3.4.   Monotonicity and invariant-region principle

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 $\phi$, Eq (3.4) (that is, $\phi$-Scheme 1) is monotone and $\smash{0\leq\phi_{i-1/2}^n \leq 1}$ for $i = 1, \ldots, N$ and $n = 1, \dots, N_T$.

Proof. We recall the assumption Eq (1.3). We first prove monotonicity of the three-point scheme for $\phi$ Eq (3.4), i.e, that $\smash{\partial\phi_{i-1/2}^{n+1}/\partial\phi_{k-1/2}^{n}\geq 0}$ for all $i = 1, \ldots, N$ and $k = i-1, i, i+1$. We have

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

We now prove that if $\smash{0\leq\phi_{i-1/2}^n\leq 1}$ for all $i$, then $\smash{0\leq\phi_{i-1/2}^{n+1}\leq 1}$ for all $i$. Clearly, Eq (1.3) implies that $\smash{0\leq\phi_{i-1/2}^0\leq 1}$ for all $i$. Since the scheme Eq (3.4) is monotone, $\mathcal H_{i-1/2}$ is non-decreasing in each argument. Since by assumption $W(1) = 0$, we get the following estimation (where we use $a^++a^- = a$):

Lemma 3.2. The function $f_{i}^{n}$ (cf. Eq (3.8)) satisfies $\|(f_{i}^{n})'\|_\infty\leq\max \{ V(0), \|V'\|_\infty \}$.

Proof. By Eq (2.7), the function $f(\varphi) = \varphi V(\varphi)$ has a single inflection point $\tilde{\omega} \in(0, 1)$ and by Lemma 3.1, $f_{i}^{n}$ has the inflection point $\tilde{\omega} \psi_{\max, i}^n \in(0, \psi_{\max, i}^n)$. We have $(f_{i}^{n})'(0) = V(0)$, $(f_{i}^{n})'(\varphi) = 0$ for $\psi_{\max, i}^n \leq \varphi \leq 1$ and the lowest (and negative) value of $(f_{i}^{n})'$ is obtained at its only critical point $\tilde{\omega} \psi_{\max, i}^n$, for which

This concludes the proof.

Lemma 3.3. There holds $\smash{G_{i}^{n} (1-\phi_{i-1/2}^{n}, 1-\phi_{i+1/2}^{n}) = 0}$ for all $i$ and $n$.

Proof. Assume that $\smash{0 < \psi_{\mathrm{max}, i}^{n} = (1-\phi_{i-1/2}^n) \wedge (1-\phi_{i+1/2}^n) \leq 1}$. By Lemma 3.1, $\hat{\psi}_{i}^n < \psi_{\mathrm{max}, i}^{n}$, hence Eq (3.10), Eq (3.11), and

imply that

Theorem 3.2. Under the assumptions of Theorem 3.1, the update formula for $\psi$, Eq (3.12) (i.e., $\psi$-Scheme 1) is monotone and along with Eq (3.4) produces approximate solutions that satisfy $\smash{0\leq\psi_{i-1/2}^n\leq 1-\phi_{i-1/2}^n}$ for all $i$ and $n$.

Proof. Assumptions Eq (1.3) and Eq (1.5) imply that $\smash{0\leq\psi_{i-1/2}^0\leq 1-\phi_{i-1/2}^0}$ for all $i$ and

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

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

Since $\smash{0 \leq \phi_{i-1/2}^n \leq 1}$ for a given $n$ and all $i$, and appealing to Eq (3.6), we have

By Eq (3.6) and Lemma 3.2 we also obtain

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

Hence, the CFL condition Eq (3.2) implies

The inequalities proved imply that $\smash{\mathcal{K}_{i-1/2}^n}$ is non-decreasing in each of its arguments. Now we use that $0\leq\psi_{i-1/2}^n\leq 1-\phi_{i-1/2}^n$ for all $i$ and Lemma 3.3 to obtain

Appealing to Eq (3.14) and the update formula for $\phi$ Eq (3.4), we get

The last equality holds since $\{ \dots \} = 0$ irrespective of whether there is a source in the cell; $Q_{i-1}^{n} - Q_{i}^{n} + Q_{\mathrm{F}, k}^n = 0$, or not; $Q_{i-1}^{n} - Q_{i}^{n} = 0$.

4.   Partial convergence analysis of Scheme 2

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 $\smash{Q_{\mathrm{F}, k}^n}$, $\smash{\phi_{\mathrm{F}, k}^n}$, and $\smash{\psi_{\mathrm{F}, k}^n}$ ($k = 1, \dots, K$) are constant and therefore do not depend on $n$. The same is assumed for the underflow volumetric flow $Q_{\mathrm{U}}$. (That is, we now study Scheme 2 suitable for Model 2.) Then Eq (3.4) and Eq (3.12) take the forms

where $q_{\mathrm{F}, k}: = {Q_{\mathrm{F}, 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 $i_{k}: = i$ if $\delta_{k, i-1/2} = 1$ (see Eq (3.1)). Then

Furthermore, we define the feed flux

such that

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

For later use we define the piecewise constant functions

where $z_{\mathrm{F}, 0}: = -\infty$, $z_{\mathrm{F}, K+1}: = \infty$, and we define the function

that allows us to write Eq (4.5) as

4.1.   Convergence of $\phi$-Scheme 2

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

with discontinuous flux

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(\phi, z)$ as given by Eq (4.9) satisfies the so-called "crossing condition" across each discontinuity

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(\phi, z^+)$ and $J (\phi, z^-)$, satisfy

which means either the graphs of $J(\cdot, z^-)$ and $J(\cdot, z^+)$ do not intersect, or if they do, there is at most one flux crossing $\phi_{\chi}$ and the graph of $J(\cdot, z^-)$ lies above that of $J(\cdot, z^+)$ to the left of $\phi_{\chi}$. For $J(\phi, z)$ as given by Eq (4.9) this condition is clearly satisfied for $z \in \{ z_{\mathrm{E}}, z_{\mathrm{U}} \}$ (considering that $j(\phi) > 0$ for $0 < \phi < 1$ implies that $J(\cdot, z^-)$ and $J(\cdot, z^+)$ do not intersect in this case), while

Thus, the crossing condition is satisfied also for $z = z_{\mathrm{F}, l}$, $l = 1, \dots, K$, since either $\phi_{\mathrm{F}, l} = 0$ and the adjacent fluxes do not intersect in $(0, 1)$, or the intersection takes place at $\phi_{\chi} = \phi_{\mathrm{F}, l}$ and Eq (4.11) holds since $q_{\mathrm{F}, 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 $z \in \mathcal{Z}$.

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

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 $\mathcal{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 $z \in \mathcal{Z}$ across which $J(\phi, z)$ is discontinuous. To state its adaptation to the situation at hand, we define the sets

Definition 4.1. A measurable function $\phi = \phi(z, t) \in L^1(\Pi_T)$ is an entropy solution of type $\mathcal{V}$ of the initial-value problem Eq (4.8), Eq (1.2a) if it satisfies the following conditions:

(ⅰ) The function $\phi$ belongs to $L^{\infty} (\Pi_T)$; for a.e. $(z, t) \in \Pi_T$ there holds $\phi(z, t) \in [0, 1]$.

(ⅱ) The function $\phi$ is a weak solution of Eq (4.8), i.e., for all smooth test functions $\zeta$ with compact support in $\Pi_T$,

(ⅲ) For all $l = 0, \dots, K+2$, for any nonnegative smooth test function $\zeta^{(l)}$ with compact support in $\smash{\Pi_T^{(l)}}$ and all $c \in [0, 1]$ there holds

(ⅳ) The following Kružkov-type ^[37] entropy inequality holds for all nonnnegative smooth test functions $\zeta$ with compact support in $\Pi_T$ and all constants $c \in \mathbb{R}$:

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 $\phi$-Scheme 2 applied to Eq (4.8), Eq (1.2a). There exists a constant $C_1$, depending on $\operatorname*{TV} (\phi_0)$, such that

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

where we define

Clearly $B_{i-1/2}^{n} \geq 0$, $C_{i+1/2}^{n} \leq 0$, and due to the CFL condition,

hence taking absolute values and summing over $i \in \mathbb{Z}$ we get, by appealing to standard arguments, that

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

which completes the proof.

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

where we define

The incremental coefficients satisfy $\smash{ C_i^n \geq 0}$ and $\smash{ D_{i}^n \geq 0}$; furthermore, the CFL condition ensures that $C_i^n + D_{i}^n \leq 1$ (in all cases for all $i$ and $n$). Notice that $\theta_i^n = 0$ with the possible exception for those indices $i$ at which $\Delta_- q_i^- \neq 0$, $\Delta_- q_i^+ \neq 0$, or $\Delta_- \gamma_i \neq 0$. According to the definition of $\gamma_i$ and that of $q_i$, see Eq (4.3), this may occur at most at a finite number of indices. Precisely, we may assert that (see Eq (4.10))

hence for all indices $i$ with the exception of finitely many indices $i$ such that $|z_j - \zeta| \leq \Delta z$ for some $\zeta \in \mathcal{Z}$, the scheme is given by the incremental form

with incremental coefficients $\smash{ C_i^n \geq 0}$, $\smash{ D_{i}^n \geq 0}$, and $C_i^n + D_{i}^n \leq 1$. 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 $V_a^b (g)$ denotes the total variation of a function $z \mapsto g(z)$ over the interval $(a, b)$.

Lemma 4.2. Consider $\phi$-Scheme 2 applied to Eq (4.8), Eq (1.2a). For any interval $[a, b]$ such that $[a, b] \cap \mathcal{Z} = \varnothing$ and any $t \in [0, T]$ there exists a total variation bound

where $C(a, b)$ is independent of $(\Delta x, \Delta t)$ and $t$ for $t \in [0, T]$.

Finally, we have shown in Theorem 3.1 that $\phi$-Scheme 1 Eq (3.4) is monotone. This applies, in particular, to the reduced $\phi$-Scheme 2 Eq (4.1) or equivalently, Eq (4.5) or Eq (4.7). Thus, $\phi$-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) ($\phi$-Scheme 2) satisfies the following entropy inequality for any $c_{i-3/2}, c_{i-1/2}, c_{ i+1/2} \in [0, 1]$:

where $h$ is defined in Eq (4.6) and the numerical entropy flux $H_i^n$ is defined by

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 $\phi^{\Delta 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 $\phi^{\Delta z}$ is defined by Eq (3.5), where the values $\smash{\phi_{i-1/2}^n}$ are defined by the scheme Eq (4.5) (that is, $\phi$-Scheme 2). Let $\Delta t, \Delta z \to 0$ with $\lambda = \Delta t / \Delta z = \mathrm{const.}$ such that the CFL condition Eq (3.2) is satisfied. Then $\phi^{\Delta z}$ converges in $L_{\mathrm{loc}}^1 (\Pi_T)$ and a.e. in $\Pi_T$ to an entropy solution of type $\mathcal{V}$ of the initial-value problem Eq (4.8), Eq (1.2a).

4.2.   $L^1$ Lipschitz continuity in time of $\psi$-Scheme 2

Next, we deal with the marching formula Eq (4.2). To this end, we define a feed flux $\smash{\tilde{h}_{\mathrm{F}, i}}$ exactly as in Eq (4.4) but with $\phi_{\mathrm{F}, i}$ replaced by $\psi_{\mathrm{F}, i}$ for $i = 1, \dots, K$. Furthermore, we recall that $\tilde{v}_{\phi \mathrm{s}} (\phi) = W(\phi)/(1-\phi)$. Thus, the scheme can be written as

Lemma 4.4 (Crandall and Tartar ^[19]). Assume that $(\Omega, \mu)$ is some measure space and that $D$ is a subset of $L^1(\Omega)$ with the property that if $u, v \in D$, then $(u \vee v) = \max\{u, v\} \in D$. Assume that $T$ is a map $T: D \ni u \mapsto T(u) \in D$ such that

Then the following statements, valid for all $u, v \in D$, are equivalent:

(ⅰ) If $u \leq v$, then $T(u) \leq T(v)$.

(ⅱ) $\int_{\Omega} ((T(u) - T(v)) \vee 0) \, \mathrm{d} \mu \leq \int_{\Omega} ((u-v) \vee 0) \, \mathrm{d} \mu$.

(ⅲ) $\int_{\Omega} |T(u) - T(v)| \, \mathrm{d} \mu \leq \int_{\Omega} |u - v | \, \mathrm{d} \mu$.

Following, for instance, ^[18], we utilize Lemma 4.4 for the following mapping. Assume that $D \subset L^1 (\mathbb{R})$ is the set of piecewise constant functions and that are constant on the intervals $I_{i-1/2}$ for $i \in \mathbb{Z}$, and that with the marching formula Eq (3.15) we associate an operator $\mathcal{K}^n: D \to D$ such that if $\psi^{\Delta z} (\cdot, t^n)$ is the piecewise constant function defined by Eq (3.13) for $t = t_n$, we may write $\psi$-Scheme 2 as

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

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

and therefore

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

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

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

5.   Convergence analysis of Scheme 3

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.,

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

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

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,

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

where we define the four-argument numerical flux

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

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

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

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

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

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

and note that

Next, we define

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

Furthermore, integration by parts yields

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

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

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

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

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

Consequently, Eq (5.5) can be written as

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

we get

which implies the inequality

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

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

we obtain from Eq (5.14)

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

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$,

Proof. We note that

We first discuss

From Eq (5.8) we get

and analogously

hence

On the other hand, we take into account that

Now

and analogously

Consequently,

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

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

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

Noticing that

we conclude that

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

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

On the other hand, in the present case

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

Lemma 3.1 (a) implies that

hence

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

Combining all possible cases we deduce that

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

where we define

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

By Lemma 3.1 (ⅱ),

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

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

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

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

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

Suppose now that

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

On the other hand,

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

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

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

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

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

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

In this case,

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

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

Furthermore, we obtain

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

and therefore

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

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

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:

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

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

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

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

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.

6.   Numerical results

6.1.   Computation of numerical error

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

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$:

with the summands defined by

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

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

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

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

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

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

6.2.   Preliminaries for numerical tests

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$.

6.3.   Simulation of a smooth solution

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.

6.4.   Illustration of the crossing condition.

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}$.

6.5.   Application 2: Counter-current fluxes

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

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.

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.

6.6.   Application 2: Co-current fluxes

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.

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.

7.   Conclusions

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.

Acknowledgements

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.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix A.   Proofs of Lemma 5.5 and Theorem 5.1

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

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

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

so we may finally write

We define the cell average

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]:

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

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

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

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

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

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}$,

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

so in the limit $\varepsilon \to 0$,

Noticing that

we get from Eq (A.3)

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

We now write

where

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

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

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

Noticing that also

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

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

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

By a computation similar to Eq (5.21) we get

where

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

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

Thus

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

and we conclude that also

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

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

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

where we define

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

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

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

hence

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

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

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