Frédéric Grognard
email: frederic.grognard@inria.fr

Andrei R Akhmetzhanov

Olivier Bernard

Optimal strategies for biomass productivity maximization in a photobioreactor using natural light

We address the question of optimization of the microalgal biomass long term productivity in the framework of production in photobioreactors under the influence of day/night cycles. For that, we propose a simple bioreactor model accounting for light attenuation in the reactor due to biomass density and we obtain the control law that optimizes productivity over a single day through the application of Pontryagin's maximum principle. The dilution rate is the main control, the input concentration being only used as secondary control to maintain the substrate concentration high. An important constraint on the obtained solution is that the biomass in the reactor should be at the same level at the beginning and at the end of the day so that the same control can be applied everyday and optimizes some form of long term productivity. Several scenarios are possible depending on the microalgae's strain parameters and the maximal admissible value of the dilution rate: bang-bang or bang-singular-bang control or, if the growth rate of the algae is very strong in the presence of light, constant maximal dilution. A bifurcation diagram is presented to illustrate for which values of the parameters these different behaviors occur. Finally, a simple sub-optimal bang-bang strategy is proposed that numerically achieves productivity levels that almost match those of the optimal strategy.

Introduction

Microalgae have recently received more and more attention in the frameworks of CO 2 fixation and renewable energy [START_REF] Huntley | CO2 mitigation and renewable oil from photosynthetic microbes: A new appraisal, Mitigation and Adaptation Strategies for[END_REF][START_REF] Chisti | Biodisel from microalgae[END_REF]. Their high actual photosynthetic yield compared to terrestrial plants (whose growth is limited by 1 CO 2 availability) leads to large potential algal biomass productions of several tens of tons per hectare and per year [START_REF] Chisti | Biodisel from microalgae[END_REF]. Also, they have a potential for numerous high added value commercial applications [START_REF] Spolaore | Commercial applications of microalgae[END_REF] as well as wastewater treatment capabilities including nutrients and heavy metal removal [START_REF] Hoffmann | Wastewater treatment with suspended and non-suspended algae[END_REF][START_REF] Oswald | Micro-algal Biotechnology[END_REF]. Microalgae are cultivated at industrial scale in devices ranging from the most simple and cheap open raceways to high-tech photobioreactors. In this paper, we focus on the industrial production of microalgae in planar systems, which will be called photobioreactors by abuse of language.

The energy balance of the process is a crucial feature when dealing with bioenergy production. Thus, the only conceivable light source (i.e the primary energy source) is natural sunlight. Since it varies during the day, we have to consider this periodic feature in the design of the optimal control law. We are thus aiming at maximizing the yield of a photobioreactor operating in continuous mode with a periodic light source that has a light phase (day) and a dark phase (night). In addition to this time-variation, we will take the auto-shading in the photobioreactor into account: the pigment concentration (mainly chlorophyll) affects the light distribution and thus the biological activity within the reactor. As a consequence, for a too high biomass, light in the photobioreactor is strongly attenuated and per-capita growth is low; this also has a major impact on the productivity.

Optimal control of bioreactors has been studied for many years whether it was for metabolites production [START_REF] Tartakovsky | Optimal control of fed-batch fermentation with autoinduction of metabolite production[END_REF], ethanol fermentation [START_REF] Wang | Optimal feed policy for fed-batch fermentation of ethanol production by <small>zymomous mobilis</small&gt[END_REF], baker yeast production (Wu, Chen & Chiou 1985) or, more generally, optimal control of fed-batch processes taking kinetics uncertainties into account [START_REF] Smets | Optimal adaptive control of (bio)chemical reactors: past, present and future[END_REF]. The control of photobioreactors is however a lot more scarce in the literature, though the influence of self-shading on the optimal setpoint (Masci, Bernard & Grognard 2010) or on an MPC control algorithm [START_REF] Becerra-Celis | Nonlinear predictive control for continuous microalgae cultivation process in a photobioreactor, Control, Automation, Robotics and Vision[END_REF] for productivity optimization have already been studied. The light-variation was mostly absent [START_REF] Masci | Microalgal biomass productivity optimization based on a photobioreactor model[END_REF][START_REF] Becerra-Celis | Nonlinear predictive control for continuous microalgae cultivation process in a photobioreactor, Control, Automation, Robotics and Vision[END_REF] or considered to be an input that could be manipulated in order to impose the physiological state of the microalgae [START_REF] Marxen | A photobioreactor system for computer controlled cultivation of microalgae[END_REF] or maximize productivity as one of the parameters of bioreactor design [START_REF] Suh | Photobioreactor engineering: Design and performance[END_REF]. The optimization we consider can also be classified with the Maximal Sustainable Yield problem that is classical eg. in fisheries and has been tackled in constant and some form of specific periodic environments [START_REF] Clark | Mathematical Bioeconomics: The Optimal Management of Renewable Resources[END_REF][START_REF] Dong | Optimal harvesting policies for periodic gompertz systems[END_REF], and distincts itself from the classical periodic optimal control where periodicity was not inherent to the system [START_REF] Guardabassi | Status of periodic optimization of dynamical systems[END_REF]. The present result is an important generalization of the result of [START_REF] Grognard | Optimization of a photobioreactor biomass production using natural light[END_REF] which otherwise has not been tackled yet in the literature.

We developed a model that presents the main features of the process, with the substrate concentration in the input (marginally) and the dilution rate (mainly) as actuators for our control design. Since we want to develop a control strategy to be used on the long run, we could choose an infinite time-horizon measure of the yield. However, we rather impose the industrially relevant constraint that the control should be identical everyday and force the state of the system to be identical at the beginning of the day and 24 hours later, while optimizing a cost over one day.

The paper is structured as follows: in Section 2 we develop the photobioreactor model; in Section 3 we identify the optimal operating mode in constant light environ-ment; in Section 4 we develop our main result, that is the form of the optimal control law in a day & night setting; in Section 5 we concentrate on stability issues and possible implementation methods; we conclude by a simulation study and bifurcation analysis in Section 6.

A photobioreactor model with light attenuation

Micro-algae growth in a photobioreactor is often modeled through one of two models, the Monod model [START_REF] Monod | Recherches sur la Croissance des Cultures Bactériennes[END_REF] or the Droop Model [START_REF] Droop | Vitamin b12 and marine ecology. iv. the kinetics of uptake growth and inhibition in monochrysis lutheri[END_REF]). Here we focus on the Monod model which, though less accurate for microalgae, is more convenient for building control laws. It writes:

Ṡ = D(S in -S) -kµ(S, I 0 (t))X , Ẋ = µ(S, I 0 (t))X -DX -rX , (2.1)
where S and X are the substrate (e.g. nitrate concentration) and biomass (measured in carbon concentration, gC.L -1 ) in the medium, which is supposed to be homogeneous through constant stirring, while D is the dilution rate, S in is the substrate input concentration, µ(S, I 0 (t)) is the microalgae biomass growth rate with I 0 (t) the incoming time-varying light intensity, and k is the substrate/biomass yield coefficient. Compared with the classical model, the additional -rX term encompasses both biomass respiration, which has been known for a long time to impact productivity in bioreactors (see e.g. [START_REF] Grobbelaar | Respiration losses in planktonic green algae cultivated in raceway ponds[END_REF]) and mortality, which can be quite high in the high-density reactors that will be considered.

Modelling the growth rate

The function µ(.) is the average microalgae biomass growth rate. At any point of the photobioreactor, growth rate is related to local light intensity. For horizontal planar high density reactors, light is assumed to be exponentially decreasing with depth because of the auto-shading effect. We consider a constant horizontal section over the height of the reactor and vertical incoming light. Light intensity follows a Beer-Lambert law [START_REF] Lambert | Photometria sive de mensura et gradibus luminis, colorum et umbrae[END_REF] where the attenuation at some depth z comes from the total biomass Xz per surface unit contained in the layer of depth [0, z]:

I(Xz) = I 0 e -aXz , (2.2)
where I 0 is the incident light and a is a light attenuation coefficient. In microalgae, chlorophyll is mostly the cause of this shadow effect and, in model (2.1), it is best represented by a fixed portion of the biomass [START_REF] Bernard | A photobioreactor model in nitrogen limited conditions[END_REF], which yields the direct dependence in X in model (2.2). With such an hypothesis on the light intensity that reaches depth z, higher light causes higher growth in the upper part of the reactor than in the bottom part. Supposing that light attenuation directly affects growth [START_REF] Huisman | Principles of the light-limited chemostat: theory and ecological applications[END_REF], the growth rate for a given depth z, when the dependency in I are taken to be of Michaelis-Menten type, can then be written as

µ z (S, I(Xz)) = I(Xz) I(Xz) + K I µ S (S) .
with K I the associated half saturation constant and the substrate related function µ S (S) is such that µ S (0) = 0, is increasing and bounded (with lim S→∞ µ S (S) = μ). Generalization for the light-related part will also be proposed later on. This results in a spatially distributed growth-rate within the reactor. However, constant stirring keeps the concentrations of S and X homogeneous in the reactor. Then, we can compute the mean growth rate in the reactor:

µ(S, I 0 , X) = 1 L L 0 I 0 e -aXz I 0 e -aXz + K I dzµ S (S) = 1 aXL ln I 0 + K I I 0 e -aXL + K I µ S (S) .
where L is the depth of the reactor. This average growth rate will be used in the lumped model that we develop.

In order to determine more precisely the model, we have now to specify the varying light pattern I 0 (t). Classically, it is considered that daylight I 0 (t) varies, during the light period, as the square of a sinusoidal function. However, the introduction of such a varying light would render the computations analytically intractable. Therefore, we approximate the light source by a step function:

I 0 (t) = Ī0 , 0 ≤ t < T -light phase , 0, T ≤ t < T -dark phase .
with T the length of the day (which will be equal to 1 if the time-unit is the day). At the equinoxes, we have that T = T 2 , but this quantity obviously depends on the time of the year. We will now provide extensions to this model so that it fits better to the production problem in high-density photobioreactors under study.

Reduction and generalization of the model

The system for which we want to build an optimal controller is therefore

   Ṡ = D(S in -S) -k aXL ln I0(t)+K I I0(t)e -aXL +K I µ S (S)X Ẋ = 1 aXL ln I0(t)+K I I0(t)e -aXL +K I µ S (S)X -rX -DX (2.3)
However, in order to maximize the productivity it is clear that the larger S the better since it results in larger growth rates of the biomass X. Hence, the control of the inflow concentration of the substrate S in should always be kept very large so as to always keep the substrate in the region where µ S (S) ≈ μ. We can then concentrate on the reduced model

Ẋ = μ aL ln I 0 (t) + K I I 0 (t)e -aXL + K I -rX -DX ,
where the only remaining control is the dilution D and which then encompasses all the relevant dynamics for the control problem. As was seen in [START_REF] Masci | Microalgal biomass productivity optimization based on a photobioreactor model[END_REF], the relevant concentration in the photobioreactor, with Beer-Lambert light attenuation, is not the volumic density X but rather the surfacic density x = XL: the evolution of this quantity is indeed independent of the reactor's depth: whether we consider a deep and very diluted reactor or a shallow high concentration reactor with identical surfacic density

ẋ = d(XL) dt = μ a ln I 0 (t) + K I I 0 (t)e -ax + K I -rx -Dx , (2.4)
which is indeed independent of the depth L. Model (2.4) is quite simple except for the only nonlinear term which directly comes from the form of µ I and from the Beer-Lambert law (2.2). In order to generalize our approach to more light responses (e.g. for high density photobioreactor with possible high light inhibition, see [START_REF] Bernard | Hurdles and challenges for modelling and control of microalgae for CO2 mitigation and biofuel production[END_REF])) and to not restrict ourselves to the Beer-Lambert law, we propose a generalized model that retains the main features of (2.4) by taking advantage of the special form of the light source: during the day of length T < T , the system is written as

ẋ = f (x) -rx -ux , (2.5)
with f (x) the total growth and u the dilution rate; during the night of length T -T , the dynamics become ẋ = -rx -ux .

(2.6)

In order to couple both these systems, we define h(t) as a step function whose value is 1 for t < T and 0 for t ≥ T that will allow to synthetize (2.5)-(2.6) in the form that is

ẋ = f (x)h(t) -rx -ux .
(2.7)

Assumption 1 The growth function satisfies

(i) f : IR + → IR is C 2 , upper-bounded, satisfies f (0) = 0, f (0) > 0, f (x) < 0 for all x ≥ 0
(ii) f (0)T > rT with r > 0 and T > T > 0.

Such a choice of f (x) yields the following somewhat trivial property:

Property 1 If f (x) satisfies Assumption 1-(i), we have that, for any x > 0:

f (0) > f (x) x > f (x) and d dx f (x) x < 0 ,
which also shows that lim inf x→+∞ f (x) ≤ 0.

Proof: We have f (x) = f (0) + x 0 f (s)ds with f (0) = 0 and f (0) > f (s) > f (x) for all 0 < s < x due to the concavity, so that f (x) < x 0 f (0)ds = f (0)x and f (x) > x 0 f (x)ds = f (x)x. The decreasing property of f (x)
x comes from the explicit computation

d dx f (x) x = f (x)x-f (x)
x 2

which is negative thanks to the previous property.

Since the function f (x) is upper-bounded, the ratio f (x)/x tends to zero as x → +∞, which means that lim inf x→+∞ f (x) ≤ lim x→+∞ f (x)

x = 0. The four properties that we impose on f (x) should be expected from the net-growth in a photobioreactor: no growth should take place in the absence of biomass, it should be increasing because additional biomass should lead to more growth (at least at low densities); it should be bounded because, when x is very large the bottom of the reactor is in the dark, so that adding new biomass simply increases the dark zone without allowing additional growth; and the growth-rate f (x) x should be decreasing because additional biomass slightly degrades the environment for all because of the shadowing it forces.

Lastly, in practice, u cannot take any value: it should be positive, but also upperbounded since its value is determined by the maximum capacity of some pumps. We will then consider that 0 ≤ u(t) ≤ ū at all times (with ū > 0).

Biomass flow rate optimization in constant light environment

In a constant light environment, the dynamics are strictly imposed by (2.5) and we will be able to exploit the system at a steady state. In this section, the best constant operating conditions are identified: the constant control u and corresponding equilibrium value x(u), when it exists, are chosen in order to maximize the biomass flow rate at the output of the reactor. This is defined as max u ux(u)A (with A the horizontal section of the reactor).

Since ux = f (x) -rx at equilibrium, maximizing uAx amounts to maximizing f (x)-rx. f (x) being concave, this is equivalent to find the unique solution of f (x) = r when it exists. Note that if r = 0, it has no solution: the smaller the positive constant u, the larger the productivity that can approach sup x>0 f (x) as much as desired. This case hints at difficulties arising when r is very small. Since lim x→+∞ f (x) ≤ 0 (see Property 1) and f (.) is a decreasing function of x because of the concavity of f , f (x) = r has a positive solution if f (0) > r is satisfied (which is a consequence of ). This condition is crucial since, without it, the only equilibrium of the system is x = 0, independently of the choice of u. Under this assumption, the biomass density that corresponds to the maximization of the biomass flow rate in a constant light environment through constant control is

x σ = (f ) -1 (r) with u σ = f (x σ ) x σ -r > 0 . (3.8)
so that that f (x) -rx is increasing for x < x σ and decreasing for x > x σ , which we will extensively use.

Taking into account the actuator constraint, if u σ ≤ ū, u σ yields the optimal productivity. Otherwise, if u σ > ū, the actuator is under-dimensioned, and it can be computed that the best control is u = ū.

The equilibrium x(u) is a decreasing function, defined only for u < f (0) -r, because it is solution of f (x(u)) x(u) -r = u with f (x) x decreasing from f (0) at x = 0. Also, x(u) is globally asymptotically stable in IR + for 0 < u < f (0)-r; indeed, since f (x)
x is a decreasing function of x, ẋ = x( f (x)

x -r -u) < 0 (resp. > 0) for x > x(u) (resp. x < x(u)), which shows global asymptotic stability of x(u). For convenience, we will define two other equilibria beyond x σ : x0 = x(0) > x(ū) = xū , the non trivial equilibria of (2.5) with u = 0 and u = ū < f (0) -r as control, respectively.

Note that the study of the present section is in line with [START_REF] Masci | Microalgal biomass productivity optimization based on a photobioreactor model[END_REF], where productivity optimization in a constant light for a Droop model with light attenuation was considered. In that study, the analysis was much complicated by the link between shading and nitrogen content of the algae, so that both S in and D had to be taken into account; this inherent complexity prevents the direct extension of that result to the day-night setting.

Productivity optimization in day/night environment

In an environment with varying light, a non constant input must be considered. Here we consider the case where the photobioreactor is operated on the long term, with a daily biomass production from the reactor outlet. The problem that we thus consider is the maximization of the biomass production over a single day

max u(.) ∈ U T 0 (u(t)A)x(t)dt ,
where U is the set of measurable functions such that 0 ≤ u(t) ≤ ū ∀t ∈ [0, T ]. We are then looking for a periodic regime, where the photobioreactor is operated identically each day. This means that the initial condition at the beginning of the day should equal the final condition at the end of the day. This then requires that we add the constraint

x(T ) = x(0) .
In actual applications, the length of the bright phase will change slightly from one day to the next. This will probably impose a slight change of biomass at the beginning of the next day but, in this preliminary study, we will consider that such a phenomenon has little effect on the qualitative form of the solutions.

We therefore are faced with the following Optimal Control Problem:

max u(t) ∈ [0,ū] T 0 u(t)x(t)dt with ẋ = f (x)h(t) -rx -ux , x(T ) = x(0) .
(4.9)

.

Dealing with the T-periodicity

In order to solve this problem, it is convenient to observe that x(T ) = x(0) = x 0 cannot be achieved for all values of x 0 even without requiring optimality. For that, we first consider u = 0 for all t, which yields the largest value of x(T ) since no biomass is taken out of the system. We first see that, for x 0 ≥ x0 (the equilibrium of system (2.5) with u = 0), we have x(T ) < x 0 because x(t) is then decreasing all along the solution (or at least constant and then decreasing for x 0 = x0 ). The initial condition for the solution of (4.9) must then satisfy x 0 < x0 . This also imposes that, in the solution of (4.9), x(t) < x0 for all t. Indeed, in the bright phase, ẋ ≤ 0 at x = x0 (which at best allows for convergence of x(t) to x0 in infinite time), while ẋ < 0 in the dark phase.

In order to guarantee that the periodicity can be achieved for some x 0 > 0, we will then concentrate on what happens for x 0 in a small neighborhood of 0. Property 1 implies that the (2.5) dynamics with u = 0 are upper-bounded:

ẋ < (f (0) -r)x ,
for x > 0, so that x(T ) < x 0 e (f (0)-r)T and x(T ) = x(T )e -r(T -T ) < x 0 e (f (0)-r)T e -r(T -T ) . In order to have the upper-bound on x(T ) larger than x 0 (for all x 0 > 0), the exponential factor needs to be larger than 1, that is Assumption 1-(ii) needs to be satisfied. Considering now a small interval (0, ] for the initial condition, x(t) ≤ x(T ) < e (f (0)-r)T = M , for all t ∈ [0, T ]. We then have, from Property 1, that f (x) > f (x)x > f ( M )x for x ∈ (0, M ]. For all t ∈ [0, T ], we then have ẋ > (f ( M )r)x and x(T ) > x 0 e (f ( M )-r)T , so that x(T ) > x 0 e (f ( M )T -rT ) . In order to have x(T ) > x 0 , it then suffices to have f ( M )T -rT > 0 for some > 0. Provided Assumption 1-(ii) is satisfied, it suffices to take small enough, and x(T ) > x 0 for all x 0 ∈ (0, ]. Assumption ( 1)-(ii) is therefore both necessary and sufficient to ensure that there exists some x 0 that has x(T ) > x 0 with u(t) = 0 for all t.

We see there the importance of this assumption, which is quite natural since it imposes that, when the population is small, that is when the growth rate is the largest, growth during the daylight period exceeds the net effect of respiration that takes place during the whole day.

Maximum principle

In this section, we will show that the solution of (4.9) can have one of three patterns. In order to solve problem (4.9), one option would be to use the maximum principle for hybrid systems [START_REF] Sussmann | A maximum principle for hybrid optimal control problems, Decision and Control[END_REF]) but, given the simplicity of our setting, we will simply use the original Pontryagin's Maximum Principle (PMP, Pontryagin, Boltyansky, [START_REF] Pontryagin | Mathematical Theory of Optimal Processes[END_REF]) in looking for a control law maximizing the Hamiltonian

H(x, u, λ, t) λ (f (x)h(t) -rx -ux) + ux , with the constraint ẋ = f (x)h(t) -rx -ux , λ = λ (-f (x)h(t) + r + u) -u .
(4.10)

In addition, we should add the constraint

λ(T ) = λ(0) .
as shown in [START_REF] Gilbert | Optimal periodic control: A general theory of necessary conditions[END_REF]. We see from the form of the Hamiltonian that

∂H ∂u = (1 -λ)x ,
so that, when λ > 1, we have u = 0, when λ < 1, we have u = ū, and when λ = 1 over some time interval, intermediate singular control might be applied. Note that, despite the discontinuity in time, no jump in the adjoint variable occurs; such a jump would only occur if using the Hamiltonian for the system extended with time as a second state, and would only concern the corresponding second adjoint state. Due to the simplicity of the time-dependence, the use of the extended system was not needed.

Of paramount importance in the proofs will be the constancy of the Hamiltonian. Indeed, it is known that the Hamiltonian is constant along optimal solutions as long as the time does not intervene into the dynamics or the payoff. Since this constancy is rooted in infinitesimal computations (in regions where u(t) is constant, it suffices to check Ḣ), this constancy is preserved in the intervals [0, T ) and (T , T ]. The Hamilto-nian presents a discontinuity at T , which would not have occurred if using the Hamiltonian for the system with time as a second state. The adjoint variables, however, remain continuous in time at that point through the maximum principle, though their derivatives can present a discontinuity.

We can first give some general statements about where and when switches can occur Proposition 1 If Assumption 1 is satisfied then, a solution of problem (4.9) satisfies the following properties (i) No switch from u = 0 to u = ū can take place in the dark phase;

(ii) Switches from u = 0 to u = ū in the bright phase take place with x ≤ x σ ;

(iii) Switches from u = ū to u = 0 in the bright phase take place with x ≥ x σ .

Proof: All these results come from the analysis of λ at the switching instant. Indeed, at a switching instant, we have λ = 1 so that:

λ = λ (-f (x)h(t) + r) .
(4.11)

The form of ∂H ∂u indicates that a switch from u = 0 to u = ū (resp. from u = ū to u = 0) can only occur if λ ≤ 0 (resp. λ ≥ 0).

In the dark phase, (4.11) becomes λ = rλ > 0; no switch from u = 0 to u = ū can therefore take place (i).

In the bright phase, λ ≤ 0 at a switching instant if f (x) ≥ r, which is equivalent to having x ≤ x σ ; this is therefore a condition for a switch from u = 0 to u = ū (this shows (ii)). Conversely for (iii).

Point (i) is in fact quite natural: if a bioreactor was closed at some moment of the night, opening it later during the night would mean that some biomass would have been uselessly respirated before that. The following theorem gives all possible solutions of the PMP Theorem 1 If Assumption 1 is satisfied then three forms of control are possible in solutions of problem (4.9):

(a) Bang-bang with u(0) = 0, a single switch to u(t) = ū taking place strictly before T and a single switch back to u(t) = 0 taking place strictly after T ;

(b) Bang-singular-bang with u(0) = 0, a switch to u(t) = u σ taking place first, a single switch to u(t) = ū taking place strictly before T and a single switch back to u(t) = 0 taking place strictly after T .

(c) Constant control at u(t) = ū

In the first two cases, the switch back to u(t) = 0 takes place with x strictly smaller than it was at the moment of switch to u(t) = ū.

Without needing an explicit form of f (x), we have been able to obtain the qualitative form of the optimal solution analytically. The proof of this result is detailed in appendix and is obtained by considering all possible situations that are in concordance with Proposition 1 and eventually eliminating all possibilities but the three cases detailed in Theorem 1. The link with the optimal solution for the day-night version of model (2.3) is then completed by imposing a very high S in during the open-reactor phases so that the substrate density is kept high and µ S (S) ≈ μ at all times.

In the bang-bang case, the optimal solution is made of four phases: (i) growth with a closed photobioreactor until a sufficient biomass level is reached -(ii) maximal harvesting with simultaneous growth -(iii) maximal harvesting, without growth, of the biomass produced during the light phase while not going below the level where the residual biomass left is sufficient to efficiently start again the next day -(iv) passive bioreactor (no harvesting, no growth, only respiration).

In the bang-singular-bang case, a (i-bis) singular phase is added between (i) and (ii) compared with the bang-bang case: during this whole phase, maximal instantaneous productivity is achieved. This solution with a singular form then seems to be naturally the most efficient one. It can however not always be achieved for two reasons:

• if u σ > ū, u σ is not an admissible control. The phase where u = ū in the bangbang case then parallels the singular phase since, in both cases, they correspond to instantaneous productivity optimization.

• if growth is not sufficiently stronger than respiration (with Assumption 1-(ii) satisfied however), there might not be enough time for x to reach x σ during the light phase.

The only remaining solution is the one with u constant at ū and is characterized by constant maximal harvesting when the actuator has been under-dimensioned and when the dark phase does not last too long, otherwise the respiration and/or the maximal dilution would wash-out the reactor. Finally, note that the phases where the reactor is closed at the beginning of the day is there to allow for the biomass to go back up to the level where productivity is high, so that, when the night duration goes to zero, the phases where u = 0 disappear. A similar reasoning can be held for the u = ū phase in case (b), which harvests the biomass at maximal rate to limit respiration at night. There is thus clear continuity between these day-night solutions as the night length goes to zero and the constant light optimal ones.

Discussion

Existence and stability of the optimal solution

Existence of an optimal solution is obvious since the achieved yield for solutions satisfying x(T ) = x(0) is bounded between 0 (obtained for u = 0) and x0 ūAT (not achievable but an upper-bound nonetheless because x(t) ≤ x0 and u(t) ≤ ū at all times). The productivity level therefore has a finite supremum, which translates into a maximum since the yield is continuous in u(t) and the set of definition of the u(t) control laws bounded by [0, ū] is compact. The optimal control problem therefore has a solution. Note that if f (x) was linear, no solution of the optimal control problem would exist because, for a given u(t) that ensures x(T ) = x(0), linearity yields that this is achieved for any x(0) with a payoff linear in x(0). Optimality then requires x(0) → +∞.

Set in a periodic framework, where it would repeat itself over time, we will show that the optimal solution is also asymptotically stable. Indeed, for given periodic u * (t) and x * (t) (defined for t ∈ IR + ) built on an optimal control law of Theorem 1 and the corresponding optimal solution, the solution of the periodic use of equation (2.7), with u = u * (t) and some x(0) > 0, generates a solution x(t) (defined for t ∈ IR + ). Convergence of x(t) to x * (t) is then shown by proving that the (x(nT )) n ∈I N sequence is converging to x * (0) > 0. For that, we consider x(0) = x * (0) (which implies x(t) = x * (t) at all times through uniqueness of solutions considerations). Building now the function

W (x(t), x * (t)) = 1 2 (ln(x(t)) -ln(x * (t))) 2 ≥ 0, its time derivative along (2.7) is Ẇ = h(t) f (x(t)) x(t) - f (x * (t)) x * (t) (ln(x(t)) -ln(x * (t)))
which is zero in the dark phase and negative in the light phase because f (x)

x is a decreasing function while ln(x) is an increasing function of x. We conclude that W (x((n + 1)T ), x * ((n + 1)T )) < W (x(nT ), x * (nT )) which, since x * is T -periodic. is equivalent to W (x((n + 1)T ), x * (0)) < W (x(nT ), x * (0)). We then have that (W (x(nT ), x * (0))) n ∈I N ≥ 0 is a decreasing, hence converging, sequence. The asymptotic behavior of (x(nT )) n ∈I N is then characterized by having (W (x(nT ), x * (0))) n ∈I N constant. Since constancy requires to have x(t) = x * (t) in the whole light phase to keep Ẇ = 0, (x(nT )) n ∈I N necessarily converges to x * (0), so that the periodic x * (t) is also globally asymptotically stable. This property is quite useful since it guarantees that, when applying the periodic optimal control over several days, the biomass level will automatically be steered to the level that ensures optimal yield.

Practical implementation

The computation of the solution of the optimal control problem (4.9) gives us two things: u * (t) (t ∈ [0, T ]), the optimal control, and x * (t) (t ∈ [0, T ]) the biomass level evolution. The latter fixes the initial and final condition x * (0) = x * (T ). Once this is set, two points are important during the implementation on a real plant: the final condition should be the desired x * (T ) so that production can start again the next day and the biomass productivity should be maximized. In the absence of perturbation and with the correct x(0) = x * (0), this is easily achieved through u * (t). If the system suffers from disturbances, inaccurate initial conditions or presents some unmodelled dynamics, neither will be achieved. With the hypothesis that the previously computed x * (T ) is desirable despite these perturbations, a simple solution to limit their consequences would be to set our problem in a setting close to a Model Predictive Control framework [START_REF] Rawlings | Tutorial overview of model predictive control[END_REF], but without receding horizon: the time interval [0, T ] is divided into m (for some fixed m ∈ IN ) smaller intervals of length T m at the beginning of which the solution of the problem

max u(t) ∈ [0,ū] T jT m (u(t)A)x(t)dt ,
with x( jT m ) given and x(T ) = x * (T )

(5.12) is computed, for j ∈ {0, • • • , m -1} and with the x dynamics given in (2.7). In detail, once a control (computed by the present MPC-like scheme) has been applied up to time jT m on the real plant, the biomass level x( jT m ) is measured and problem (5.12) is either solved online through an efficient algorithm (e.g. shooting method [START_REF] Betts | Practical Methods for Optimal Control and Estimation Using Nonlinear Programming, Advances in design and control[END_REF])) or uses the pre-computed solution of (5.12) that would have been achieved offline for all initial conditions. The obtained control u(t) is then applied between times x( jT m ) and x( (j+1)T m ), where the procedure is started again. Note that the solution of (5.12) is linked to the value function of the original problem since it gives the cost-to-go at time jT m ; this is a link between this scheme and a construction of an explicit feedback built on the Hamilton-Jacobi-Bellman equation, which hinges on the explicit construction of the value function.

Note that, due to perturbations having occured or unmodelled dynamics, there might be some j for which problem (5.12) has no solution because the terminal constraint x(T ) = x * (T ) cannot be satisfied for any control. In that case, the productivity objective is dropped and problem (5.12) is replaced by min u(t) ∈ [0,ū] |x(T ) -x * (T )| for the given x( jT m ). Monotonicity of x(T ) with respect to u(t) implies that the solution of this problem is trivially u(t) = ū (resp. u = 0) if the constraint x(T ) = x * (T ) in (5.12) cannot be satisfied because x(T ) > x * (T ) (resp. x(T ) < x * (T )) independently of u(t). The procedure is then started again at time (j+1)T m by first evaluating if (5.12) has a solution. This strategy is chosen because the terminal constraint is deemed more important than the yield for the perenity of the process.

Simulations

In this section, we concentrate on model (2.4), which means that the high substrate level hypothesis is supposed to be valid through a high value of S in in model (2.3). Simulations of the solution of (4.9) are then performed with f (x) = μ a ln Ī0 + K I Ī0 e -ax + K I .

Results for the microalgae Isochrysis galbana

We will now analyze the forms of solutions (a), (b), and (c) of Theorem 1. For that, we start with a dynamical model (2.4) for the growth of Isochrysis galbana with the parameters taken as in [START_REF] Masci | Microalgal biomass productivity optimization based on a photobioreactor model[END_REF]).

With such parameters, the optimal solution, represented by the blue curves in Figure 1 is of the bang-singular-bang type and achieves a daily surfacic productivity of 6.33g[C]/m 2 for a total cumulated flow 1 0 u(τ ) dτ equal to 0.453, that is 45% of the medium has been renewed during the 24 hours. We then considered the application of a constant control during the 24 hours and optimized the level of this control numerically. The optimum was achieved for û = 0.461, which yields a daily productivity

Parameter Value Units μ 1.7 day -1 a 0.5 m 2 /g[C] Ī0 1500 µmol quanta m -2 s -1 K I 20 µmol quanta m -2 s -1 r 0.07 day -1 T 0.5 day ū 2 day -1
Table 1: Growth and bioreactor parameters for Isochrysis galbana equal to 6.26g[C]/m 2 and a cumulated flow equal to 0, 461 also. Two points need to be noticed from this comparison: (i) the productivity improvement with the optimal solution, compared to the best constant control, is very weak (1.11%); (ii) the total flows required to attain both optimal solutions are very similar. The fact that the improvement of the productivity is small is not surprising: the necessity of shutting down the chemostat at night is linked to the respiration that would consume the biomass. In the present case, the respiration is weak so that, during one night, only 3.44% of the biomass is consumed. This phenomenon is here marginal so that the optimal control approach developed to limit it provides little gain and what really matters is the total flow that goes through the photobioreactor. Almost any control that achieves large values of the biomass will yield near-optimal productivities because this large biomass is not penalized by respiration.

In order to cover all cases of Theorem 1, we considered two other situations: ū = 0.8 and ū = 0.1, which prevent the existence of a singular phase because ū < u σ . For ū = 0.8, the optimal solution is bang-bang with a productivity level very close to the previous case (6.30g[C]/m 2 ) and the total flow also varies little (0.457). With ū = 0.1, a total flow close to 0.46 cannot be achieved and a significant performance degradation occurs: u(t) is constant at u = 0.1 and the productivity drops to 4.35g[C]/m 2 .

Importance of the respiration

In this section we explore the impact of a larger value for parameter r, which can be due to increased respiration, or to a high mortality as often is the case in high density cultures (however, we will stick to the respiration terminology). If we now consider Isochrysis galbana in an environment where mortality brings r to 0.7day -1 , we expect the optimal strategy to have a much bigger impact on the outcome. Indeed, such a respiration level has much more importance at night than in the previous case since it consumes 29.5% of the biomass, hence the importance of limiting the biomass level at night. We see in Figure 2 that the optimal solution is here bang-bang with a short opening window at the end of the day and at the beginning of the night to harvest the produced biomass (u = ū for t ∈ [0.479, 0.527]). The optimal production is here 0.607g[C]/m 2 for a total flow of 0.096, that is a very little daily medium renewal while the best constant control û = 0.095 yields 0.519g[C]/m 2 . We see here that the daily total flows are again almost equivalent but that, because respiration is limited by the early harvest, the productivity is here improved by 17% through the bang-bang approach, which is far from being negligible, especially since it is made with almost x Figure 1: Bang-singular-bang optimal solution (in black) for model (2.4) confronted with the most productive constant dilution rate scenario (in grey) for the microalgae Isochrysis galbana and the parameters of Table 1. At the top is the evolution of the biomass and at the bottom that of the control. The black dash-dotted lines represent the values of x σ and u σ respectively the same cumulated through-flow, hence consumed energy, as the best constant control case.

Near-optimal strategies

We have seen that the daily flow that goes through the chemostat has great importance for the productivity level of a solution: in the small respiration case, the bang-bang controls with ū = 2 or ū = 0.8 and the constant control with û = 0.461 have almost the same productivity, while having very similar total flow going through the reactor. In order to exploit that property that, we propose a strategy that is of the bang-bang type and consists in having the dilution equal to ū in the interval

[T - ũ 2ū , T + ũ 2ū
] and 0 outside of this interval. That way, the total flow that goes through the reactor is equal to ũ, so that we will be able to compare the obtained productivity level between that strategy and constant control strategies that have the same daily total flow. Being bang-bang, this strategy has the form of an optimal solution defined in Theorem 1; however, the switching times are most probably not the ones that the PMP would prescribe, so that only near-optimality is achieved. We did the computations for ū = 2 and values of ũ that did not lead to the wash-out of the reactor for species of subsection 6.2 (Figure 3). We see that, the constant control at value ũ (in red) yields a less productive process than the corresponding suboptimal bang-bang strategy. Therefore, when the exact optimal control law is not computed, it is advantageous to choose a bang-bang one rather than a constant control. This strategy is strongly advisable since the optimal productivity level is represented through a dotted level in the figure and the near-optimal strategy achieves it almost for an appropriate daily flow. This was expected in the high respiration case where the optimal solution is bang-bang, as would also be the case in a low-respiration case where the optimal solution is bang-singular-bang, where respiration is not too detrimental if x(t) goes above x σ . Such would probably not be the case when confronting the near-optimal strategy to a bang-singular-bang optimal pattern that would occur in a high respiration case.

Conclusions

In this paper, we have shown how the day-night constraint influences the optimal control strategy that achieves maximal biomass daily productivity in a photobioreactor and reduces the optimal productivity level. We have identified three families of strategies that can achieve optimality: bang-singular-bang, bang-bang, and constant maximal control; the first two are characterized by an harvesting of the biomass at the end of the light phase and beginning of the night to limit the negative effects of the respiration, while the latter leads to permanent harvesting because the maximal dilution rate is too weak compared with the growth rate of the biomass. These families of control : Productivity levels of species of Section 6.2: in grey, the productivity attained with the constant control u = ũ and in black the one with the near-optimal strategy of Section 6.3. The dashed-line represent the optimal productivity levels with ū = 2 strategies have been built for a large set of nonlinear one-dimensional models with light-dark phases so that they can be applied beyond their motivation in this paper: that is a photobioreactor with Monod-like growth and Beer-Lambert light attenuation. Such a formalism even allows to go beyond chemostat applications and consider populations that grow in more general periodic environments (see [START_REF] Grognard | Periodic optimal control for biomass productivity maximization in a photobioreactor using natural light[END_REF] for an example).

Through simulations, we have shown that the necessity of applying an optimal control strategy strongly depends on the respiration rate: if the latter is weak, constant control can achieve almost the same performance since the night phase consumes little biomass. However, we have also shown that a better choice is probably to apply a bangbang control law which, in the presented simulations, always yields better productivity than comparable strategies with constant dilution. This is particularly supported by the fact this better productivity is achieved with a very similar cumulated effort as the one necessary for a constant dilution rate. Wu, W. T., Chen, K. C. & Chiou, H. W. (1985). On-line optimal control for fedbatch culture of bakers yeast production., Biotechnology and Bioengineering 27(5): 756-760.

A Proof of Theorem 1

We will detail the different forms of solutions and show that no other solution than the one exposed in Theorem 1 can occur, based on a discussion on λ 0 λ(0). Bang-bang with λ 0 > 1: With λ 0 > 1, we have u = 0 at times 0 and T . At least, a switch of u from 0 to ū and back to 0 then needs to occur between time 0 and T ; otherwise, u(t) = 0 at all times and the payoff T 0 (u(t)A)x(t)dt would be 0, since the argument of the integral that defines the productivity would always be 0. Since the productivity is non-negative by constructive, bang-bang solutions that switch to ū can do no worse than that and are worth examining. As have seen in , the switch from 0 to ū cannot take place in the dark phase. For the solution that we study, a switch then needs to take place at time t 0ū in the (0, T ) interval and for x(t 0ū ) = x 0ū < x σ and λ(t 0ū ) = 1 (Proposition 1-(ii)).

The switch from u = ū to u = 0 then needs to occur at some time t ū0 > t 0ū (and some value of x(t ū0 ) = x ū0 ). Two possibilities then need to be considered:

• t ū0 ≤ T : from Proposition 1, we can deduce that x ū0 ≥ x σ . Applying u = 0 then forces convergence toward x0 for t ū0 ≤ t ≤ T , so that x(t) increases, and so stays above x σ (so that x(T ) x T ≥ x σ > x 0 ). No switch back to u = ū can then take place neither before T (because of (ii) of Proposition 1) nor after T (because of (i) of Proposition 1).

• t ū0 > T : from Proposition 1-(i), we can deduce that no other switch takes place afterward. Two cases will then be considered: either x 0ū ≥ xū so that x(t) decreases and x 0ū ≥ x T > x ū0 or x 0ū < xū so that x 0ū < x T , and we do not necessarily know if x 0ū > x ū0

We have then shown that any bang-bang solution with u(0) = 0 presents a single switch to u = ū (before T ) and a single switch to u = 0 (before or after T ). We are then left with two things to show: first that t ū0 ≤ T leads to a contradiction and then that, when x 0ū < xū , we indeed have x 0ū > x ū0 .

We will now show that the switch from ū to 0 needs to occur strictly after T . For that, we will consider that t ū0 ≤ T and show that this leads to a contradiction. Using the constancy of the Hamiltonian in [0, T ), we have, by continuity of the variables at time T , that

λ 0 (f (x 0 ) -rx 0 ) = λ T (f (x T ) -rx T ) ,
and, in (T , T ], with x(T ) = x 0 and λ(T ) = λ 0 :

λ 0 (-rx 0 ) = λ T (-rx T ) .
Taking the differences between the equalities, we obtain

λ 0 f (x 0 ) = λ T f (x T ) .
Taking the quotient of the last two equalities then yields

f (x 0 ) x 0 = f (x T ) x T ,
Remembering from Property 1 that f (x)

x is a strictly decreasing function, this leads to a contradiction because x T > x 0 since the biomass only decreases at night. The switch from u = ū to u = 0 then takes place with t ū0 > T .

We will now show that, even in the case where x 0ū < xū , the second switch takes place with x ū0 < x 0ū . For that, we first note that, since no from u = ū to u = 0 takes place in the [t 0ū , T ] interval, f (x T ) -rx T > f (x 0ū ) -rx 0ū . Indeed, since x(t) is increasing along the considered solutions from x 0ū < x σ , f (x(t)) -rx(t) first increases from f (x 0ū ) -rx 0ū and then decreases (once x(t) gets above x σ ); if we had f (x T ) -rx T ≤ f (x 0ū ) -rx 0ū , there would be a time t belonging to the interval (t 0ū , T ] where f (x( t)) -rx( t) = f (x 0ū ) -rx 0ū . A switch then necessarily would have taken place at that instant because of constancy of the Hamiltonian

(f (x 0ū ) -rx 0ū ) = λ( t) f (x( t)) -(r + ū)x( t) + ūx( t) results, using f (x( t)) -rx( t) = f (x 0ū ) -rx 0ū , in (1 -λ( t)) f (x( t)) -rx( t) = (1 -λ( t))ūx( t) .
This can only be achieved if λ( t) = 1, which corresponds to a switching instant (and a contradiction) or f (x( t)) -rx( t) = x( t), which would mean that x( t) = xū , which also is impossible since, in that phase, x(t) was converging in infinite time toward xū (and x 0ū < xū ). We then conclude that, as announced, f (x T )-rx T > f (x 0ū )-rx 0ū .

We then use the constancy of the Hamiltonian in both to show that x ū0 < x 0ū . In the bright phase, we have

λ 0 (f (x 0 ) -rx 0 ) = (f (x 0ū ) -rx 0ū ) ,
and, in the dark phase λ 0 (-rx 0 ) = (-rx ū0 ) .

Taking the difference between these equalities yields λ 0 f (x 0 ) = (f (x 0ū ) -rx 0ū ) + rx ū0 . Finally we note that, after t ū0 , the dynamics become λ = rλ , so that, using λ(t ū0 ) = 1 and λ(T ) = λ 0 , we have λ 0 = e r(T -tū0) , which yields f (x 0 )e r(T -tū0)rx ū0 = f (x 0ū ) -rx 0ū . Concavity implies that f (x 0 )e r(T -tū0) > f (x 0 e r(T -tū0) ) = f (x ū0 ), so that the constancy of the Hamiltonian conditions can only be satisfied if

f (x ū0 ) -rx ū0 < f (x 0ū ) -rx 0ū .
This can only be achieved for x ū0 < x 0ū , since x ū0 < x T and f (

x T ) -rx T > f (x 0ū ) -rx 0ū > f (x ū0 ) -rx ū0 .
Bang-singular-bang with λ 0 > 1:

We will first look at what a singular arc could be. For that, we see that ∂H ∂u = (1 -λ)x should be 0 over a time interval which, since a biomass level of x = 0 does not make sense when optimizing the productivity, amounts to imposing λ = 1. We then compute its time derivatives. d dt (λ -1) |λ=1 = -f (x)h(t) + r .

When h(t) = 0, that is in the dark phase, no singular arc is thus possible. When h(t) = 1, this derivative is equal to zero when x = x σ defined in (3.8). The singular control is then the control that maintains this equilibrium, that is u σ defined in (3.8). This control is positive thanks to Assumption 1-(ii) but is only admissible if u σ < ū.

No singular control can exist otherwise. In fact, this implies that xū < x σ < x0 because x(u) is decreasing. The case where u σ = ū has been handled through the bang-bang case.

When a singular branch appears in the optimal solution, it has to be kept among local optimality candidates because the second order Legendre-Clebsch condition [START_REF] Kelley | Singular extremals[END_REF] ∂ ∂u

d 2 dτ 2 ∂H ∂u = -f (x σ ) ≥ 0
is satisfied on the singular arc because of the concavity.

The construction of the solution is very similar to that in the purely bang-bang case. Similarly, a first switch needs to occur in the interval (0, T ). This switch can be from u = 0 to u = ū or from u = 0 to u = u σ and should occur with x ≤ x σ (Proposition 1-(ii)). In fact, if a switch first occurs to u = ū, x then converges toward xū , so that it does not go on or above x σ , which prevents any other switch before T (Proposition 1-(iii)). No singular arc can then appear. These solutions have been handled earlier.

In bang-singular-bang cases, a switch from 0 to u σ then directly takes place once λ = 1 at (t 0σ , x σ ).

From there, λ(t) = 1 and x(t) = x σ for some time. This could be until t σ0 ≤ T , followed directly by u = 0 or the singular arc could end at time t σ ū < T , where a switch occurs toward u = ū; note that the strict inequality is due to the fact that, if we had t σ ū = T , we would then have λ(T ) = 1 and λ > 0 for t ≥ T , so that u directly goes to 0: no actual switch to u = ū has taken place. Using the constancy of the Hamiltonian, we can conclude that a direct switch from u σ to u = 0 is not possible at t σ0 ≤ T . Indeed, if such an early switch occurred, we would have λ 0 (f (x 0 ) -rx 0 ) = λ T (f (x T ) -rx T ) , λ 0 (-rx 0 ) = λ T (-rx T ) .

Taking the differences between the equalities, we obtain

λ 0 f (x 0 ) = λ T f (x T ) .
Taking the quotient of the last two equalities then yields f (x 0 )

x 0 = f (x T )

x T , which we have shown earlier be impossible. The optimal solution then leaves the singular arc with u = ū strictly before T , and switches to 0 strictly after T .

From then on, things are unchanged with respect to the bang-bang case.

Solution with λ 0 ≤ 1:

Since λ > 0 in λ = 1 in the dark phase, such a solution would mean that harvesting takes place during the whole dark phase because no transition from u = 0 to u = ū can take place in this phase (Proposition 1-(i)). Two possibilities then occur: either u = ū all the time or switches from u = ū to u = 0 or u σ and then back to u = ū take place in the interval (0, T ).

In the latter case, the first switch from u = ū to u = 0 can only take place with x > x σ (Proposition 1-(iii)). Then, when the control u = 0 is applied for some time, the solution x(t) is increasing. No switch back to u = ū can then take place before T because such a switch would require x(t) < x σ (Proposition 1-(ii)), which cannot occur. An optimal solution of this form cannot exist.

We can also show that no strategy in the (0, T ) interval can have the form u = ū → u σ → u = 0 or ū. Indeed, since 0 < u σ < ū, x0 > x σ > xū , the singular arc can only be reached with u = ū when coming from above it. If the switch that takes place at the end of the singular phase is from u σ to 0, x(t) will increase and there should be a subsequent switch from 0 to ū which is impossible with x(t) > x σ . If the switch that takes place at the end of the singular phase is from u σ to ū, x(t) will decrease all the time between t σ ū and T , which is in contradiction with the fact that we had x(0) > x σ . The only potential optimal control in that family is therefore u(t) = ū for all times.

Figure 2 :

 2 Figure2: Bang-bang optimal solution (in black) for model (2.4) confronted with the most productive constant dilution rate scenario (in grey) for a high-respiration species (r = 0.7) and ū = 2

  Figure3: Productivity levels of species of Section 6.2: in grey, the productivity attained with the constant control u = ũ and in black the one with the near-optimal strategy of Section 6.3. The dashed-line represent the optimal productivity levels with ū = 2

Acknowledgments

We are very grateful to the three reviewers who carefully read the manuscript and strongly contributed to its improvement.


Michaela Blessing
email: m.blessing@brgm.fr

Alain Saada

Stéfan Colombano

Jennifer Hellal

Marc Crampon

MONITORED NATURAL ATTENUATION OF A COMPLEX MIXTURE OF ORGANOCHLORINES USING COMPOUND-SPECIFIC ISOTOPE ANALYSIS

Monitored natural attenuation of a complex mixture of organochlorines using compound-specific isotope analysis.

Compound-specific isotope analysis (CSIA) by online-coupling of gas chromatography and isotope ratio mass spectrometry (GC-IRMS) offers a versatile tool to study the origin of a pollution and/or to decipher multiple contamination sources, to verify if biodegradation of organic pollutants is occurring, to identify degradation mechanisms, and to determine the rate and extent of in-situ biodegradation. In the last decade, compound-specific isotope analysis (CSIA) has evolved as a diagnostic tool for contaminated site investigation and for monitoring the efficiency of intrinsic natural or enhanced remediation processes, and is now available for studying an increasing number of environmental pollutants [START_REF] Elsner | Current challenges in compound-specific stable isotope analysis of environmental organic contaminants[END_REF]. Chemical compounds contain heavy and light stable isotopes in a certain ratio, i.e. 13 C/ 12 C, D/H, and 15 N/ 14 N. In general, molecules containing the lighter isotope react faster than those containing the heavier isotope. Bond-breaking processes induced by microbial metabolism, can cause a large isotope fractionation leading to an enrichment of heavy isotopes in the residual fraction of the contaminant. CSIA allows interpreting these changes in the isotopic composition of individual organic compounds, thus providing a direct proxy for chemical transformations. CSIA permits qualitative and quantitative information on in situ (bio)degradation, can distinguish different reaction mechanisms, and provides additional means for source allocation and differentiation [START_REF] Schmidt | Compound-specific stable isotope analysis of organic contaminants in natural environments: a critical review of the state of the art, prospects, and future challenges[END_REF], Blessing et al. 2009[START_REF] Braeckevelt | Field applicability of Compound-Specific Isotope Analysis (CSIA) for characterization and quantification of in situ contaminant degradation in aquifers[END_REF][START_REF] Negrel | Isotopic methods give clues about the origins of trace metals and organic pollutants in the environment[END_REF].

We applied compound-specific carbon isotope analysis (CSIA) in a remediation study to evaluate the potential of monitored natural attenuation (MNA) at an industrial site. The study site is an organochlorine production site with a long operational history and leaking facilities or accidental spills released high amounts of contaminants to soil and groundwater. A control and monitoring program of DNAPL at the study site is installed since the 80s. A hydraulical confinement separated the DNAPL source from the contamination plume that consists of lighter, more soluble chlorinated compounds, mainly PCE, TCE, DCE, TeCA, DCA, TCP and DCP. Biostimulation using acetate and H2 amendment was performed in a limited area within the plume.

The aim of the present work is to demonstrate the potential of the CSIA approach for monitoring in-situ biodegradation and for characterizing and quantifying transformation processes of chlorinated pollutants in groundwater towards the ultimate goal of an implementation of natural and/or enhanced bioremediation at the contaminated site. To this end, a total of four monitoring wells in the source zones and nine groundwater wells within the downgradient contaminant plume over several kilometres length were sampled and stable carbon isotope ratios of the aforementioned chlorinated compounds were analyzed. The use of headspace solid-phase micro extraction as sample pre-concentration and extraction technique allowed for compound-specific carbon isotope analysis for most of the compounds even at concentrations in the low ppb-range (≥ 3µg/L). Significant 13 C isotope enrichments with shifts of up to 10‰ were observed along the anoxic flow path of the plume. Observed fractionation is compound-as well as pathway-specific (aerobic and anaerobic biodegradation). Field data evaluated using the Rayleigh equation demonstrate that dispersion plays a major role, esp. within the oxic zones of the plume. For a reliable quantification of the extent of in-situ biodegradation within the aquifer, the selection of appropriate isotope enrichment factors representative of the specific biogeochemical conditions at the site is crucial [START_REF] Pooley | Aerobic Biodegradation of Chlorinated Ethenes in a Fractured Bedrock Aquifer: Quantitative Assessment by Compound-Specific Isotope Analysis (CSIA) and Reactive Transport Modeling[END_REF]. Isotope enrichment factors for the aforementioned organochlorine compounds that have been determined in laboratory microcosm studies under aerobic and anaerobic conditions, respectively, will be presented and compared to the extent of isotope fractionation observed at our field site. In our contribution we will demonstrate the usefulness of the CSIA approach to determine both the nature and the extent of natural attenuation of chlorinated contaminants in the field. In addition, we will present new developments with respect to analytical techniques and demonstrate, based on experience from field studies, how isotope signatures can help revealing different pollution sources and/or spill events.


Control of the Dissolution Rate of an Active Pharmaceutical Ingredient by Using Melt Extrusion Coupled With Supercritical CO 2

INTRODUCTION

Hot-melt extrusion (HME) is a process for converting a raw material into a product of uniform shape and density by forcing it through a die under controlled conditions [START_REF] Rauwendaal | Polymer Extrusion[END_REF]. It represents an efficient manufacturing technology to disperse drugs in a melt up to a true molecular solution of a pharmaceutical ingredient into a polymeric matrix [START_REF] Breitenbach | [END_REF]. Most of the applications describe the preparation of solid dispersions by HME either to increase the aqueous solubility and oral bioavailability of the active substance or to control its release. The advantages of HME are: (i) avoidance of organic solvent; (ii) reduction of processing steps (one pot method); (iii) elimination of the good compressibility requirement for the active ingredients and excipients; (iv) high level of drug content uniformity; and (v) improved bioavailability through drug solubilization or molecular level dispersion in water soluble matrix [START_REF] Repka | Encyclopedia of pharmaceutical technology[END_REF].

One of the major drawbacks of HME is that it can be applied only for producing thermally stable pharmaceuticals. For a number of drugs, proteins and other excipients, HME cannot be used unless a plasticizer is added in order to reduce the viscosity of the mixture lowering thus the processing temperature in the extruder. However, these additives remain in the product and thus alter its properties and performance.

In this context, supercritical carbon dioxide (scCO 2 ) represents an interesting alternative. Indeed, injection of scCO 2 in extrusion process modifies the rheological properties of the polymer in the barrel of the extruder and scCO 2 acts as a blowing agent when flowing through the die 4 . The pressure drop induces a thermodynamic instability in the polymer matrix, generating a large number of bubbles. The cell growth continues until the foam is rigidified (when T<T g ). The reduction of viscosity decreases the mechanical constraints and the operating temperature within the extruder. Thus, combining extrusion with scCO 2 allows using fragile or thermally sensitive molecules, like pharmaceutical molecules, with the absence of residue in the final material.

For the extrusion of pharmaceuticals, copolymers derived from acrylic acid and methacrylates, such as Eudragit types, are promising matrices. Due to their multifunctional nature, unique properties and good biocompatibility, these polymers are considered to be of common grades for coating of pharmaceuticals [START_REF] Dittgen | [END_REF]. In particular, Eudragit E, due to its tertiary amine groups, forms coating soluble in gastric fluid in acidic pH 6 . Eudragit E has already been successfully implemented for HME. Six et al. [7] described the HME of the poorly water soluble itraconazole with Eudragit E100 polymeric carrier. Melt extrusion of this drug with Eudragit E polymer resulted in significantly enhanced dissolution rate. More recently, Verreck et al. showed that pressurized carbon dioxide, acting as a reversible plasticizer and foaming agent, reduces the melt extrusion temperature of several pharmaceutically acceptable polymers, including Eudragit E100 [8].

Carvedilol is a nonselective -adrenergic blocking agent with -blocking activity [9]. It is a lipid soluble compound, practically insoluble in water and poorly absorbed from the gastrointestinal tract. The slow absorption of Carvedilol was attributed to its poor water solubility. Its absolute bioavailability is 25-35%. Eudragit E100 has been already used to improve the bioavailability of Carvedilol. Indeed, E100 nanoparticles containing Carvedilol were prepared by nanoprecipitation [10]. Finally, high encapsulation efficiency was obtained with a faster rate of release of Carvedilol from the nanoparticles.

In this work a scCO 2 coupled extrusion process has been used to generate solid dosage form based on Eudragit E100 with Carvedilol drug according to the method previously elaborated for producing foams of this polymer [11,12].

MATERIALS AND METHODS

Materials

Carvedilol (CAR) was kindly provided by EGIS Plc (Hungary). Eudragit E (E100) was purchased from IMCD (France) and CO 2 from Air Liquide (France).

Extrusion coupled with supercritical CO 2

Figure 1 shows the experimental set up, which has previously been detailed elsewhere [9,10]. The single-screw extruder has a 30 mm-screw diameter and a length to diameter ratio (L/D) of 35 (Rheoscam, SCAMEX). A large L/D ratio generally allows increasing the capacity of mixing and melting but also energy consumption. The first conical part of the screw has a length to diameter ratio of 20 and allows the transport of solid polymers and then, their melting and plasticizing. Then, the screw has a cylindrical geometry from the first gastight ring to the die with a length to diameter ratio of 15. Two restriction rings have been fitted out in order to prevent scCO 2 from backflowing. The diameter of the die was 1 mm and its length was 11.5 mm.

The temperature inside the barrel is measured at five locations: T a and T b before the CO 2 injection, T c and T d after the injection and T e in the die. In this study, the zone temperatures were fixed as follows: T a =125°C, T b =T c =T d =135°C, and T e =130°C. The temperature and the pressure of the polymer inside the extruder were measured by means of three pressure (P 1 , P 2 and P 3 ) and two temperature sensors (T 1 and T 2 ). Errors associated to pressure and temperature sensors were about 0.2 MPa and 3.3°C respectively. CO 2 was pumped from a cylinder by a syringe pump (260D, ISCO) regulated at 7°C. It was then introduced at a constant volumetric flow rate of 0.25 cm 3 .min -1 at a pressure which was kept slightly higher than the pressure P 1 . The CO 2 injector was positioned at a length to diameter ratio of 20 from the feed hopper. The pressure and the volumetric CO 2 flow rate were measured within the syringe pump.

At the beginning of the extrusion process, pure E100 was fed in the extruder. Before the injection of scCO 2 , the torque was controlled by the rotating speed of the screw (higher rotating speed resulted in higher torque) and was kept at high values close to the established alarm limit to ensure high production rate and low residence time. Once steady state conditions have been reached with E100, drug polymer mixture was introduced through the hopper. Two drug mass contents were tested: 5 and 20%. The scCO 2 was introduced into the barrel once the drug polymer mixture had reached the die. During the scCO 2 injection, the torque decreased and the rotation speed could be increased to reach higher production rates and lower residence time. When steady state conditions have been reached again, extrudates were collected and air-cooled at ambient temperature. The samples were obtained under the following operating conditions (Table 1). The sampling time and the mass of the collected samples were measured in order to determine the mass-flow. Extrudates were then ground by an IKA MF10 microfine grinder drive with a MF 10.2 impact grinding head. The applied sieve was MF 2.0 mm and the rotation speed was 3000 rpm.

Scanning electron microscopy (SEM)

Morphology of the samples was investigated by a JEOL 6380LVa type scanning electron microscope. Each specimen was fixed by conductive double sided carbon adhesive tape and gold-coated by JEOL 1200 in order to avoid electrostatic charging.

Fourier transformed infrared spectrometry (FTIR)

KBr pastille sample preparation method was used for the FTIR measurement. The obtained pastilles were analyzed by Bruker Tensor 37 type spectrometer equipped with DTGS detector in the range of 4000-400 cm -1 with a resolution of 4 cm -1 .

Determination of drug content

Samples equivalent to 10 mg of CAR were placed into volumetric flask (1 dm 3 ) and were diluted with 0.1 N HCl. After the complete dissolution and homogenization the concentration of CAR was determined by UV spectrophotometry at 242 nm (Hewlett-Packard HP 8452A).

In vitro dissolution measurement

The dissolution studies were performed by an Erweka DT6 dissolution tester (USP II apparatus). Samples were placed in the dissolution vessel containing 900 ml 0.1 N HCl maintained at 37±0.5°C and stirred at 50 rpm. Samples (5 ml) were collected periodically and the concentration of CAR was determined by UV spectrophotometry at 242 nm (Hewlett-Packard HP 8452A).

RESULTS AND DISCUSSION

Supercritical CO 2 aided extrusion

The application of scCO 2 in the field of pharmaceutical extrusion has two remarkable advantages from the stability point of view. Processing temperature can be decreased [8] and/or the mass flow rate can be increased. If higher mass-flow is achieved, the residence time of the material can be shortened significantly; furthermore, the productivity is higher. A lower temperature and an increased flow rate result in lower risk of decomposition and total amount of impurities (decomposition products) can be kept under the maximum limit.

The aim of this work was to achieve amorphous drug morphology with good homogeneity. Therefore, to be sure that API is melted, the temperature of the process was not decreased (comparing to the reference experiments without scCO 2 ). The melting point of Carvedilol is 117°C [13] while the temperatures of the five zones were kept slightly higher (125-135°C) than this temperature. The residence time, however, was shortened in order to preserve the drug stability. The viscosity of the material in the barrel is correlated with the torque. As the rotation speed was increased, the measured torque increased. Without the injection of scCO 2 , a maximum rotation speed of 5 rpm could be used without alarm stop. In contrast, after introducing scCO 2 , the torque decreased significantly and the rotation speed could be increased stepwise up to 10 rpm. Due to the plasticization effect of scCO 2 , the torque at 10 rpm was lower than at 5 rpm without scCO 2 (Table 1). This gain in the maximal throughput demonstrates the valuable plasticizing efficiency of scCO 2 .

Beside the scCO 2 , the addition of the active ingredient provided remarkable plasticizing effect too, which could be detected by measuring the torque (Table 1). Carvedilol, being a relatively small molecule proved to play the role of a plasticizer. E100-CAR mixtures exhibit lower torque than pure E100 owing to the decreased viscosity. Scanning electron microscopy (SEM)

The surface of extrudates, shown in the SEM images of Figure 3, is quite different to the one made by traditional extrusion process, which is smooth compared to the foamed extrudates. The specific surface area is clearly larger in the case of the foamed samples providing extended availability to the dissolving medium, which may increase the dissolution rate.

(A) (B) Fourier transform infrared spectrometry (FTIR)

The extruded products were analyzed by FTIR spectrometry in order to check their stability under the processing conditions. The results are shown in Figure 4. Compared to the untreated pristine components appearance, no new peak could be detected in the spectra of the extrudates. It confirms that no detectable decomposition occurred during the extrusion. CAR has got two common polymorphic crystalline forms. These are CAR Form I and Form II [14,15]. In the absorption region of O-H and N-H stretching vibration CAR Form II has got one strong band at ~3340 cm -1 and CAR Form I has got two strong bands at ~3250 cm -1 and 3450 cm -1 . E100 has got only two very slight peaks in this region.

In the spectra of E100-CAR extrudates, one broader peak was observable from ~3300 cm -1 to 3450 cm -1 which demonstrates the interaction between the E100 and CAR. It suggests that the CAR is not in crystalline form because the characteristic peaks of Form I and Form II are not present in the spectra. The broad peak indicates probably that CAR is present in amorphous form. This deduction corresponds to the results of Pokharkar et al. [16].

(A) (B) (D) (C)
Further analyzes were made by X-ray diffraction and differential scanning calorimetry (not presented here) and both techniques verified that CAR is in amorphous phase whatever the drug content is (5% or 20%).

In vitro dissolution and drug content determination

The in vitro dissolution of the different samples exhibited significant differences. The dissolution curves are shown in Figure 5. Crystalline CAR (curve E) shows moderate dissolution rate which could be improved by HME with E100 (curves A, B C and D). It was found that the E100 matrix governs the dissolution rate of the extrudates. When molecules of E100 matrix were dissolved, the CAR molecules dissolved simultaneously because CAR molecules were in amorphous form (as shown above) and there was no need to overcome the lattice energy. It confirms that CAR was molecularly dispersed or nearly molecularly dispersed in the E100 matrix.

As E100 dissolves in acidic media rapidly the dissolution times of the extrudates were found to be short as well. The dissolution rate depends on the surface area of the material exposed to the dissolution media [17]. Though the foamed extrudates exhibit an increased surface area, the unmilled foams (curve D) had significantly slower dissolution than common extrudates and ground extrudates. The reason for that result lays on the floating of the unmilled foamy samples on the surface of the dissolution media resulting in limited contact between the phases. In contrast the common and the ground extrudates sank into the dissolution media allowing the dissolution process to proceed faster through the fully accessible surfaces. Consequently the dissolution times of ground samples (5 min) is significantly improved compared to that of crystalline CAR.

(D) (C) (B) (E) (A)

The drug contents of the extrudates, determined with UV spectrophotometer, corresponded in all cases to the initially introduced CAR and the relative standard deviation (1.31 % RSD) of the measured contents was under the established limit of U.S. Pharmacopeia (4.6 % RSD) [START_REF]Uniformity of Dosage Units[END_REF].

CONCLUSION

In this work, fast dissolving formulation of a poorly water soluble drug against hypertension was produced by conventional and supercritical CO 2 assisted extrusion process using Eudragit E as a polymer matrix. The manufactured fine solid dispersions of Carvedilol in Eudragit E increased the dissolution rate of Carvedilol. It shows the dissolution improving efficiency of Eudragit E matrix, produced by HME, and its applicability as fast dissolving drug delivery system. Moreover, the mild operating conditions and the productivity of the technology could be improved by decreasing the viscosity with scCO 2 as a temporary plasticizer.

Figure 1 .

 1 Figure 1. Experimental device

Figure 2 Figure 2 .

 22 Figure2shows the photographs of E100-CAR extrudates. The extrudates made by traditional extrusion process (A) were glassy and transparent. The extrudates made by scCO 2 aided extrusion (B) were opaque due to the newly formed internal surfaces.

Figure 3 .

 3 Figure 3. SEM images of E100-CAR extrudates, (A) exp. 6 , (B) exp. 9

Figure 4 .

 4 Figure 4. FTIR spectra of samples, (A) Exp. 9 ground, (B) E100 powder, (C) CAR Form II, (D) CAR Form I

Figure 5 .

 5 Figure 5. Dissolution profiles, (A) Exp. 9 ground, (B) Exp. 6 ground, (C) Exp. 6, (D) Exp 9, (E) Crystalline CAR

Table 1 . Compositions and operating conditions of the samples

 1

	Exp	CAR (mass %)	CO 2 (cm 3 .min -1 )	Mean torque (Nm)	Screw rotation speed (rpm)
	1	0	0	104	5
	2	5	0	96	5
	3	5	0.25	66	5
	4	5	0.25	75	7
	5	5	0.25	93	10
	6	20	0	83	5
	7	20	0.25	55	5
	8	20	0.25	64	7
	9	20	0.25	76	10