Tewfik Sari
email: tewfik.sari@uha.fr

Frederic Mazenc
email: mazenc@supagro.inra.fr

GLOBAL DYNAMICS OF THE CHEMOSTAT WITH DIFFERENT REMOVAL RATES AND VARIABLE YIELDS

Keywords: chemostat, competitive exclusion principle, Lyapunov function, global asymptotic stability, variable yield model AMS subject classifications. 92A15, 92A17, 34C15, 34C35

In this paper, we consider a competition model between n species in a chemostat including both monotone and non-monotone growth functions, distinct removal rates and variable yields. We show that only the species with the lowest break-even concentration survives, provided that additional technical conditions on the growth functions and yields are satisfied. We construct a Lyapunov function which reduces to the Lyapunov function used by S. B. Hsu [SIAM J. Appl. Math., 34 (1978), pp. 760-763] in the Monod case when the growth functions are of Michaelis-Menten type and the yields are constant. Various applications are given including linear, quadratic and cubic yields.

Introduction.

In this paper we study the global dynamics of the following model of the chemostat in which n populations of microorganisms compete for a single growth-limiting substrate:

S ′ (t) = D[S 0 -S(t)] - n i=1 f i (S(t))x i (t) x ′ i (t) = [p i (S(t)) -D i ]x i (t), i = 1 • • • n, (1.1)
where S(0) ≥ 0 and x i (0) > 0, i = 1 • • • n and S 0 , D and D i are positive constants. In these equations, S(t) denotes the concentration of the substrate at time t; x i (t) denotes the concentration of the ith population of microorganisms at time t; f i (S) represents the uptake rate of substrate of the ith population; p i (S) represents the per-capita growth rate of the ith population and so the function y i (S), defined by y i (S) = pi (S) fi(S) is the growth yield; S 0 and D denote, respectively, the concentration of substrate in the feed bottle and the flow rate of the chemostat; each D i represents the removal rate of the ith population. For general background on model (1.1), in the constant yield case y i (S) = Y i , the reader is referred to the monograph of Smith and Waltman [START_REF] Smith | The Theory of the Chemostat, Dynamics of Microbial Competition[END_REF].

The global analysis of this model was considered by Hsu, Hubbell and Waltman [START_REF] Hsu | A mathematical theory for single nutrient competition in continuous culture of micro-organisms[END_REF], in the Monod case [START_REF] Monod | La technique de culture continue. Théorie et applications[END_REF] when the growth functions are of Michaelis-Menten form,

p i (S) = a i S b i + S , (1.2)
and the yields are constant y i (S) = Y i , and D i = D for i = 1 • • • n. The authors showed that only the species with the lowest break-even concentration survives. Thus the competitive exclusion principle (CEP) holds: only one species survives, namely the species which makes optimal use of the resources. Hsu [START_REF] Hsu | Limiting behavior for competing species[END_REF] applied a Lyapunov argument to give a simple and elegant proof of the result in [START_REF] Hsu | A mathematical theory for single nutrient competition in continuous culture of micro-organisms[END_REF] for the case of different removal rates D i . The Lyapunov function V H used by Hsu is

V H = S λ1 σ -λ 1 σ dσ + c 1 x1 x * 1 ξ -x * 1 ξ dξ + n i=2 c i x i , (1.3) where c i = 1 Yi ai ai-Di , i = 1 • • • n, x * 1 = DY 1 S 0 -λ1 D1
and λ 1 = b1D1 a1-D1 is the lowest break-even concentration of the species.

Wolkowicz and Lu [START_REF] Wolkowicz | Global dynamics of a mathematical model of competition in the chemostat: general response functions and differential death rates[END_REF] extended the results of [START_REF] Hsu | Limiting behavior for competing species[END_REF] by allowing more general growth functions. These authors used the Lyapunov function

V W L = S 0 -λ 1 D 1 S λ1 p 1 (σ) -D 1 S 0 -σ dσ + 1 Y 1 x1 x * 1 ξ -x * 1 ξ dξ + n i=2 α i Y i x i , (1.4)
where α i , i = 2 • • • n are positive constants to be determined. They identified a large class of growth functions, including many prototypes of growth functions often found in the literature, where the constant α i in (1.4) can always be found. Despite the fact the α i cannot be found for all growth functions, the work of Wolkowicz and Lu [START_REF] Wolkowicz | Global dynamics of a mathematical model of competition in the chemostat: general response functions and differential death rates[END_REF] represents a major step in the extension of the result of Hsu [START_REF] Hsu | Limiting behavior for competing species[END_REF] to general growth functions.

The CEP has also been proved under a variety of hypotheses by Armstrong and McGehee [START_REF] Armstrong | Competitive exclusion[END_REF], Butler and Wolkowicz [START_REF] Butler | A mathematical model of the chemostat with a general class of functions describing nutrient uptake[END_REF], Wolkowicz and Xia [START_REF] Wolkowicz | Global asymptotic behavior of a chemostat model with discrete delays[END_REF] and Li [START_REF] Li | Global asymptotic behavior of the chemostat : general response functions and differential removal rates[END_REF]. The hypotheses used in [START_REF] Armstrong | Competitive exclusion[END_REF][START_REF] Butler | A mathematical model of the chemostat with a general class of functions describing nutrient uptake[END_REF][START_REF] Hsu | A mathematical theory for single nutrient competition in continuous culture of micro-organisms[END_REF][START_REF] Hsu | Limiting behavior for competing species[END_REF][START_REF] Li | Global asymptotic behavior of the chemostat : general response functions and differential removal rates[END_REF][START_REF] Wolkowicz | Global dynamics of a mathematical model of competition in the chemostat: general response functions and differential death rates[END_REF][START_REF] Wolkowicz | Global asymptotic behavior of a chemostat model with discrete delays[END_REF] are summarized in Table 1 of [START_REF] De Leenheer | Competition in the chemostat: some remarks[END_REF]. However, the problem is not yet completely solved: the CEP holds for a large class of growth functions but an important open question remains: is the CEP true assuming only that the f i are monotone with no restriction on the D i ? This major open problem remains unresolved, see in particular [START_REF] De Leenheer | Competition in the chemostat: some remarks[END_REF][START_REF] Wolkowicz | Microbial dynamics in a chemostat : competition, growth, implication of enrichment[END_REF]. For other studies and complements on the use of Lyapunov techniques in the chemostat, see [START_REF] Gajardo | Competitive exclusion principle in a model of chemostat with delays[END_REF][START_REF] Lobry | Effect on persistence of intra-specific competition in competition models[END_REF][START_REF] Mazenc | Harmand Stabilization in a two-species chemostat with Monod growth functions[END_REF][START_REF] Mazenc | Harmand Further results on stabilization of periodic trajectories for a chemostat with two species[END_REF][START_REF] Mazenc | On the stability of periodic solutions in the perturbed chemostat[END_REF].

The variable yield case was considered, for n = 1, 2 by Pilyugin and Waltman [START_REF] Pilyugin | Multiple limit cycles in the chemostat with variable yields[END_REF], with a particular interest to linear and quadratic yields, and by Huang, Zhu and Chang [START_REF] Huang | Limit cycles in a chemostat with variable yields and growth rates[END_REF]. The model (1.1), with variable yields, was considered by Arino, Pilyugin and Wolkowicz [START_REF] Arino | Considerations on yield, nutrient uptake, cellular growth, and competition in chemostat models[END_REF]. For biological motivations concerning the dependence of the yields on the substrate, see [START_REF] Arino | Considerations on yield, nutrient uptake, cellular growth, and competition in chemostat models[END_REF][START_REF] Pilyugin | Multiple limit cycles in the chemostat with variable yields[END_REF] and the references therein.

Notice that, in the case when the growth functions are of Michaelis-Menten form (1.2), the Lyapunov function (1.4) does not reduce to the Lyapunov function (1.3). Our aim in this paper is to extend the Lyapunov function (1.3) of Hsu [START_REF] Hsu | Limiting behavior for competing species[END_REF] to the chemostat with a more general class of growth functions and variable yields. Our Lyapunov function is given by

V = S λ1 p 1 (σ) -D 1 f 1 (σ) dσ + x1 x * 1 ξ -x * 1 ξ dξ + n i=2 α i x i (1.5)
where α i , i = 2 • • • n are positive constants to be determined. This Lyapunov function is just a multiple of the Lyapunov function (1.3) that Hsu used in [START_REF] Hsu | Limiting behavior for competing species[END_REF] in the Monod case, see Section 3.1. It is also a multiple of the one used in [START_REF] Wolkowicz | Global asymptotic behavior of a chemostat model with discrete delays[END_REF] Page 1039 or [START_REF] Wolkowicz | Microbial dynamics in a chemostat : competition, growth, implication of enrichment[END_REF] Section 3.3, in the case of one species, general growth function and constant yield, see Section 3.2.

The paper is organized as follows. In section 2 we prove our main result (see Theorem 2.2) and we compare it with the main result in [START_REF] Wolkowicz | Global dynamics of a mathematical model of competition in the chemostat: general response functions and differential death rates[END_REF] (see Theorem 2.3), where the yields are assumed to be constant. It should be noticed that, in the case when the yields are constant, our result follows from the result in [START_REF] Wolkowicz | Global dynamics of a mathematical model of competition in the chemostat: general response functions and differential death rates[END_REF]. Actually, both theorems 2.2 and 2.3 are corollaries of a more general result, which is valid in the case when the yields are variable [START_REF] Sari | Global dynamics of the chemostat with variable yields[END_REF]. In Section 3.1 we consider the Monod model with constant yields. In Section 3.2 we consider the one species case and we show that our Lyapunov function can be used to obtain the same result as in [START_REF] Arino | Considerations on yield, nutrient uptake, cellular growth, and competition in chemostat models[END_REF]. In Section 3.3 we show that for the Monod model with constant yields replaced by either linear or quadratic functions of S, under certain additional technical assumptions, the CEP still holds (see Corollary 3.1). In Section 3.4 we consider the model of Pilyugin and Waltman [START_REF] Pilyugin | Multiple limit cycles in the chemostat with variable yields[END_REF] which was used to demonstrate that a periodic orbit was possible in the case of variable yield model. In this model, with two species, where one yield is constant and the other is cubic in S, we show that our Lyapunov function can be used to prove that for some values of the parameters the CEP holds (see Corollary 3.2). In Section 3.5 we identify a class of growth functions, including Lotka-Volterra and Michaelis-Menten growth functions where our Lyapunov function works. Concluding remarks are given in Section 4.

2. Global asymptotic stability. We make the following assumptions on the functions p i and f i :

• p i , f i : R + → R + are continuous,
• p i (0) = f i (0) = 0 and for all S > 0, p i (S) > 0 and f i (S) > 0. Following Butler and Wolkowicz [START_REF] Butler | A mathematical model of the chemostat with a general class of functions describing nutrient uptake[END_REF], we make the following assumptions on the form of the growth functions p i : there exist positive extended real numbers λ i and µ i with

λ i ≤ µ i ≤ +∞ such that p i (S) < D i if S / ∈ [λ i , µ i ],
and

p i (S) > D i if S ∈]λ i , µ i [.
Hence there are at most two values, S = λ i and S = µ i , called the break-even concentrations, satisfying the equation p i (S) = D i . We adopt the convention µ i = ∞ if this equation has only one solution and λ i = ∞ if it has no solution.

It is known (see Theorem 4.1 [START_REF] Arino | Considerations on yield, nutrient uptake, cellular growth, and competition in chemostat models[END_REF]) that the non-negative cone is invariant under the flow of (1.1) and all solutions are defined and remain bounded for all t ≥ 0. System (1.1) can have many equilibria: the washout equilibrium E 0 = (S 0 , 0, • • • , 0), which is locally exponentially stable if and only if for all i = 1 • • • n, S 0 / ∈ [λ i , µ i ] and the equilibria E * i and E * * i where all components of E * i and E * * i vanish except for the first and the (i + 1)th, which are

S = λ i , x i = x * i = F i (λ i ), for E * i and S = µ i , x i = x * * i = F i (µ i ), for E * * i respectively, where F i (S) = D S 0 -S f i (S) .
The equilibrium E * i lies in the non-negative cone if and only if λ i ≤ S 0 . If λ i < λ j for all i = j and F ′ i (λ i ) < 0 then it is locally exponentially stable. It coalesces with E 0 when λ i = S 0 . The equilibrium E * * i lies in the non-negative cone if and only if µ i ≤ S 0 and is locally exponentially unstable if it exists. Its coalesces with E 0 when µ i = S 0 . Besides these equilibria, the system (1.1) can have a continuous set of non-isolated equilibria in the non-generic cases where two or more of the break-even concentrations are equal. In what follows we assume, that λ 1 < λ 2 ≤ • • • ≤ λ n , and λ 1 < S 0 < µ 1 . Hence E 0 is locally exponentially unstable and the equilibrium

E * 1 = (λ 1 , x * 1 , 0, • • • , 0), where x * 1 = F 1 (λ 1 ) = D S 0 -λ1 f1(λ1)
, lies in the non-negative cone. It is locally exponentially stable if and only if F ′ 1 (λ 1 ) < 0. We consider the global asymptotic stability of E * 1 .

Before presenting the results, we need the following lemma, Lemma 2.1. The solutions S(t), x i (t), i = 1 • • • n of (1.1) with positive initial conditions are positive and bounded, and if λ i < S 0 < µ i for some i = 1 • • • n, then S(t) < S 0 for all sufficiently large t.

Proof. The proof is similar to the proof of Lemma 2.1 in [START_REF] Wolkowicz | Global dynamics of a mathematical model of competition in the chemostat: general response functions and differential death rates[END_REF] obtained for the model (1.1) in the case where the yields are constant. We have the following result. Theorem 2.2. Assume that 1.

λ 1 < λ 2 ≤ • • • ≤ λ n , and λ 1 < S 0 < µ 1 . 2. There exist constants α i > 0 for each i ≥ 2 satisfying λ i < S 0 , such that max 0<S<λ1 g i (S) ≤ α i ≤ min λi<S<ρi g i (S), (2.1)
where g i (S) = fi(S) f1(S) p1(S)-D1 pi(S)-Di and ρ i = min(µ i , S 0 ). 3. The function F 1 (S) = D S 0 -S f1(S) satisfies F 1 (S) > F 1 (λ 1 ) if S ∈]0, λ 1 [, and
F 1 (S) < F 1 (λ 1 ) if S ∈]λ 1 , S 0 [. Then the equilibrium E *
1 is globally asymptotically stable for system (1.1) with respect to the interior of the positive cone.

Proof. From Lemma 2.1 it follows that there is no loss of generality in restricting our attention to 0 ≤ S < S 0 . Consider the function V = V (S, x 1 , • • • , x n ) given by (1.5), where α i are the positive constants satisfying (2.1). The function V is continuously differentiable in the positive cone and positive except at the point E * 1 , where it is equal to 0. The derivative of V along the trajectories of (1.1) is given by

V ′ = [p 1 (S) -D 1 ] [F 1 (S) -x * 1 ] + n i=2 x i θ i (S),
where

θ i (S) = [p i (S) -D i ] [α i -g i (S)]
. First, note that, using hypotheses 1 and 3, the first term of the above sum is always non-positive for 0 < S < S 0 and equals 0 for

S ∈]0, S 0 [ if and only if S = λ 1 . If S ∈ [λ 1 , λ i ] then p i (S) < D i and p 1 (S) > D 1 so that g i (S) < 0 < α i for any choice of α i > 0. Similarly if µ i < S 0 and S ∈ [µ i , S 0 ]
then p i (S) < D i and p 1 (S) > D 1 so that g i (S) < 0 < α i for any choice of α i > 0.

On the other hand, if S ∈ [0, λ 1 ] then p i (S) < D i and, using (2.1), g i (S) ≤ α i so that θ i (S) < 0. Finally, if S ∈ [λ i , ρ i ] then p i (S) > D i and g i (S) ≥ α i so that θ i (S) < 0. Thus θ i (S) < 0 for every S ∈]0, S 0 [, provided that the numbers α i satisfy (2.1). Hence V ′ ≤ 0 and V ′ = 0 if and only if S = λ 1 and x i = 0 for i = 2 • • • n. By the Krasovskii-LaSalle extension Theorem, the ω-limit set of the trajectory is E * 1 . In the case when the yields are constant, y i (S) = Y i , (1.1) takes the form

S ′ = D[S 0 -S] - n i=1 p i (S) Y i x i , x ′ i = [p i (S) -D i ]x i , i = 1 • • • n. (2.2)
Using the Lyapunov function (1.4), Wolkowicz and Lu (see Theorem 2.3 in [START_REF] Wolkowicz | Global dynamics of a mathematical model of competition in the chemostat: general response functions and differential death rates[END_REF]) proved the following result Theorem 2.3. Assume that 1.

λ 1 < λ 2 ≤ • • • ≤ λ n , and λ 1 < S 0 < µ 1 . 2. There exist constants α W L i > 0 for each i ≥ 2 satisfying λ i < S 0 such that max 0<S<λ1 g W L i (S) ≤ α W L i ≤ min λi<S<ρi g W L i (S), (2.3) where g W L i (S) = pi(S) D1 p1(S)-D1 pi(S)-Di S 0 -λ1
S 0 -S and ρ i = min(µ i , S 0 ). Then the equilibrium E * 1 is globally asymptotically stable for system (2.2) with respect to the interior of the positive cone. It should be noticed that, in the case when the yields are constant, conditions (2.3) are consequences of conditions (2.1) in Theorem 2.2. Indeed, we have

g W L i (S) = (S 0 -λ 1 )Y i D 1 f 1 (S) S 0 -S g i (S) = (S 0 -λ 1 )Y i D D 1 g i (S) F 1 (S) .
Thus, hypotheses 2 and 3 of Theorem 2.2 imply hypothesis 2 of Theorem 2.3. Hence, in the case when the yields are constant, Theorem 2.2 follows from Theorem 2.3. It is of interest to identify classes of growth functions where conditions (2.1) are satisfied, and hence Theorem 2.2 can be applied. We give below a result which will be used in the following section to verify easily that conditions (2.1) are satisfied. This proposition is similar to Corollary 2.4 in [START_REF] Wolkowicz | Global dynamics of a mathematical model of competition in the chemostat: general response functions and differential death rates[END_REF].

Proposition 2.4. Suppose that for each i ≥ 2, w i (S) = g i (S) S-λi S-λ1 satisfies max 0<S<λ1 w i (S) ≤ min λi<S<ρi w i (S). (2.4) Then conditions (2.1) are satisfied.

Proof. The function h i (S) = S-λ1 S-λi is decreasing on [0, λ 1 ] and on ]λ i , +∞[ and admits 1 as a horizontal asymptote. Thus (see Figure 3.1, right)

max 0<S<λ1 h i (S) = h i (0) < 1 < h i (ρ i ) = min λi<S<ρi h i (S). (2.5)
By (2.4) and (2.5) the functions g i (S) = w i (S)h i (S) satisfy (2.1).

Applications.

In this section we show how Theorem 2.2 can be fruitfully used to analyze the stability properties of systems whose yield functions depend on the variable S. We begin with the classical Monod case where the yields are constant and the growth functions are of Michaelis-Menten form.

3.1. The Monod case. Consider the particular case where the growth functions p i (S) are given by (1.2) and the yields are constant. System (1.1) takes the form

S ′ = D(S 0 -S) - n i=1 a i S b i + S x i Y i , x ′ i = a i S b i + S -D i x i , i = 1 • • • n. (3.1)
We consider the case where, for all i = 1 • • • n, a i > D i . In that case:

λ i = b i D i a i -D i , µ i = ∞, g i (S) = w i S -λ 1 S -λ i , where w i = a i Y 1 (a 1 -D 1 ) a 1 Y i (a i -D i ) .
By Proposition 2.4, the conditions (2.1) are satisfied. Since

F 1 (S) = Y 1 D(S 0 -S) b 1 + S a 1 S and F ′ 1 (S) = -Y 1 D S 2 + b 1 S 0 a 1 S 2 ,
the first derivative of the function F 1 (S) is negative. Hence, hypothesis 3 in Theorem 2.2 is satisfied. The global stability of the equilibrium E * 1 of (3.1) follows from Theorem 2.2. This result was obtained by Hsu [START_REF] Hsu | Limiting behavior for competing species[END_REF], using the Lyapunov function (1.3). Notice that, in this case, the Lyapunov function (1.5) is simply V = a1-D1 a1Y1 V H where V H is the Lyapunov function (1.3) used by Hsu [START_REF] Hsu | Limiting behavior for competing species[END_REF].

S gi λ1 λi Ai λ 1 λ i Ai c 1 c i S hi λ1 λi λ 1 λ i 1 Fig. 3.1.
On the left, the graph of the function g i (S) for p i (S) = a i S b i +S and y i (S) = Y i (1 + c i S): hypothesis 2 in Corollary 3.1 is not satisfied. On the right the graph of the function h i (S) = S-λ 1 S-λ i .

3.2. One species. In the case n = 1, (1.1) takes the form

S ′ = D(S 0 -S) -x 1 f 1 (S), x ′ 1 = [p 1 (S) -D 1 ]x 1 (3.2)
If λ 1 < S 0 < µ 1 and hypothesis 3 in Theorem 2.2 is satisfied then the equilibrium E * 1 = (λ 1 , x * 1 ) of (3.2), where x * 1 = D S 0 -λ1 f1(λ1) is globally asymptotically stable with respect to the interior of the positive quadrant. This result follows from Theorem 2.2 since in the case where n = 1 the condition (2.1) is obviously satisfied. The global asymptotic stability of E * 1 was obtained previously by Arino, Pilyugin and Wolkowicz [START_REF] Arino | Considerations on yield, nutrient uptake, cellular growth, and competition in chemostat models[END_REF]. These authors used the following Lyapunov function

V AP W = S 0 -λ 1 f 1 (λ 1 ) S λ1 p 1 (σ) -D 1 S 0 -σ dσ + x1 x * 1 ξ -x * 1 ξ dξ.
They proved (see [START_REF] Arino | Considerations on yield, nutrient uptake, cellular growth, and competition in chemostat models[END_REF], Theorem 2.11) that if 1 -f1(S)(S 0 -λ1) f1(λ1)(S 0 -S) has exactly one sign change for S ∈ (0, S 0 ) then E * 1 is globally asymptotically stable. The condition on the change of sign is equivalent to hypothesis 3 in Theorem 2.2. Notice that the Lyapunov function we obtain is not proportional to the Lyapunov function V AP W considered in [START_REF] Arino | Considerations on yield, nutrient uptake, cellular growth, and competition in chemostat models[END_REF]. However, in the case when the yields is constant the global asymptotic stability of the equilibrium E * 1 of the system

S ′ = D(S 0 -S) -x 1 p 1 (S) Y 1 , x ′ 1 = [p 1 (S) -D 1 ]x 1 ,
can be obtained, using the following Lyapunov function (see [START_REF] Wolkowicz | Global asymptotic behavior of a chemostat model with discrete delays[END_REF] page 1039 or [START_REF] Wolkowicz | Microbial dynamics in a chemostat : competition, growth, implication of enrichment[END_REF], Section 3.3)

V W BL = S λ1 p 1 (σ) -D 1 p 1 (σ) dσ + 1 Y 1 x1 x * 1 ξ -x * 1 ξ dξ.
In this case we simply have V = Y 1 V W BL , where V is our Lyapunov function (1.5).

3.3. Michaelis-Menten growth functions and linear or quadratic yields. Consider the particular case of (1.1), where the growth functions p i (S) are given by (1.2), and the yields y i (S) = p i (S)/f i (S) are linear

y i (S) = Y i (1 + c i S) (3.3) or quadratic y i (S) = Y i (1 + c i S 2 ). (3.4)
where Y i > 0 and c i ≥ 0. System (1.1) takes the form

S ′ = D(S 0 -S) - n i=1 a i S b i + S x i y i (S) , x ′ i = a i S b i + S -D i x i , i = 1 • • • n. (3.5)
Corollary 3.1. Consider system (3.5) where the yields are given by (3.3) or (3.4). Assume that 1.

λ 1 < λ 2 ≤ • • • ≤ λ n and λ 1 < S 0 , where λ i = biDi ai-Di . 2. For each i ≥ 2 satisfying λ i < S 0 we have c i λ 1 ≤ c 1 λ i . 3. The function F 1 (S) = D S 0 -S a1S (b 1 + S)y 1 (S) satisfies F 1 (S) > F 1 (λ 1 ) if S ∈ ]0, λ 1 [, and F 1 (S) < F 1 (λ 1 ) if S ∈]λ 1 , S 0 [. Then the equilibrium E *
1 is globally asymptotically stable for (3.5) with respect to the interior of the positive cone.

Proof. For linear yields (3.3) we have

g i (S) = w i (S) S -λ 1 S -λ i , where w i (S) = A i 1 + c 1 S 1 + c i S , A i = a i Y 1 (a 1 -D 1 ) a 1 Y i (a i -D i ) .
Two cases can be distinguished. If By Proposition 2.4 the functions g i (S) = w i (S)h i (S) satisfy (2.1). If c 1 < c i then from the expression

g ′ i (S) = A i c 1 -c i (1 + c i S) 2 S -λ 1 S -λ i + A i 1 + c 1 S 1 + c i S λ 1 -λ i (S -λ i ) 2 ,
we deduce that g ′ i (S) < 0 for all 0 ≤ S ≤ λ 1 and S > λ i . Hence (see Fig. 3.1, left) max

0<S<λ1 g i (S) = g i (0) = A i λ 1 λ i and min λi<S<S 0 g i (S) = g i (S 0 ) > g i (+∞) = A i c 1 c i .
Under hypothesis 2 there exists α i satisfying (2.1). The result follows by Theorem 2.2. For quadratic yields (3.4) we have

g i (S) = w i (S) S -λ 1 S -λ i , where w i (S) = A i 1 + c 1 S 2 1 + c i S 2 , A i = a i Y 1 (a 1 -D 1 ) a 1 Y i (a i -D i ) .
Thus

g ′ i (S) = A i 2(c 1 -c i )S (1 + c i S 2 ) 2 S -λ 1 S -λ i + A i 1 + c 1 S 2 1 + c i S 2 λ 1 -λ i (S -λ i ) 2 .
Next, the proof is mutatis mutandis the same as the proof given above for the case of linear yields (3.3). This result contains as a particular case the result of Hsu [START_REF] Hsu | Limiting behavior for competing species[END_REF] which corresponds to the case where the yields are constant. Indeed, for constant yields c i = 0, so that hypotheses 2 and 3 in Corollary 3.1 are satisfied. Remark. For linear or quadratic yields the function F 1 (S) is not monotone in general on the interval ]0, S 0 [, and it is not easy to give a condition on the parameters for which hypothesis 3 in Corollary 3.1 holds. However, in each example, the graphical depiction of this hypothesis is very simple as shown in Fig. 3.2.

3.4. Pilyugin-Waltman's example. This system was given in [START_REF] Pilyugin | Multiple limit cycles in the chemostat with variable yields[END_REF] as a model of the competition in the chemostat exhibiting limit cycles. The existence of the limit cycles is a consequence of the variable yield in the model. The model takes the form

S ′ = 1 -S -2S 0.7+S x1 1+cS 3 -m2S 6.5+S x2 120 x ′ 1 = [ 2S 0.7+S -1]x 1 x ′ 2 = [ m2S 6.5+S -1]x 2 .
(3.6)

In their study Pilyugin and Waltman [START_REF] Pilyugin | Multiple limit cycles in the chemostat with variable yields[END_REF] fixed c = 50 and considered m 2 as a bifurcation parameter. They showed that for m 2 ≥ 9.85 the system exhibits sustained oscillations. In this section we fix m 2 = 10 and we consider c ≥ 0 as a bifurcation parameter. In this case we have

λ 1 = 0.7 2 -1 = 0.7, λ 2 = 6.5 10 -1 ≈ 0.72, F 1 (S) = (1 -S)(0.7 + S)(1 + cS 3 ) 2S .
Straightforward computations lead to the formula

F ′ 1 (λ 1 ) = 49 2000 c - 17 14 .
Hence F ′ 1 (λ 1 ) > 0 if and only if c > c 2 where c 2 = 17000 343 . An analysis of the behavior of the function F 1 (S) shows (see Fig. 3.3) that there exist two bifurcation values c 0 and c 1 , 0 < c 0 < c 1 < c 2 , such that the function F 1 (S) is decreasing on ]0, S 0 [ if and only if 0 ≤ c ≤ c 0 and the function

S 0 0, 2 0,4 0,6 0,8 1,0 0 0,5 1,0 1,5 2,0 S 0 0,2 0,4 0,6 0,8 1,0 0 0,5 1,0 1,5 2,0 S
F 1 (S) has two extrema S 1 , S 2 ∈]0, λ 1 [ satisfying F 1 (S 1 ) < F 1 (λ 1 ) < F 1 (S 2 ) if and only if c ∈]c 1 , c 2 [. S 0 0,2 0,4 0,6 0, 8 1,0 0 0,5 1,0 1,5 2,0 S 0
g 2 (S) = w 2 (S) S -λ 1 S -λ 2 where w 2 (S) = 1 + cS 3 216 .
For c ≥ 0, the function w 2 (S) is non-decreasing. By Proposition 2.4, the condition (2.1) with i = 2 holds (see Fig. 3.4, right), and the result follows from Theorem 2.2.

Pilyugin and Waltman showed by numerical simulations that their system exhibits limit cycles in the case where c = 50 and m 2 ≥ 9.85 (see Fig. 4 in [START_REF] Pilyugin | Multiple limit cycles in the chemostat with variable yields[END_REF]). The example was revisited by Huang, Zhu and Chang [START_REF] Huang | Limit cycles in a chemostat with variable yields and growth rates[END_REF] who claimed that the limit cycle of the system should remain only on the face x 2 = 0 (see [START_REF] Huang | Limit cycles in a chemostat with variable yields and growth rates[END_REF], Remark 2). We do not agree with this claim. We performed ourselves numerical simulations and actually the limit cycle is contained within the positive cone as shown in Fig. 4 in [START_REF] Pilyugin | Multiple limit cycles in the chemostat with variable yields[END_REF] and not in the face x 2 = 0 as claimed in [START_REF] Huang | Limit cycles in a chemostat with variable yields and growth rates[END_REF]. On the center, the magnification of the neighborhood of λ 1 = 0.7 shows that F ′ 1 (λ 1 ) > 0. On the right, the magnification of the neighborhood of λ 1 = 0.71 shows that F ′ 1 (λ 1 ) < 0.

Huang, Zhu and Chang [START_REF] Huang | Limit cycles in a chemostat with variable yields and growth rates[END_REF] made a simple modification by replacing 2S/(0.7 + S) with 2S/(0.71 + S) in (3.6) and obtained an example exhibiting competitive exclusion. The model takes the form

S ′ = 1 -S -2S 0.71+S x1 1+50S 3 -m2S 6.5+S x2 120 x ′ 1 = [ 2S 0.7+S -1]x 1 x ′ 2 = [ m2S 6.5+S -1]x 2 . (3.7)
It is claimed, without proof, in [START_REF] Huang | Limit cycles in a chemostat with variable yields and growth rates[END_REF] that the equilibrium E * 1 is globally asymptotically stable. Hypothesis 3 in Theorem 2.2 is not satisfied (see Fig. 3.5, left) and we cannot prove the global asymptotic stability of E * 1 . However an explanation of the high sensitivity when 0.7 is replaced by 0.71 is easy to find. Actually the plots of the function F 1 (S) in the case of (3.6), where c = 50 and (3.7) are very similar (see Fig. 3.5, left), but a magnification of the neighborhood of the value S = λ 1 shows the differences (see Fig. 3.5, center and right). In (3.6), F ′ 1 (λ 1 ) > 0. Hence the equilibrium E * 1 is locally exponentially unstable. In (3.7), F ′ 1 (λ 1 ) < 0. Hence the equilibrium E * 1 is locally exponentially stable.

3.5. Further applications. In this section we describe a class of growth functions p i (S) and yields y i (S) for which constants α i satisfying (2.1) exist and hence Theorem 2.2 can be applied. It is convenient to use the notation

P i (S) = S -λ i p i (S) -D i p i (S) =⇒ p i (S) = D i P i (S) P i (S) + λ i -S . (3.8)
Remark. We can take any functions P i (S) that are positive for 0 < S ≤ S 0 and satisfy P i (0) = 0 and use the righthand side of formulas (3.8) to define the functions p i (S). The function P i (S) must satisfy the condition P i (S) > S -λ i , so that le denominator in p i (S) remains positive. If we find a class of yield functions y i (S) such that the conditions (2.4) hold, where the functions w i (S), considered in Proposition 2.4, are given by w i (S) = y 1 (S) y i (S)

P i (S) P 1 (S)
then we can use Proposition 2.4 to obtain the global asymptotic stability of the equilibrium E * 1 . The Holling type II (Michaelis-Menten or Monod) growth functions p i (S) = D i m i S (m i -1)S + λ i (3.9) correspond to the choice P i (S) = m i and m i > 1. The Holling type III (or sigmoidal) growth functions

p i (S) = D i m i S 2 (a i + S)(b i + S) , with m i = (a i + λ i )(b i + λ i ) λ 2 i (3.10)
correspond to the choice P i (S) = (ai+λi)(bi+λi)S 2 (ai+bi)λiS+aibi(S+λi) . Here µ i = +∞. The prototype for a non-monotone growth function

p i (S) = D i m i S (a i + S)(b i + S) , with m i = (a i + λ i )(b i + λ i ) λ i (3.11) corresponds to the choice P i (S) = (ai+λi)(bi+λi) aibi-λiS
. Here µ i = aibi λi . The growth functions (3.9-3.11) were considered by Wolkowicz and Lu [START_REF] Wolkowicz | Global dynamics of a mathematical model of competition in the chemostat: general response functions and differential death rates[END_REF] who indicated for each combination of them that it is always possible to find appropriate constants α W L i satisfying the criterion (2.3).

Hereafter we define two new classes of functions, which are not considered in the literature, for which our results apply. A class of monotone growth functions of the form (3.8) is obtained with P i (S) = α i S + αS 1+βS , where α > 0, β > 0 and α i ≥ 1. In this case we have

p i (S) = D i α i S(1 + α + βS) αα i S + (1 + βS)(α i S -S + λ i ) .
For constant yields y i (S) = Y i the functions w i (S) = Y1αi Yiα1 are constant and hence conditions (2.4) are satisfied. For linear yields (3.3) we have w i (S) = Y1αi Yiα1 1+c1S 1+ciS . If c 1 ≥ c i then w i (S) is non-decreasing and hence conditions (2.4) are satisfied. For quadratic yields (3.4) we have w i (S) = Y1αi Yiα1 1+c1S 2 1+ciS 2 . If c 1 ≥ c i then w i (S) is nondecreasing and hence conditions (2.4) are satisfied.

A class of non-monotone growth functions of the form (3.8) is obtained with P i (S) = α i S 2 and α i > 1 4λi . In this case we have

p i (S) = D i α i S 2 α i S 2 -S + λ i .
For constant yields y i (S) = Y i the functions w i (S) = Y1αi Yiα1 are constant and hence conditions (2.4) are satisfied. For linear yields (3. 4. Discussion. In this paper we considered a mathematical model (1.1) of n species of microorganisms in competition in a chemostat for a single resource. The model incorporates both monotone and non-monotone growth functions, distinct removal rates and variable yields. We demonstrated that the CEP holds for a large class of growth functions and yields.

In the case where the yields are constant, it is known [START_REF] Butler | A mathematical model of the chemostat with a general class of functions describing nutrient uptake[END_REF] that the CEP holds provided that D i = D for all i, the set Q = i∈N ]λ i , µ i [ is connected, and S 0 ∈ Q, where N = {i : λ i < S 0 }. Wolkowicz and Lu [START_REF] Wolkowicz | Global dynamics of a mathematical model of competition in the chemostat: general response functions and differential death rates[END_REF] conjectured that this result can be extended to the case of different removal rates. Under hypothesis 1 in Theorem 2.2, it is clear that the set Q is connected, and S 0 ∈ Q. The condition λ 1 < λ i for i = 1 in hypothesis 1 can be stated without loss of generality, by labelling the populations such that the index i = 1 corresponds to the lowest break-even concentration, but the condition λ 1 < S 0 < µ 1 in hypothesis 1 cannot be stated without loss of generality. If µ 1 < S 0 , it is not possible to show the CEP by the methods that we used. To the best of our knowledge, in the case of different removal rates and non-monotone growth functions, the CEP has been proved only under the assumption S 0 < µ 1 [START_REF] Li | Global asymptotic behavior of the chemostat : general response functions and differential removal rates[END_REF][START_REF] Wolkowicz | Global dynamics of a mathematical model of competition in the chemostat: general response functions and differential death rates[END_REF][START_REF] Wolkowicz | Global asymptotic behavior of a chemostat model with discrete delays[END_REF]. However, Rapaport and Harmand [START_REF] Rapaport | Biological control of the chemostat with nonmonotone response and different removal rates[END_REF] considered the case of two populations and proposed conditions on the growth functions such that the CEP holds under the condition µ 1 < S 0 . It should be interesting to extend their methods to more general cases. We leave this problem for future investigations.

In the case of constant yields, numerical simulations of model (1.1) have only displayed competitive exclusion. Our results concern also the case of variable yields, for which it is known [START_REF] Arino | Considerations on yield, nutrient uptake, cellular growth, and competition in chemostat models[END_REF][START_REF] Huang | Limit cycles in a chemostat with variable yields and growth rates[END_REF][START_REF] Pilyugin | Multiple limit cycles in the chemostat with variable yields[END_REF]] that more exotic dynamical behaviors, including limit cycles and chaos, are possible. Thus in the case of variable yields, it is of great importance to have criteria ensuring the global convergence to an equilibrium with at most one surviving species. Under certain technical restrictions, we extended the result of Hsu [START_REF] Hsu | Limiting behavior for competing species[END_REF] to the case of linear or quadratic yields.

Our proof relies on the construction of non-strict Lyapunov functions, i.e. Lyapunov functions whose derivatives along the trajectories are non-positive. We conjecture that the strictification techniques of Chapter 5 of [START_REF] Malisoff | Constructions of strict Lyapunov functions[END_REF] can be used to construct strict Lyapunov functions, i.e. Lyapunov functions whose derivative along the trajectories are definite negative, which next can be used to establish some robustness properties. This can be the subject of further research.

Fig. 2 . 1 .

 21 Fig. 2.1. Graphical depictions of the hypotheses 2 and 3 in Theorem 2.2. On the left, hypothesis 2. On the right, hypothesis 3.

c 1 ≥

 1 c i then the function w i (S) is non-decreasing over [0, +∞[. Thus max 0<S<λ1 w i (S) ≤ min λi<S<ρi w i (S).

Fig. 3 . 2 .

 32 Fig. 3.2. The graph of the function F 1 (S) for p 1 (S) = a 1 S b 1 +S (where a 1 = 2.1/1.1, b 1 = 0.5/1.1 and λ 1 = 0.5) and y 1(S) = Y 1 (1 + c 1 S). On the left, the case c 1 = 1.8 for which hypothesis 3 in Corollary 3.1 is not satisfied. On the right, the case c 1 = 1 for which this hypothesis is satisfied.

cFig. 3 . 3 .

 33 Fig. 3.3. Graphs of the function F 1 (S) for p 1 (S) = 2S/(0.7 + S) and y 1 (S) = 1 + cS 3 in the cases c = c 0 , c 1 , c 2 .

Fig. 3 . 4 .

 34 Fig. 3.4. Graphical verification of hypotheses 2 and 3 in Theorem 2.2 for (3.6) with c = 8 and m 2 = 10. On the left, the graph of the function F 1 (S). On the right the graph of the function g 2 (S).

Fig. 3 . 5 .

 35 Fig.3.5. On the left, the plot of F 1 (S) for (3.6) or(3.7) where c = 50 and m 2 = 10. On the center, the magnification of the neighborhood of λ 1 = 0.7 shows that F ′ 1 (λ 1 ) > 0. On the right, the magnification of the neighborhood of λ 1 = 0.71 shows that F ′ 1 (λ 1 ) < 0.

2 1+ciS 2 .

 22 3) we have w i (S) = Y1αi Yiα1 1+c1S 1+ciS . If c 1 ≥ c i then w i (S) is non-decreasing and hence conditions (2.4) are satisfied. For quadratic yields (3.4) we have w i (S) = Y1αi Yiα1 1+c1S If c 1 ≥ c i then w i (S) is nondecreasing and hence conditions (2.4) are satisfied.

  Corollary 3.2. The equilibrium E * 1 is locally exponentially unstable if and only if c > c 2 . If 0 ≤ c < c 1 then the equilibrium E * 1 is globally asymptotically stable . Proof. Since λ 1 < λ 2 and F ′ 1 (λ 1 ) > 0 if and only if c > c 2 the equilibrium E * 1 is locally exponentially stable if and only if c > c 2 . Hypothesis 3 in Theorem 2.2 is satisfied if and only if 0 ≤ c < c 1 (see Fig. 3.4, left). The function g 2 (S) is defined by

Acknowledgments. We thank Alain Rapaport and Jérôme Harmand for stimulating discussions during the preparation of this work and Gail S.K. Wolkowicz for valuable comments and suggestions.


Florent Di Meglio

Glenn-Ole Kaasa

Nicolas Petit

Vidar Alstad

Reproducing slugging oscillations of a real oil well

This paper addresses the problem of reproducing oscillations generated by the well-known slugging phenomenon in multiphase flow. Reported investigations show how to determine the parameters of a recently proposed ordinary differential equations system, so that it captures the characteristics of actually observed slugging oscillations. A tuning procedure based on the mathematical properties of the model is presented. It is then applied to a test case consisting of a real oil well located in the North Sea. An observer using only topside measurements allows asymptotic reconstruction of critical variables such as the downhole pressure.

I. Introduction

For years, the slugging flow regime has been a concern for the oil industry. This phenomenon takes place in the long pipes connecting underground sources (reservoir or wellhead) to sub sea or surface production facilities. These pipes contain oil, water and gas forming a two-phase medium. Among all the possible two-phase flow regimes (bubbly, annular, churn, etc.) that can occur in this process of oil production, slugging is the less desirable one. In this regime, inhomogeneous spatial distribution of the liquid (usually oil and water) and gas phases coming from the sources generate oscillations of the flow rates and pressures inside the pipes, constituting an intermittent flow that is particularly visible at the outlet. These oscillations are referred to as the 'slugging' phenomenon, or slug flow. Their characteristics vary according to the pressure and temperature conditions, the geometry of the pipes or the composition of the gasliquid mixture. When their magnitude becomes too large, oil production must be stopped. A malicious effect is that, even when the production is not stopped, the oscillations reduce the overall oil production rate by increasing the average pressure in the pipes.

An important literature has been dedicated to the modeling of this phenomenon, in view of finding ways to suppress it. In the 1960s, hydrodynamic models [START_REF] Kordyban | A flow model for two-phase slug flow in horizontal tubes[END_REF], [START_REF] Fernandes | Hydrodynamic model for gas-liquid slug flow in vertical tubes[END_REF], [START_REF] Dukler | A model for gas-liquid slug flow in horizontal and near horizontal tubes[END_REF] were developed to describe in details the characteristics (speed, length, frequency) of the 'slugs' of liquid traveling along the pipes. These models were meant to help predict the behavior of certain systems, by formulating necessary conditions for the existence of the slugs. Such conditions are presented in [START_REF] Ruder | Necessary conditions for the existence of stable slugs[END_REF], and are based on the physical behavior of one slug1 . In the 1980s, transient mechanistic models have been proposed. They were very soon embedded into commercial multiphase flow simulators, such as TACITE [START_REF] Ferre | TACITE Physical Model Description Manual[END_REF] and OLGA [START_REF] Bendiksen | The dynamic twofluid model olga: Theory and application[END_REF]. These simulation-oriented models are composed of nonlinear coupled sets of partial differential equations (PDE). Their steady-state behavior is very accurate when compared against experimental or even real field data, for certain classes of systems (in terms of size, pressures, velocities, etc.). They also predict with a relatively good accuracy the transient behavior of certain slugging systems, especially risers with low-points2 . Yet, they are very difficult to analyze from a mathematical standpoint, besides trial anderror sensitivity analysis. The mere existence of solutions is by no means evident, let alone their stability.

More recently, simpler transient models have been developed for analytic studies and control design purposes. Jansen et al. describe in [START_REF] Jansen | Automatic control of unstable gas lifted wells[END_REF] the physical principles of a model further detailed in [START_REF] Imsland | Output Feedback and Stabilization and Control of Positive Systems[END_REF]. Their model, which is able to reproduce the slugging phenomenon for gas-lifted wells, is composed of three ordinary differential equations (ODE) with lumped states. They correspond to mass conversation laws for two volumes of gas, and one for the liquid. The relative simplicity of this model, compared to more involved PDE ones, allows for greater analysis possibilities. Sinègre [START_REF] Sinègre | Dynamic study of unstable phenomena stepping in gas-lift activated systems[END_REF] for example, has proven the existence of a periodic orbit for the model, under simplifying assumptions. More importantly, the model can be used to design control laws in view of suppressing the slugging phenomenon [START_REF] Havre | Active feedback control as the solution to severe slugging[END_REF], [START_REF] Imsland | Output Feedback and Stabilization and Control of Positive Systems[END_REF], or estimators for unmeasured pressures [START_REF] Eikrem | Stabilization of gas lifted wells[END_REF], [START_REF] Scibilia | Stabilization of gas-lift oil wells using topside measurements[END_REF], [START_REF] Sinègre | Dynamic study of unstable phenomena stepping in gas-lift activated systems[END_REF]. Similarly, Storkaas et al. [START_REF] Storkaas | Stabilization of severe slugging based on a low-dimensional nonlinear model[END_REF] have proposed a similar third order ODE model able to reproduce the behavior of slugging risers presenting a low-point. His study of the linear properties of the model gives insight into the controllability of the slugging phenomenon. Also, the model was used by Siahaan et al. [START_REF] Siahaan | Suppressing riser-based slugging in multiphase flow by state feedback[END_REF] to propose a nonlinear control structure.

The downside of having such ODE models is that they may be too simple to reflect all the subtle characteristics of the slugging phenomenon. As a result, some of the parameters of the model need to be carefully tuned, in order to allow the model to fit experimental data. Storkaas e.g. describes in [START_REF] Storkaas | Control solutions to avoid slug flow in pipeline-riser systems[END_REF] a tuning procedure used to find the values of five of the parameters of the model. This, of course, requires the a priori knowledge of data regarding the system under consideration. The upside is that, despite their relative simplicity, the models are able to compete with much more involved models, in terms of matching the behavior of real systems.

In this article, we show that the model presented in [START_REF] Di Meglio | A first principle model for multiphase slugging flow in vertical risers[END_REF] can reproduce with a very good accuracy the behavior of a real well. This is the contribution of this article. In a first step, the parameters of the model are estimated by simple tuning rules based on the physical properties and a mathematical analysis of the model. Then, an observer is designed to reconstruct unmeasured states from topside pressure measurements. These results must be put into perspective of the stabilization results presented in [START_REF] Di Meglio | Model-based control of slugging flow: an experimental case study[END_REF]. There, the proposed model-based controller was able to outperform PI controllers when tested on an experimental multiphase flow loop, which is a down-scaled version of a real system. The ability of the model to match the behavior of a real well, and to improve upon the performances of PI controllers on experimental facilities, suggest that a control design based on it has some potential to stabilize the flow on actual slugging wells.

The paper is organized as follows. In Section II, we review the physical principles and the equations of the model proposed in [START_REF] Di Meglio | A first principle model for multiphase slugging flow in vertical risers[END_REF]. Then, in Section III, we introduce the tuning rules and list the information needed. Then, in Section IV, we focus on the design of an observer. It is experimentally validated on real-well data. We give some conclusions and perspective in Section V.

II. Review of the model

Consider the pipe schematically depicted in Figure 1. It is filled with gas and liquid inflows from the reservoir at constant rates. The outflow can be regulated by a valve (production choke). The model proposed in [START_REF] Di Meglio | A first principle model for multiphase slugging flow in vertical risers[END_REF] stresses the possible existence of an obstructing liquid interface restricting the flow of gas, modeled as a virtual valve. In the case of a pipe presenting a low-point angle, this precisely corresponds to the idea originally proposed in the Storkaas model [START_REF] Storkaas | Stabilization of severe slugging based on a low-dimensional nonlinear model[END_REF], where the height of liquid at the low-point determines the opening of the virtual valve. The main modeling assumption in [START_REF] Di Meglio | A first principle model for multiphase slugging flow in vertical risers[END_REF] consists in considering the existence of such a valve, even when no low-point exists. The virtual valve separates the pipe into two parts. Upstream the valve, the existence of an elongated bubble of gas, of constant volume V eb , is assumed. The section downstream the valve is referred to as the riser. It has a length L, and the flow there is assumed to be distributed. Only a fraction of the gas (determined by 0 < < 1) accumulates in the elongated bubble, causing a build-up of pressure at the roots of the instability. When the pressure in the bubble gets high enough, gas flows through the virtual valve into the riser. This increases the pressure there, and pushes the liquid out of the system. In this model, the state variables are the mass

Production Choke

Virtual downhole valve Liquid Gas Elongated bubble Fig. 1. Schematic view of an inclined riser carrying an inhomogeneous blend of liquid and gas of gas in the elongated bubble m g,eb , the mass of gas in the riser m g,r , and the mass of liquid in the riser m l,r . The first principle equations read

ṁg,eb (t) = (1 -)w g,in -w g (t) (1) ṁg,r (t) = w g,in + w g (t) -w g,out (t) (2) ṁl,r (t) = w l,in -w l,out (t) (3)
where the g and l subscripts stand for gas and liquid respectively, in and out refer to inflow and outflow, and w g (t) is the flow of gas through the virtual valve. Apart from the inflows, which are assumed constant, the flow rates depend on the pressure drop over the choke under consideration. A linear relation of the form

w g = C g max(0, p eb -p r,b ) (4)
where p r,b is the pressure at the bottom of the riser, is chosen for the virtual valve. The max(0, •) function guarantees that no back flow can be admitted through the valve. The flow through the production choke, which is diphasic, requires a more accurate relation given by

w out = uC c √ p r,t -p s
where u ∈ (0, 1) is the controlled opening of the choke, p r,t and p s the pressure upstream and downstream the choke, respectively. The dependence with respect to density is neglected, assuming that the mixture density at the choke is constant. Then, the flows of gas and liquid through the production choke are given by their mass fractions, which yields

w l,out = m l,r m l,r + m g,r w out ≈ w out w g,out = m g,r m l,r + m g,r w out ≈ m g,r m l,r w out
Finally, the pressures (topside, and in the elongated bubble) are given by the ideal gas law, under isothermic assumptions.

The volume of the elongated bubble V eb is assumed to be constant, whereas the gas downstream the virtual choke is assumed to be compressible. Its volume V g,r depends on the mass of liquid in this part

V g,r = LA - m l,r + m l,still ρ l
where L is the length of the riser and A the cross-section area. m l,still represents the mass of liquid present at all times in the riser, which is not concerned by the dynamics. Besides, the pressure drop over the whole riser is assumed to be gravitydominated (friction being neglected)

p eb = RT MV eb m g,eb p r,b = p r,t + g sin θ A (m l,r + m l,still ) p r,t = RT M LA - m l,r +m l,still ρ l m g,r
where θ is the mean inclination of the pipe.

III. Tuning procedure

In this section, we describe how to choose the parameters in model ( 1)-( 2)-( 3) to match the behavior of a given system. These parameters are: L, θ, A, g, ρ l , R, T , M, p s , w g,in , w l,in , C c , m l,still , V eb and . The procedure requires a priori information on the system. First, the geometry of the pipe is required to set the values of L, A, and θ. The equilibrium values of the topside pressure ȳr,t (ū) and downhole pressure ȳr,b (ū) at one operating point ū are also required. These values can be measured, when corresponding sensors are available, or otherwise estimated by an advanced multiphase flow simulator (such as OLGA TM ). They determine the correct values of the choke constant C c and the still mass of liquid m l,still . Also, the value of the production choke opening that causes oscillations to appear on the real system u * is needed. The analytic study of the bifurcation point of the model then determines the volume V eb of the elongated bubble. Eventually, the magnitude of the pressure oscillations determines the fraction of gas that is not trapped in the elongated bubble.

A. Length of the riser

The riser is defined, in the model, as the part of the pipe downstream the virtual valve. The location of this valve is unknown a priori. In some cases, the geometry suggests a natural location, where an irregularity in the flow is likely to occur. It usually takes the form of an angle, when the well has a near-horizontal part followed by an inclined one. Once again, in the case of a low-point angle (when the near-horizontal part is declining), this assumption is consistent with the works of Storkaas [START_REF] Storkaas | Control solutions to avoid slug flow in pipeline-riser systems[END_REF]. When nothing in the geometry suggests the existence of an irregularity, other considerations, such as the location of the sensors can determine the most suited location for the virtual valve. The geometry of the real well considered in this article, depicted on figure 5, clearly indicates that the virtual valve must be placed at the end of the near-horizontal section. This location coincides with that of the downhole pressure sensor.

B. Straightforward parameters

Numerous parameters in the model can be directly inferred from the geometry of the system, or the nature of the liquidgas mixture. For example, the inclination parameter θ is the mean inclination of the pipe, and the cross-section area A is known. Similarly, the separator pressure and the temperature inside the pipe are measured and almost constant in practice. As a result, values for θ, A, g, ρ l , R, T , M and p s can be obtained from the geometry of the well under consideration and PVT data sheets for the fluid characteristics.

C. Inflow rates

Next, values must be given to the inflow rates of liquid (w l,in ) and gas(w g,in ). The difficulty is that they are scarcely measured. Besides, they are not constant: the assumption that they are is a simplification, which is formulated to ease the computations. We use in the model the values of the time-average outflow rates, commonly measured by Multiphase Flow Meters (MPM). At steady-state, the inflow and outflow rates must be constant and equal. Even though the average values of the outflow rates do not match their equilibrium values 3 , they lie in the same range. This reveals accurate enough for the model to reproduce the behavior of the original system. When no flow meter can be used, an alternative solution is to use a commercial multiphase flow simulator like OLGA TM to compute the steady-state values of the inflow rates from the system characteristics.

D. Choke constants

As previously mentioned, the model contains two valve equations ( 4)-( 5): one for the virtual valve, and one for the production valve. The choke constant for the production valve is usually provided by its manufacturer. Yet, equation [START_REF] Dukler | A model for gas-liquid slug flow in horizontal and near horizontal tubes[END_REF] assumes that the density of the mixture is constant. Therefore, the equilibrium value of the topside pressure cannot perfectly fit the system for all values of the choke opening. For our model, it is given by pr,t = p s + w 2 l,in ρ l C 2 c u 2 . In order to find a for C c , we choose to match it with the one of the real system ȳr,t (ū) for one value of the choke opening ū. To measure the equilibrium topside pressure y r,t (ū), ū has to lie in the stable region ū < u * . On the other hand, the chosen equilibrium of the model should be as close as possible to the real one in the unstable region, where the flow must be stabilized. As a result, ū is taken as close as possible to the bifurcation point, but strictly in the stable region. The value of C c must be chosen so that

C c = w l,in ū √ ȳr,t -p s
The second choke constant appearing in ( 4) is that of the virtual valve. Because it is difficult to mathematically study its impact on the dynamics, a rough estimate has to be picked manually at this stage of the tuning procedure. It may then be adjusted afterwards, once the parameters have been chosen, to slightly reshape the simulated oscillations. A typical value for this parameter is C g = 10 -4 kg.s.m -1 .

E. Still mass of liquid

As mentioned in Section II, the parameter m l,still represents the (minimum) mass of liquid present in the riser at all times. Its value is chosen using the equilibrium value of the bottom hole pressure for one operating point ū, given by pr,b (ū) = pr,t (ū) + ( ml,r (ū) + m l,still ) g sin θ A

= pr,t (ū) + pr,t (ū) pr,t (ū) + b

x (ρ l LAm l,still ) + m l,still g sin θ A

= pr,t (ū) + pr,t (ū)LAρ l + b xm l,still pr,t (ū) + b x

g sin θ A (7)
where x = w g,in w l,in is the equilibrium gas-liquid mass ratio, and b = ρ l RT M . All the parameters except m l,still in this equation are known. Therefore, knowing (or, once again, estimating it with OLGA TM ) the value of the steady-state bottom hole pressure ȳr,b (ū) of the system for one operating point ū, one should set the parameter m l,still to

m l,still = ȳr,b -ȳr,t g/A sin θ 1 + ȳr,t b x - ȳr,t b x V r ρ l
The value of m l,still highly depends on the system under consideration: it can be 0, e.g. for short flowline risers, and can reach 1 2 V r ρ l , which corresponds to at least half the riser being permanently filled with liquid.

F. Elongated bubble

A value for the volume of the elongated bubble V eb has yet to be determined. This parameter is used to set the position of the bifurcation point of the model to that of the real system, which we note u * . More specifically, the bifurcation point is characterized by the appearance of two purely imaginary eigenvalues. To find the best value for V eb , one could write Routh's criterion for the Jacobian matrix J(u * ) (expressions for the Jacobian matrix coefficients are reported in Appendix), and solve the obtained equation with respect to V eb . Unfortunately, the criterion does not take the form of an analytically tractable expression, and V eb must be determined with a numerical solver (e.g. Newton's method). One can achieve this by solving the following problem

{λ(u * , V eb )} = 0 (8)
with respect to V eb , where λ is one of the two complex conjugate eigenvalues of the system (the third eigenvalue being always real). Figure 2 represents the variation, for u = u * , of the real part of λ as V eb increases. The corresponding locus of the eigenvalues (parametrized by V eb ) are plotted in Figure 3.

λ(u * , V eb ) λ(u * , V eb ) V eb = 100 Optimal V eb value u * = 0.4 Veb = 1
Fig. 3. Locus of the eigenvalues λ(u * ) and λ(u * ) when V eb varies. For low values of V eb , the eigenvalues are in the LHP. As the parameter increases, the eigenvalues cross the imaginary axis, which corresponds to the criterion plotted in Figure 2 being 0.

G. Split inflow of gas

Eventually, one must choose a value for , which determines the fraction of the amount of gas that accumulates into the elongated bubble. directly impacts on the magnitude of the pressure oscillations. In particular, it can be mathematically proven 4 that when the virtual valve is closed, i.e. when p eb < p r,b , the mass of liquid m l,r converges asymptotically to m l,r = pr,t pr,t + b x (ρ l LAm l,still ), and the topside pressure converges to pr,t . Yet, these values are only approached because the virtual valve does not remain closed in the observed oscillations. Therefore, the pressure at the base of the riser p r,b converges to a "pseudo-equilibrium" p r,b = pr,t + (m l,r + m l,still ) g sin θ A When the virtual valve opens again, the mass of liquid soon drops. This point is illustrated in Figure 4 (y max r,bȳr,t ) A g sin θm l,still Otherwise, when no measure of the downhole pressure is available, has to be tuned manually, with a trial-and-error approach to match the magnitude of the topside pressure oscillations. As indicated by the expression of m l,r , the magnitude of the oscillations is a decreasing function of .

H. Summary

Table I gives a summary of the tuning procedure, along with the corresponding values for the case we study here. These values, reproduced courtesy of Statoil, correspond to an actual 7731 m-long well in the North Sea, schematically depicted on Figure 5. As mentioned in Section III-A, the location of the virtual valve coincides with that of the downhole pressure sensor, at the end of the near-horizontal section. This assumption is validated by the resultant value of volume V eb = 48 m 3 . Indeed, the total volume of the pipe is, in the model, LA+V eb = 140.8 m 3 , which is very close to the actual total volume of the well V = 140.4 m 3 . One should notice that nothing, in the computation of V eb , guarantees this equality a priori. In general, it would be possible to obtain a total volume LA + V eb greater or smaller than the total well volume. A greater volume would correspond to an elongated bubble originating upstream the well (in the case A smaller value implies that the gas pocket occupies only a part of the volume of the well-section upstream the riser. This is also a valid assumption. Too large or inconsistent discrepancies should lead to reconsider the candidate value for L at the beginning of the identification process.

Variable

IV. Resulting dynamics and observer design

Figure 6 shows a comparison between the model dynamics resulting from the aforedescribed tuning procedure and actual well measurements. Because the length L coincides with the depth of the sensor, the reconstructed downhole pressure and the measurements can be directly compared. Should this length differ from the sensor depth, some compensation would have been considered, to provide fair comparisons. The results show the relevance of the model as the magnitude and frequency of the simulated oscillations almost match that of the real system. One should notice the scarce sampling rate of the downhole pressure measurements. This is frequent with downhole sensors which are often unreliable, tend to deviate and sometimes break down due to the harshness of their environment. We now show how to use the model in an observer structure to synchronize it with the actual system, and reconstruct a robust estimate of the downhole pressure.

A. Observer design

In this part, we consider the problem of estimating the whole state of system (1)-( 2)-( 3) when only the topside pressure p r,t is measured. This is an important topic from an engineering point of view, because this measurement is often available. First, we re-write the model in the new set of coordinates x = (x 1 , x 2 , x 3 ) = (m g,eb , p r,t , m l,r )

ẋ1 = f 1 (x) = (1 -)w g,in -C g max(0, ax 1 -x 2 -c(x 3 + m 3,still )) ẋ2 = f 2 (x) = b m ∆ 3 -x 3 w g,in + C g max(0, ax 1 -x 2 -c(x 3 + m l,still )) +w l,in x 2 b - m ∆ 3 x 3 x 2 b uC c ρm (x 2 -p s )       ẋ3 = f 3 (x) = w l,in -uC c ρm (x 2 -p s )
where a = RT MV eb and c = g/A sin θ. The measured output is now p r,t = x 2 . The observer equations read

ẋ1 = f 1 ( x) ẋ2 = f 2 ( x) -k( x2 -x 2 ) ẋ3 = f 3 ( x)
where k is the observer gain. Even though this resembles a high-gain structure, we were not able to come up with a proof of the convergence yet. This point will be the subject of future investigation. We only present here an experimental validation of this approach.

B. Simulations and experimental validation a) Simulation results:

The observer was first tested in simulation, on the model itself. Figure 7 shows the convergence of the trajectory of the observer to that of the plant. The observer was turned on 12 000 s after the start of the simulation.

b) Real well data: The observer was then tested against real measurements coming from the above mentioned well. Figure 8 shows a comparison between the measured pressures and the estimates computed by the observer using only the topside pressure measurement. Not only this method allows to synchronize the model with the system, but it also provides a good estimation of the downhole pressure. As presented in [START_REF] Di Meglio | Model-based control of slugging flow: an experimental case study[END_REF], these estimates can then be used in a control algorithm to stabilize the flow and suppress the slugging phenomenon. In particular, the estimated downhole pressure can replace the measurements, compensating for their large sampling time.

V. Conclusion

A method to tune the parameters of the model presented in [START_REF] Di Meglio | A first principle model for multiphase slugging flow in vertical risers[END_REF] was proposed. The method is based on physical considerations and a mathematical analysis of the model. The efficiency of this technique is illustrated by comparison with real slugging well data. It highlights the relevance of the model itself and its capability to reproduce the behavior of slugging wells. An observer structure, allowing to reconstruct unavailable or unreliable pressure measurements has been presented, along with its experimental validation. As presented in [START_REF] Di Meglio | Model-based control of slugging flow: an experimental case study[END_REF], a controller can be derived from the model and allow suppression of the slugging phenomenon. The controller was tested on an experimental multiphase flow loop in Porsgrunn, Norway, with good success. In particular, it was able to stabilize the flow around relatively high values of the choke opening. Transposed on an actual offshore well, this result is expected to yield significant production increases. Fig. 8. Comparison between the observer estimates for k = 0.05 and the real measurements. After only one oscillation, the system is perfectly synchronized in terms of frequency and amplitude. The observer allows to reconstruct the downhole pressure from topside measurements.

3 VFig. 2 .

 32 Fig.2. Value of {λ(u * , V eb )} as V eb increases. When the eigenvalues cross the imaginary axis, the criterion is exactly 0. The corresponding value of V eb should be used in the model.

Fig. 4 .

 4 Fig. 4. Convergence to a pseudo-equilibrium depending on a good approximation of the maximum riser base pressure over the cycle. When a measure of the bottom hole pressure is available, can be chosen to match the peak of the oscillations of p r,b with that of the system ȳmax r,b by setting = ȳr,t b x ρ l V r -(y max r,bȳr,t ) A g sin θ

Fig. 5 .

 5 Fig. 5. Schematic view of the considered well

  Information used Value R -R = 8.314 J.K -1 mol -1 g g = 9.81 m.s -2 T Topside temperature sensor T = 363 K L Geometry of the well L = 5200 m A Geometry of the well A = 1.77 × 10 -2 m 2 θ Geometry of the well θ = 0.78 rad ρ l PVT data sheet ρ l = 900 kg.m -3 M PVT data sheet M = 22 × 10 -3 kg.mol -1 p s Separator pressure sensor p s = 6.6 bar w l,in Multiphase flow meter w l,in = 11.75 kg.s -1 w g,in Multiphase flow meter w g,in = 0.82 kg.s -1 C c Equilibrium topside pressure C c = 8.32 × 10 -2 m at ū = 35%: ȳr,t = 8.23 bar m l,still Equilibrium downhole pressure m l,still = 3.73 × 10 4 kg at ū = 35%: ȳr,b = 170 bar ≈ 1 3 ρ l V r V eb Bifurcation point u * ≈ 20% V eb = 48 m 3 Max value of the downhole = 0.78 pressure ȳmax r,b = 176 barTABLE I List of parameters used in the model

Fig. 6 .

 6 Fig. 6. Topside and downhole pressure oscillations. Comparison between the (open-loop) model and well dynamics at u = 35%.

Fig. 7 .

 7 Fig. 7. Convergence of the observer to the model trajectory. The observer only needs one oscillation to be synchronized with the plant, and a few more to converge completely.

  F. Di Meglio (corresponding author) is PhD candidate in Mathematics and Control at Centre Automatique et Systèmes, Unité Mathématiques et Syst'emes, MINES Paristech, 60, Bd St-Michel, 75272 Paris, Cedex 06, France G.O. Kaasa is Research Engineer at Statoil ASA, Research Center Porsgrunn, Heroya Forksningspark 3908 Porsgrunn, Norway N. Petit is Professor at Centre Automatique et Systèmes, Unité Mathématiques et Systèmes, MINES Paristech, 60, Bd St-Michel, 75272 Paris, Cedex 06, France V. Alstad is Research Engineer at Statoil ASA, Research Center Porsgrunn, Heroya Forksningspark 3908 Porsgrunn, Norway

  . Therefore, p

												r,b	is
	7.5 8 8.5 9 4 4.1 4.2 1.65 Riser top pressure Riser base pressure 1.7 1.75 Mass of liquid	1 x 10 5 1 x 10 4 1 7 x 10	1.1 1.1 1.1	1.2 1.2 1.2	1.3 1.3 1.3	1.4 1.4 1.4	1.5 1.5 1.5	1.6 1.6 1.6	1.7 1.7 1.7	1.8 1.8 1.8	1.9 1.9 1.9	2 4 x 10 2 x 10 4 2
												x 10 4
		1	1.1	1.2	1.3	1.4	1.5	1.6	1.7	1.8	1.9	2
												x 10

Similarly, Bøe[START_REF] Bøe | Severe slugging characteristics ; part i:flow regime for severe slugging ; part ii: Point model simulation study[END_REF] and Taitel[START_REF] Taitel | Stability of severe slugging[END_REF] derived, from force balance equations expressed on the body of one slug, conditions on the superficial velocities of gas and liquid to predict whether a given system will experience slugging or not. Much later, these criteria were proved to correspond to a Hopf bifurcation by Zakarian[START_REF] Zakarian | Analysis of two-phase flow instablilities in pipe-riser systems[END_REF], who analyzed a simpler DAE model and showed that, when varying certain parameters, the stability properties of the system switch at a certain critical (bifurcation) point.

Unfortunately, they fail to predict the instability of some systems (e.g. purely vertical wells), in the sense that, in certain cases, they may predict a steady flow regime for wells that are actually experiencing slugging.

thanks to an analytic integration of (1)-(2) and its asymptotics.

Appendix

The equations linearized around an operating point u read d dt (m g,eb , m g,r , m l,r ) T = J(u)(m g,eb , m g,r , m l,r ) T , with