A CONSERVATION LAW WITH MULTIPLY DISCONTINUOUS FLUX MODELLING A FLOTATION COLUMN

. Flotation is a unit operation extensively used in the recovery of valuable minerals in mineral processing and related applications. Essential insight to the hydrodynamics of a ﬂotation column can be obtained by study-ing just two phases: gas and ﬂuid. To this end, the approach based on the drift-ﬂux theory, proposed in similar form by several authors, is reformulated as a one-dimensional non-linear conservation law with a multiply discontin- uous ﬂux. The unknown is the gas volume fraction as a function of height and time, and the ﬂux function depends discontinuously on spatial position due to several feed inlets. The resulting model is similar, but not equivalent, to previously studied clariﬁer-thickener models for solid-liquid separation and therefore adds a new real-world application to the ﬁeld of conservation laws with discontinuous ﬂux. Steady-state solutions are studied in detail, including their construction by applying an appropriate entropy condition across each ﬂux discontinuity. This analysis leads to operating charts and tables collect- ing all possible steady states along with some necessary conditions for their feasibility in each case. Numerical experiments show that the transient model recovers the steady states, depending on the feed rates of the diﬀerent inlets.


(Communicated by Kenneth Karlsen)
Abstract. Flotation is a unit operation extensively used in the recovery of valuable minerals in mineral processing and related applications. Essential insight to the hydrodynamics of a flotation column can be obtained by studying just two phases: gas and fluid. To this end, the approach based on the drift-flux theory, proposed in similar form by several authors, is reformulated as a one-dimensional non-linear conservation law with a multiply discontinuous flux. The unknown is the gas volume fraction as a function of height and time, and the flux function depends discontinuously on spatial position due to several feed inlets. The resulting model is similar, but not equivalent, to previously studied clarifier-thickener models for solid-liquid separation and therefore adds a new real-world application to the field of conservation laws with discontinuous flux. Steady-state solutions are studied in detail, including their construction by applying an appropriate entropy condition across each flux discontinuity. This analysis leads to operating charts and tables collecting all possible steady states along with some necessary conditions for their feasibility in each case. Numerical experiments show that the transient model recovers the steady states, depending on the feed rates of the different inlets.

Introduction.
1.1. Scope. Flotation is a unit operation that is extensively used in the recovery of valuable minerals and coals in mineral processing but also in many other applications in environmental and chemical engineering [12,21,33,36,42]. It is a physico-chemical separation process that utilizes the difference in surface properties of the valuable hydrophobic minerals and the unwanted hydrophilic gangue material. The theory of froth flotation is complex and involves three phases (solids, water, and froth or gas) with many subprocesses [42]. The principle of the conventional flotation process is roughly as follows. Gas is introduced close to the bottom   [21,38]), including heights of singular sources z G , z F and z W , the underflow level z U , and the effluent level z E . Right: Corresponding conceptual model of the flotation column used in this work, indicating the volumetric feed flows Q G , Q F and Q W , the underflow rate Q U , the effluent rate Q E , and the spatially piecewise constant bulk velocity q = q(z, t).
of a flotation column (see Figure 1), and the bubbles generated then rise upward throughout the pulp that contains hydrophobic and hydrophilic solid particles. The hydrophobic particles in the pulp attach to the bubbles. Since the overall density of the bubble-particle aggregates is less than that of the medium, the aggregates then float to the top of the column, where the desired product, the foam or froth carrying the valuable material (the concentrate in mining) is removed, usually through a launder. Additionally, close to the top wash water is injected to assist with the rejection of entrained impurities [39] and to increase the froth stability and improve recovery [21,31]. Once the hydrophobic particles have attached to the air bubbles, flotation can be considered as a separation between relatively large low-density entities, called air bubbles, and a suspension of liquid and gangue. Consequently, flotation can be described as a gas-liquid separation process by buoyancy analogous to the solid-liquid separation by gravity sedimentation in clarifier-thickeners [8,10,15].
Well-established spatially one-dimensional models of clarifier-thickeners can be formulated as a scalar conservation law for the local solids concentration as a function of depth and time, where the flux is discontinuous as a function of spatial position due to upward-and downward-directed bulk flows, transitions to overflow and underflow transport, and a singular source term marking the feed [8,10,15]. Clarifier-thickener models have motivated in part the mathematical research on conservation laws with discontinuous flux [3,4,6,10,[14][15][16][17][18][19][20].
It is the purpose of this paper to formulate, partially analyze, solve for steady states, and numerically simulate a related model for a flotation column, where we follow [13] and limit ourselves to a one-dimensional two-phase system of gas bubbles dispersed in a fluid, or rather a suspension of liquid and gangue. Hence, we do not model any sedimentation of solid particles in the suspension. The final form of the model (cf. Figure 1) is the conservation law with multiply discontinuous flux where φ is the sought volume fraction of gas bubbles as a function of height z and time t, and F (z, t, φ) is a flux function, made precise in later parts of the paper, that is nonlinear in φ and discontinuous in z at five different positions. Three of these positions, z G , z F and z W , are associated with singular feed sources of gas, feed slurry, and wash water, respectively, at given rates q S and volume fractions φ S , where δ(·) denotes the Dirac symbol. The model (1.1) is posed for z ∈ R without boundary conditions, and is therefore supplied solely with initial data While the time-dependent partial differential equation (PDE) (1.1) describes transient variations of φ as a function of position and time, a property of practical interest in applications are the stationary solutions to the model that correspond to undisturbed normal states of operation. A steady-state solution of (1.1) is generally a piecewise constant function, with possible jump discontinuities both within the four zones of Figure 1, and across the five spatial discontinuities z = z E , etc.
The main outcomes of this work are to a classification of all steady-state solutions by means of diagrams and tables, and numerical simulations of dynamic behaviour. The variety of real-world applications of conservation laws with discontinuous flux is hereby widened to include flotation.

1.2.
Related work. Our model formulation is based on the description of a flotation column by Stevenson et al. [38], Dickinson and Galvin [13], and Galvin and Dickinson [22] that is based on algebraic expressions for the gas and liquid fluxes, velocities and volume fractions. The description of one-dimensional two-phase flows based on the continuity equations for both phases and closed by defining a relative flux, or drift flux, between both phases as a function of volume fraction was introduced by Wallis [40], as is elaborated in [35]. Treatments that invoke this drift-flux analysis to describe flotation processes include [25,27,31,39,43,44]. However, all these works utilize these variables for steady-state analyses, but do not incorporate the drift-flux variables into one solvable PDE model for transient simulations, which is precisely the main contribution of the present paper.
As stated above, the theory of conservation laws with discontinuous flux has seen a vast amount of interest in recent years, where the typical model equation is or equivalently, in terms of the Heaviside step function H(z), The basic difficulty associated with (1.3) is as follows. Suppose, for simplicity, that φ = φ(z, t), z ∈ R, t > 0 is a piecewise constant solution of (1.3) having traces φ − (t) := lim z↑0 φ(z, t) and φ + (t) := lim z↓0 φ(z, t) at z = 0. Then the fluxes to both sides of z = 0 must be equal at any time, which compels the jump condition This single equation does not define the two traces uniquely and one needs to specify a selection principle or jump entropy condition to single out pairs that besides satisfying (1.4) are admissible. This selection principle usually depends on the particular physical reality (1.3) is supposed to model. For instance, applications of (1.3) also include traffic flow with discontinuously changing road surface conditions, ion etching, two-phase flow in heterogeneous porous media, and medical applications (see [5], [24,Ch. 8], and [28] for overviews and references). We use here the admissibility condition from [14], which has proved to be the natural one for the related problem of continuous sedimentation [15]. Furthermore, its generalization [18] to the case of a scalar convection-diffusion equation with spatial discontinuity in both the flux and diffusion functions implies the physically relevant solution in the case of the well-established model of continuous sedimentation with compression [10].
As is stated in [24, p. 426], there are different "recipes" to select unique solutions of the Riemann problem of (1.3), and all of them eventually lead to uniqueness of the initial value problem for (1.3), according to the unified treatment in [1].
1.3. Outline of the paper. In Section 2, we derive the model equations for the local fraction of gas bubbles, detailing the definition of the flux density functions in each zone of the spatial domain and the treatment of the feed inlets. Some notation necessary for the description of the steady-state solutions is also introduced. In Section 3, we focus on the characterization of all possible steady states for the flotation model previously defined, providing a detailed study of the derivation process at the different spatial discontinuities introduced by the feed inlets. Some steady states exist only under certain conditions on the injection rates and to get an overview of all possibilities, we present operating charts and tables for the categorization of all steady states. In Section 4, we briefly review the numerical method to approximately solve the flotation model. Some numerical examples are provided in Section 5 and, finally, we present some conclusions in Section 6.

2.1.
Phase velocities and drift flux. Assume that we consider a region of space that is free of sources and sinks, that φ is the local fraction of gas bubbles, and that v f and v g are the phase velocities of the gas and fluid, respectively. Then the conservation of mass equations for both can be written in local form as ∂φ ∂t where we assume that the gas bubbles are incompressible and do not coalesce. Then, defining the volume average velocity, or bulk flux of the suspension, and the gas-fluid relative velocity v r := v g − v f , we may rewrite the first equation in (2.1) and replace the second by the sum of both to obtain ∂φ ∂t It is assumed that v r = v r e z , where e z is the upward-pointing unit vector and v r is given as a function of φ so that the gas drift-flux function is 3) The drift-flux function j g (φ) gives the gas flux in a closed column relative to the column. Here, v term is the terminal velocity of a single bubble in an unbounded fluid. As is stated in [38], there exists a number of methods to calculate v term , and Wallis' generalized correlation [41] leading to v term is recommended, see [38, Appendix A] for details. This correlation involves additional quantities such as equilibrium surface tension and the viscosity of the fluid. Its detailed discussion is beyond the focus of this paper; for the present analysis it suffices to assume that v term is a set constant for a given material and equipment. Furthermore, in (2.3), V (φ) is a dimensionless hindered-bubbling function that can, for instance, be given by the Richardson-Zaki expression [34] The maximum possible volume fraction of bubbles is φ max = 1. Realistic values of n RZ range from n RZ = 2 to n RZ = 3.2 (cf., e.g., [13,22,25,31,39]). Finally, in one space dimension (in the z-direction) and away from sinks or sources, ∇ · q = 0 reduces to ∂q/∂z = 0, so q will depend on t only, and the only equation that needs to be solved is the nonlinear first-order conservation law Hence, j(φ, t) is the total flux of gas. If we denote the total fluid flux by ϕ := (1 − φ)v f , then (2.2) yields the simple relation q = j + ϕ between the three fluxes, which all generally may have any sign.

2.2.
Volumetric flows, bulk velocities and flux functions. It is assumed that the unit has a constant cross-sectional area A, and that concentrations are horizontally constant so that all variables depend on height z and time t only. The unit is fed at heights z = z W , z = z F and z = z G , with wash water, fluid (feed slurry), and gas, respectively (see Figure 1), where we assume that z W > z F > z G . The corresponding volumetric flows Q W ≥ 0, Q F ≥ 0 and Q G ≥ 0 are assumed to be given functions of time, as is the volumetric underflow rate Q U ≥ 0. Furthermore, Q E is the resulting effluent volumetric overflow, which is assumed to be nonnegative, i.e., the mixture is conserved and the vessel is always completely filled with mixture. The spatially piecewise constant bulk velocity q = q(z, t) is defined by the values q 1 (t) to q 4 (t) in the corresponding four zones in the vessel; see Figure 1. To simplify notation, we define the velocities q S (t) := Q S (t)/A for S ∈ {E, F, G, U, W}.
Starting from the bottom of the vessel, we have q 1 = −q U , q 2 = q 1 + q G , etc. and we obtain the total bulk velocity function Hence, we always have q E = q 4 ≥ 0 and q 1 = −q U ≤ 0, whereas q 2 and q 3 may have any sign.
We denote the intervals z > z E and z < z U by the effluent and underflow zone, respectively. During normal operation, zone 1 contains only liquid. Above z = z G there is a region of bubbles rising with a high velocity up to an upper region of froth with a high volume fraction of gas; see Figure 1. The large discontinuity between the bubbly and froth regions is usually located in the interval [z F , z W ], either at the location of one injection point or within a zone. As we will see, discontinuities in φ are possible at every injection point and inside some zones. Usually, zone 4 is small, and sometimes z G = z F , i.e., gas and feed slurry are injected at the same location.
We assume that in the effluent and underflow zones, the gas and the fluid move at the same velocity, so we set v r := 0, and therefore j g := 0 in these zones. The total flux function is then defined as (2.7) 2.3. Feed sources and governing equation in final form. The conservation law (2.5) is completed by including the feed of material at levels z W , z F and z G at volume rates Q W , Q F and Q G . The feed mechanisms give rise to singular source terms that extend (2.5) to the final governing model equation (1.1). The given gas volume fraction of the three time-dependent feed streams is denoted by φ W , φ F and φ G , respectively. In agreement to common practice, we assume that either pure gas or pure liquid is injected through the respective singular sources, so we assume from now on φ G ≡ 1 and φ F = φ W ≡ 0.

Entropy solutions.
Within each zone, the governing equation (1.1), (2.7) reduces to (2.5), and we consider the Cauchy problem of this equation. A piecewise smooth function φ = φ(x, t) is defined to be an entropy solution of the problem if φ is continuously differentiable everywhere with the exception of a finite number of and the jump entropy condition It is well known that entropy solutions in the sense of Oleinik [30] are also the unique entropy solutions in the sense of Kružkov-type [26] integral inequalities (cf., e.g., [24]). On the other hand and as mentioned in Section 1.2, at the five spatial discontinuities of problem (1.1), a generalized entropy is needed [14,18]. Since we only construct steady-state solutions, we review that condition in Section 3.1 for this purpose, which means less notation than for the dynamic case.

2.5.
Properties of the zone flux functions. We assume that the drift-flux function j g (φ) ≥ 0 is continuously differentiable, satisfies j g (0) = j g (1) = 0, and is concave-convex with an inflection point φ infl ∈ (0, 1); see  These values are also utilized in [13].) The form of j g singles out certain distinguished values of the volume fraction that appear in the steady-state solution (see Figure 2). The same values appear in the related problem of continuous sedimentation [16,17].
Since we will mostly refer to steady-state situations, we often supress the dependence on t in j 1 , . . . , j 4 . The functions j 1 , . . . , j 4 and j g have the same inflection point φ infl , since they only differ by a linear term. If j k (φ) has a zero in (0, 1], which can happen only for k = 1, 2, 3 and q k < 0, then we denote it by φ kZ . If j k (φ) < 0 for all φ ∈ (0, 1], we set φ kZ := 0. We definē which are the bulk velocities such that the slope of j k (φ) is zero at φ max = 1 and φ infl , respectively. To reduce the number of cases to investigate, we assume in this work thatq To obtain a definition for all values of q k , we note that the restriction (j g | (φ infl ,1) ) is a decreasing function so that we can define Given φ kM and q k ≥ 0, we define φ km as the unique value satisfying For q k <q, j k (φ) assumes a local maximum point at φ M k ∈ [0, φ infl ). Let q neg := −j g (0) be the value below which j k is a decreasing function. For q k ≤ q neg , the local maximum is φ M k := 0. For q k ≥q, we set φ M k := φ infl . In some instances it is convenient to write out the dependence on q of the flux function, i.e., j k (φ; q k ), and of the specific concentrations, e.g., φ M k (q k ). The following properties follow directly from the definitions above (cf. [16,Lemma 2]).
Lemma 2.1. The following properties hold: 3. Steady-state solutions. In order to completely describe all possible steady states of model (1.1) (under the assumptions of Section 2.5), we here extract from the theory of conservation laws with discontinuous flux function [14,18] what is necessary to construct steady-state solutions in a neighbourhood of a spatial discontinuity.
To start with, we investigate the case when the steady-state solution is constant φ k in zone k = 1, 2, 3, 4, and constant φ U and φ E in the underflow and effluent zones, respectively. In Section 3.9, we describe the general solution, having also discontinuities within one or several zones. First, we review the necessary theory and notation, and then go through the possible couplings between zones 1 to 4 and the effluent and underflow zones.
3.1. Construction of steady-state solutions for a conservation law with discontinuous flux function. We consider the conservation law with discontinuous flux function (1.3). The equation should be interpreted in the weak sense and we seek steady-state solutions of the form where φ ± are constants. The conservation law across z = 0 implies the jump condition g(φ − ) = f (φ + ) (see (1.4)). This single equation has two unknowns. The generalized entropy condition [14] selects the physically relevant solution in a neighbourhood of z = 0 for a dynamic solution of (1.3) for given initial data at t = 0. We define the auxiliary functionŝ Sinceǧ(·; φ − ) is non-increasing andf (·; φ + ) is non-decreasing, the intersection of the graphs of these functions occurs at a unique flux value, if there exists an intersection. For the model of flotation, this is always the case; the proof of this statement can be made in the same way as for the problem of continuous sedimentation; see [15]. We define the set of possible φ-values of the intersection as and the corresponding unique flux value η(φ + , φ − ) :=f (Φ; φ + ). Since we are here only interested in stationary solutions, the generalized entropy condition can be stated asf

3.2.
Couplings at the underflow location z = z U . In a neighbourhood of z = z U , the PDE (1.1) becomes We now suppress the time dependence and seek possible constant solutions φ 1 and φ U in the two neighbouring zones so that the entropy condition (3.1) is satisfied. With the notation of Section 3.1, we have . The task is to find a pair φ U and φ 1 so that the entropy condition (3.1) is satisfied. This condition now readŝ where η(φ 1 , φ U ) denotes the flux value of the intersection of U (·; φ U ) and 1 (·; φ 1 ). As j U (φ) = −q U φ is a linear decreasing function, we have U (·; φ U ) = j U for any φ U ∈ [0, 1]. The function j 1 has two monotonicity intervals separated by the maximum point φ M 1 , which leads to the following cases: Fig. 3(a). The only possible intersection between 1 and U is φ 1 = 0; hence φ U = 0. Thus, the zero solution on both sides is the only possible steady-state coupling in this case.
), the middle plot shows that the only possible intersection is, as in the previous case, φ 1 = 0, but is outside the interval of definition of φ 1 , i.e., there is no possible steady state We conclude this subsection by stating the possible steady-state values in the underflow and the first zone: 3.3. Couplings at the effluent location z = z E . Here, j E (φ) = q E φ is an increasing linear function. The procedure is analogous to the coupling at z = z U . In fact, this case is the same as the one at the underflow level in the problem of continuous sedimentation. Details can be found in [15,Section 9], and we state here directly the possible steady-state values in zone 4 and the effluent zone:  Figure 3. The decreasing function U (·; φ U ) = j U and three possible cases of graphs of 1 (·; φ 1 ) depending on φ 1 . The intersection of 1 (·; φ 1 ) and j U defines the possible values in a steady-state solution.

Couplings at
The flux functions to consider for the steady-state coupling are thus which intersect only at φ = 1, and the entropy condition iŝ We have q 1 = −q U ≤ 0, but q 2 may have any sign. We will make a main division depending on the sign of q 2 . From (3.3), we know that a steady state in zone 1 satisfies φ 1 ∈ {0} ∪ [φ 1Z , 1]. These two intervals should be coupled with the three monotonicity intervals of j 2 , which are separated by φ M 2 and φ 2M . We will get conditions on q G in the subcases.

Remark 1.
In the division into subcases further on, we will sometimes let such overlap in the following way. Instead of having disjoint intervals defining two subcases, e.g., . This overlap will reduce the number of conditions in terms of inequalities. Then the same steady-state solution can occur in two subcases, but the final number of steady-state solutions is not affected. Fig. 4(a). A necessary condition for a steadystate solution is Fig. 4(b). A necessary condition for a steadystate solution is (G).   Fig. 4(e). A steady-state coupling is possible with φ 2 ≤ φ 1 with equality if and only if q 1 = q 2 , i.e., q G = 0. Fig. 4(f). The only possible coupling is Fig. 5(a). The same as Case G1(a) with condition (G). Fig. 5(b1) and (b2). As in Case G1(b) a necessary condition is (G). The largest value in zone 2 is φ 2 = φ 2Z and then q 1 = q 2 , i.e., q G = 0, see plot (b2). Fig. 5(c). The same as Case G1(d): Solution exists with either a positive, see Fig. 5(d1) or a negative flux (d2).

3.5.
Couplings at z = z F . We derive the possible constant steady states in zones 2 and 3 considering their coupling at z = z F . This is the most complicated case where both bulk velocities, hence both fluxes, can be either positive or negative. Since there is no injection of bubbles, φ F = 0, the fluxes to consider are We require that the entropy condition (3.1) holds, which now reads We now study the possible intersections between 2 (·; φ 2 ) and 3 (·; φ 3 ) and make a division into three main cases F1-F3 depending on the signs of the two bulk velocities q 2 and q 3 . Subdivisions are then made according to the intervals of monotonicity of the flux functions j 2 and j 3 , which give qualitatively different intersections of 2 (·; φ 2 ) and 3 (·; φ 3 ) because of their strictly monotone parts and plateaus. Some cases will be empty, i.e., no steady-state solution is possible and this may depend on q 2 and q 3 . When a steady-state solution is possible, it always satisfies Fig. 6(a). In this case there is a possible intersection point between 2 (·; φ 2 ) and 3 (·; φ 3 ). We have The flux value of the intersection of 2 and . Since this case requires φ M 2 < φ 2 , (3.5) is not satisfied and there is no possible stationary solution. Fig. 6(c). The only possibility for the plateau of 2 (·; φ 2 ) to intersect the increasing part of 3 (·; φ 3 ) is that the plateau lies above the value of the local maximum of j 2 . In other words, this case is empty unless the following condition is satisfied: Fig. 6(d). For (3.5) to be satisfied, two plateaus are involved in the intersection. A necessary condition for this is that the flux value of the local maximum of j 2 is larger than or equal to the flux value of the local minimum of j 3 , i.e., ]. An intersection occurs between the decreasing part of 2 (·; φ 2 ) and the plateau of 3 (·; φ 3 ) only if (FII) holds.
Two plateaus are necessarily involved in the intersection, which occurs only if the following condition is satisfied: The intersection occurs at the flux level is not allowed in this subcase, so (3.5) is not satisfied and this case is therefore empty.
. This case is empty with the same motivation as in subcase (g).
Similarly to subcase (f), a steady-state solution is possible only if (FIII) holds. We have q 2 ≤ q 3 ⇔ φ 2 ≤ φ 3 . 2. Case F2. q 2 ≤ 0 ≤ q 3 . We use the same principle of subdivision as in Case F1. However, when q 2 ≤ 0, the local minimizer of j 2 (·; q 2 ) is φ 2M = 1, so the domain of j 2 (·; q 2 ) consists of two disjoint intervals where j 2 is monotone. Together with the three intervals of j 3 , we get six subcases.
. This case is empty with analogous arguments as in Case F1 The intersection of 2 (·; φ 2 ) and 3 (·; φ 3 ) is qualitatively the same as in Case F1(d) and (e) and the steady state represented in Fig. 7 By analogy with subcase (c), (FII) has to be satisfied.
. This case is analogous to Case F1(g), and is consequently empty.
. Analogous to Case F1(h) and empty case. 3. Case F3. q 2 ≤ q 3 ≤ 0. In this case φ 2M = φ 3M = 1, so there are only four subcases with the maximum point of each flux function as the point of division of two monotonicity intervals.
. This is qualitatively the same as Case F2(a) with steady-state solution possible.
. This case is empty with analogous arguments as in Case F1(b) or F2(b) (no plot is shown).  Figure 6. Case F1: 0 ≤ q 2 ≤ q 3 . Possible steady-state values for zones 2 and 3. In the special case q 2 = q 3 , the diagonal plots (a), (e) and (i) are the only ones where φ 2 = φ 3 occurs.
There is a possible intersection only at a non-negative flux value. Then bubbles flow upwards and we have Intersection may occur at a positive or negative flux value (bubbles moving downwards or upwards), see Fig. 8(d1) and (d2), respectively. 3.6. Couplings at z = z W . Since φ W = 0, the fluxes to consider for the entropy condition are 3 (·; φ 3 ) and 4 (·; φ 4 ). From (3.4) we know that only φ 4 ∈ [0, φ 4m ] ∪ [φ 4M , 1] are possible steady states in zone 4. In these two intervals, j 4 is increasing and we combine these intervals with the three monotone parts of j 3 . Fig. 9(a).
] is an empty case (no plot is shown). (c) φ 3 ∈ [φ 3M , 1], φ 4 ∈ [0, φ 4m ] is an empty case (no plot is shown).  Necessarily, φ 4 = φ 4M holds. Another necessary condition for this steady state is that the flux value of the local maximum of j 3 is larger than or equal to the flux value of the local minimum of j 4 , i.e., Fig. 9(e). In this case φ 4 = φ 4M and the necessary condition is (WI). Fig. 9(f). This case is empty unless the following condition holds: The conclusions are the same as in Case W1.

3.7.
Operating charts for bulk velocities. The different necessary conditions on the fluxes that appear in the derivation of possible steady states can be visualized in operating charts. These are two-dimensional diagrams involving the bulk velocities at an injection point. Condition (G) can be written as where the dependence on q 2 is written out as in Lemma 2.1. This lemma gives that j 2 (φ M 2 (q 2 ); q 2 ) is an increasing function of intermediate (normal) values of q 2 and otherwise constant (for q 2 < q neg and q 2 >q). Its graph, and consequently the region in (q 2 , q G )-space where (G) is satisfied, are shown in Fig. 10(a). In Fig. 10(b), the corresponding region in (q U , q G )-space is shown. This has been obtained by the linear mapping of the curve q G = j 2 (φ M 2 (q 2 ); q 2 ) (since q U = q G − q 2 ): We will now do the same for the coupling at z = z F and conditions (FI)-(FIII). These lead to overlapping regions in the (q 2 , q 3 )-plane leading to different numbers of possible steady states in different regions. In Fig. 11, the regions are shadowed in which (FI) (red), (FII) (blue) and (FIII) (grey) hold.
To obtain these regions, we define the following functions with respect to conditions (FI)-(FII): These functions are continuously differentiable by Lemma 2.1. The following lemma gives the qualitative properties of the boundaries of the regions. Lemma 3.1. There exists a uniqueq ∈ (0,q) and a unique continuously differentiable functionh II such thatq = j 3 (φ 3M (q);q) and the following hold: Proof. Using Lemma 2.1, we get  where the last inequality follows from the fact that j 2 (·;q) is an increasing function. Hence, h I is a continuous increasing function and (3.6) follows. For condition (FII), we have (by means of Lemma 2.1) Since the latter derivative is always non-zero, the implicit function theorem implies the existence of a continuously differentiable functionh II satisfying For q 2 ≤ q neg and q 3 = 0, we have and for q 3 = q 2 ≥q, we have Condition (FIII) can directly be rewritten by using the identity j 2 (1; q 2 ) = q 2 . Hence, the boundary of the region where condition (FIII) holds is given by the curve j 3 (φ 3M (q 3 ); q 3 ), and its properties follow from Lemma 2.1 and the definition of φ 3M (q 3 ). Finally, the definition ofq gives that As we did for the gas injection point, we now transform the boundary curves of conditions (FI)-(FIII) from the (q 2 , q 3 )-plane to the (q G , q F )-plane with the linear mapping This mapping depends on the chosen steady-state value of q U , denoted by q SS U , from the chart in Fig. 10(b). For a given value of q SS U , we therefore use the scale q G − q SS U on the horizontal axis; see Fig. 11(b).
For the coupling at z = z W , we note that conditions (WI) and (WII) are the same as conditions (FII) and (FIII) except for the zone index increased by one. Hence, the analogous statements of Lemma 3.1 hold, except for those involvingq. We can draw the operating chart in the (q 3 , q 4 )-plane, and with the mapping the chart in the (q F , q W )-plane, where the steady-state values q SS U and q SS G have been chosen from the previous operating charts; see Fig. 12.
It seems logical that when working with steady states, it is necessary to restrict the amount of gas, fluid or water pumped into the tank. Moreover, from the desliming point of view, situations as huge quantities of fluid leaving at the top of the column, loosening the froth or washing it in excess, or gas bubbles flowing down and out of the tank through the bottom tailings underflow, are not convenient. Conditions (G), (FI)-(FIII) and (WI)-(WII) set a theoretical upper limit for the values of q G , q F and q W , respectively, for a fixed given value of q SS U , as Figures 10-12 show. For instance, Fig. 11(b) shows that if we want any condition (FI)-(FIII) to be satisfied for q G − q SS U ∈ [−0.1, 0.9] cm/s, then we should set q F < 0.7 cm/s. 3.8. Tables and visualization of steady states. Whereas the operating charts in Section 3.7 give an overview on how to choose the bulk velocities with respect to conditions for certain steady-state couplings at the points of injection and outlets, we here collect in Tables 1 and 3 all possible steady-state combinations between these couplings for the cases q 3 ≥ q 2 ≥ 0 and q 2 < 0, respectively. Some of the couplings made locally, for example for φ 2 and φ 3 in Section 3.5, do not appear in the tables because they are not possible when taking into account also the couplings with φ 1 and φ 4 . The tables should be read as follows. A possible steady-state solution with constant concentration in each zone is obtained by connecting adjacent rectangles passing only over vertical lines (and no corners), from the left column φ U to the right φ E . In Table 2 we show some examples of admissible and inadmissible paths, in green and red respectively, for steady states for the case q 3 ≥ q 2 ≥ 0 in Table 1. As it can be seen, the inadmissible path passes over a horizontal line in the φ 2column, yielding two different possible intervals of definition for φ 2 , which is not allowed. Table 3 shows the possible steady states when q 2 < 0 and the sign restrictions of q 3 are given in the table. The appearance of conditions (G), (FI), etc. in the tables means that the corresponding steady state is possible only if the corresponding conditions are satisfied. As can be seen in the first and the last columns of the tables, the values of φ U and φ E are uniquely given by the chosen value of φ 1 and φ 4 , respectively. Also notice that although in the derivation of the cases we allow the intervals of definition of φ i to overlap, we do not do this in the construction of the tables, to avoid the appearance of a possible steady state more than once and keep the tables as simple as possible.
3.9. Steady states having discontinuities within zones. Despite the diversity of possible steady states with a constant value in each zone, which are categorized  Table 2. Admissible (green) and inadmissible (red) paths for steady-state construction in Table 1. in Tables 1 and 3, there exist further steady states with possible stationary discontinuities within the zones. For example, in zone 1, the constant solutions φ 1 = 0 and φ 1Z have the same flux value, i.e., j 1 (0) = j 1 (φ 1Z ). This means that there may be a stationary discontinuity from φ 1 = 0 below to φ 1Z above the discontinuity, which may be located anywhere in the interval (z U , z G ), since the entropy condition (2.9) is satisfied. If such a discontinuity exists in zone 1, then the value φ 1Z can be coupled with admissible values of φ 2 according to Table 1 or 3. Table 3. Collection of possible steady-states for the flotation column when In the same way, stationary discontinuities from a lower to a higher value are possible in zones 2, 3 and 4 for all bulk velocities q k that imply that the drift-flux function in the zone is not monotone; hence, q k > q neg should hold. Analogously, there may be discontinuities from higher gas volume fraction below a standing discontinuity than above, if the drift flux has its local minimum point φ kM (q k ) ∈ (φ infl , 1), which is equivalent to q k ∈ (q,q). We provide examples of such discontinuities in Section 5.
4. Numerical method. For the discretization of the model, we follow the procedure of [3] for the sedimentation process in a clarifier-thickener.
We subdivide the tank into N internal layers, or cells, and four external layers, two at the top and two at the bottom, corresponding to the overflow and underflow zones, respectively, each of them with depth ∆z = 1/N . We let φ i (t) ≈ φ(z i , t) denote the average of the exact solution over layer i, i.e. (z i−1 , z i ), at time t, i.e., For each layer i, we use the balance law corresponding to (1.1) to obtain where the flux F (z, t, φ) is defined by (2.7). This is approximated by Godunov's numerical flux [23], which is where k ∈ {E, 1, 2, 3, 4, U} denotes the present zone of the vessel (k depends on i). We define the indices i G , i F , i W to indicate that the gas inlet is located in the layer (z iG−1 , z iG ), and the fluid and wash water inlets in (z iF−1 , z iF ) and (z iW−1 , z iW ), respectively. For instance, the numerical fluxes F i between layers (z iG−1 , z iG ) and (z iF−1 , z iF ) will be computed using the function j 2 (φ).
Substituting the numerical fluxes into the exact version of the conservation law (4.1), we get the following method-of-lines formula: we obtain the following fully discrete method for φ i,n ≈ φ i (t n ), where the upper index n stands for evaluation at time t = t n : With these values, the velocities in the column are and hence q U = q 1 = −0.1 and q E = q 4 = 0.2353. Using Table 1 we know that, for these values satisfying conditions (G), (FII), (FIII), (WI) and (WII), there are eight possible steady states with a constant volume fraction in each zone. In Table 4, the admissible paths for the case φ 1 = 0 can be seen, corresponding to steady states (a)-(d) in Figure 13. The paths with φ 1 ∈ [φ 1Z , 1] can be obtained analogously, starting from the corresponding rectangle in Table 1. A detailed description of all steady states is provided by Figures 13-15. Table 4. Examples 1 and 2: admissible paths for steady states with φ 1 = 0, corresponding to steady states (a)-(d) in Figure 13.  In addition to the volume fractions φ k , k ∈ {U, 1, 2, 3, 4, E}, Figure 15 shows the gas and fluid fluxes (j k and ϕ k ) in each zone k represented by red and blue arrows, respectively. In all eight cases, the gas flux is zero below z = z G , since j 1 (φ 1 ) = 0, and equal to the constant value j 4 (φ 4M ) above. The flux of water in zone k is then Let us now run two simulations.
Example 1. Initially, we consider a column filled with only fluid, i.e., φ(z, 0) = 0 for all z, and we start pumping gas, fluid and wash water into it, with the velocities stated above. The gas rises fast and reaches the top of the vessel in a short time; see Figure 16, which shows the time evolution of the gas concentration. As can be seen in Figure 17(a), at t = 150, the first steady state, see Figure 15(a), has been reached. Then we close the top of the vessel by setting q E = 0 and q U = −0.3353 and let the gas accumulate at the top until t = 300, see Figure 17(b), when we open the top again (the previous values are used). As can be seen in Figures 17(c) and (d), after opening the top, a discontinuity arises within zone 3 between the fluid and wash water inlets, and becomes stationary at z ≈ 55; see the explanation in Section 3.9. The new steady state is thus not among the eight ones shown in Figure 15, but can be seen as combination of cases (b) and (c) of both Figures 13  and 15.
Example 2. Figures 18 and 19 show the results when we close the top of the tank for a longer period, in this case until t = 400. Then, still another steady state is reached, namely the intermediate one of Figure 15(d) and (h) having a discontinuity in zone 1. The flux of gas is zero in zone 1, which means that there is a region of bubbles standing still with the concentration φ ≈ 0.65 just below the gas inlet z G = 25.  5.2. Example 3: Overloaded tank. We now define an initial set of velocities for the inlets and the outlets for which no steady state with φ 1 = 0 is possible. Consequently, the system is overloaded and gas has to leave also through the underflow, so that φ 1 > 0. The velocities in the column chosen are and hence q U = q 1 = −0.1 and q E = q 4 = 0.17. A detailed description of the possible steady states is provided by Figure 20.
As in Examples 1 and 2, we start with a column filled with only fluid, i.e., φ(z, 0) = 0 for all z, and we start pumping gas, fluid and wash water into it, with the velocities stated above. In Figure 21, we can see how the gas rises and reaches the top of the column in a short time. A steady state, corresponding to case (d) in Figure 20, has been reached by t = 2000, as it can be seen in Figure 21(d). Clearly this is not an acceptable set of parameters for the desliming process: an unnecessary excess of gas is flowing downwards and leaving the tank through the bottom.

Example 4: Desliming.
A flotation column is generally operated so that there are two regions: one with bubbles (intermediate gas volume fraction) and one with froth. In the bubbly region, usually in zone 2 and sometimes the lower part of zone 3, the hydrophobic particles of the pulp slurry attach to the air bubbles. In the froth region, located above the bubbly region, further enrichment takes place and the foam is efficient for promoting water rejection. The injection of wash water several (but finite number of) discontinuities. In [10], a Kružkov-type of entropy condition was used together with a crossing condition for the fluxes for the proof of uniqueness. It is worth noting that the flux discontinuities within the present model do satisfy this condition. The entropy condition used here implies, however, the Kružkov-type and uniqueness is obtained without the crossing condition [18].
One-dimensional models such as the one treated herein are easier to solve than multi-dimensional multiphase flow models, and may be useful to model a flotation cell within plant-wide simulators [2]. Nevertheless, for practical use, the present model should be improved and refined in future work. Some suggested directions of future work are as follows. Starting with a property of the model in its present formulation, we recall that the present analysis is limited to drift flux functions j g that satisfy (2.10), i.e., that have a horizontal tangent at φ = 1, for which the Richardson-Zaki expression (2.3), (2.4) is the most common example. Results of the steady-state analysis will be more complicated, but can be obtained by similar techniques to those of Section 3, wheneverq = −j g (1) > 0. To underline that this case is relevant for investigation we mention that Pal and Masliyah [31] come to the conclusion that the Richardson-Zaki expression (2.4), withq = 0, is suitable with n RZ = 2.39. However, their Figures 10 and 11, displaying drift flux data obtained by their own experiments and from the literature, including those obtained within the froth region, strongly indicate that a function j g withq > 0 is more suitable. In fact, in their next paper [32] they modify the expression for V to V (φ) = 0.8(1 − φ) exp(−2.9φ 2.1 ), such that indeed −j g (1) = 0.8 exp(−2.9) > 0.
With respect to the hydrodynamical setup of the flotation column and its conceptual counterpart as drawn in Figure 1, we mention that Vandenberghe et al. [39] propose an interesting recirculation: according to their Figure 1, mixture is sucked from the column at a determined level, aerated, and re-injected at another position. The extraction of material at a given rate but whose composition is part of the solution gives rise to a singular sink term whose mathematical treatment is more involved than that of a singular source term (as those appearing in (1.1)). The basic difficulty is that the sink term cannot be incorporated into the flux function; rather, the sink is represented by a new non-conservative transport term (see [6]).
In several instances, our analyses invoke available mathematical and numerical results for clarifier-thickener models. If sediment compressibility is included in a clarifier-thickener model, an effect that arises if the solid-liquid suspension is flocculated, then the governing equation for such models features an additional strongly degenerating diffusion term [10] whose appropriate treatment, roughly speaking, arises from handling it as part of the convective flux. Such a term may also be motivated in the application to flotation in future works: for instance, Narsimhan [29] derives such a term to account for gradual compressibility of the foam layer that is caused by flow of liquid through a network of plateau borders due to gravitational and capillary forces [29]. On the other hand, Stevenson et al. [37] propose a convection-diffusion model for the transport of gangue in flotation froth.
Clearly, an obvious shortcoming of the present description, although it is in line with the cited treatments of literature [13,22,25,27,31,32,38,39,43,44], is that it does not explicitly model the transport and settling of solid particles. It would be highly desirable to extend the model by solids phases, for example of hydrophobic and hydrophilic particles (of minerals and gangue material), and to include their attachment to bubble and transport via the liquid and gas components. The likely outcome of such a description is a convection-diffusion-reaction system with discontinuous flux akin to a recently advanced model of continuous sedimentation with reactions [4]. On the other hand, several gas or solid phases representing size classes that segregate and form areas of different composition can be included, and lead to first-order hyperbolic systems of conservation laws with discontinuous flux, under determined circumstances (see, e.g., [11] and the references cited in that paper).
Within the present work the emphasis has been on the construction of stationary solutions. The numerical scheme utilized, Godunov's scheme with suitable modifications to handle the flux discontinuities (see Section 4), is monotone provided that ∆t and ∆z satisfy the CFL condition (4.4), and therefore delivers approximate solutions that converge to the unique entropy solution of the problem (1.1), (1.2). However, formally second-order accurate solvers are possible for instance by techniques of variable extrapolation, see e.g. [7,9].