B Polizzi

O Bernard

M Ribot

A time-space model for the growth of micro-algae biofilms for biofuel production

Keywords: Biofilm growth, Photosynthetic micro-algae biofilm, Front propagation, Fluid dynamics model, Numerical simulations

We present in this paper a spatial model describing the growth of a photosynthetic microalgae biofilm. In this model we consider photosynthesis, extracellular matrix excretion, and mortality. These mechanisms are described precisely using kinetic laws that take into account some saturation effects which limit the reaction rates and involve different components that we treat individually. In particular, to obtain a more detailed description of the microalgae growth, we consider separately the lipids they contain and the functional part of microalgae (proteins, RNA, etc ...), the latter playing a leading role in photosynthesis. We also consider the components dissolved in liquid phase as CO 2 . The model is based on mixture theory and the behavior of each component is described on the one hand by mass conservation, which takes into account biological features of the system, and on the other hand by conservation of momentum, which describes the physical properties of the components. Some numerical simulations are displayed in the one-dimensional case and show that the model is able to estimate accurately the biofilm productivity.

Introduction

The term "biofilm" was introduced by Costerton et al. [START_REF] Costerton | How bacteria stick[END_REF] in 1978 to design a mixture of organisms, embedded in an extra-cellular matrix and attached to a surface in contact with water. Biofilm is the prevailing mode of micro-organisms life. They are ubiquitous in nature and they appear commonly in various domains, such as medical infections, reprocessing of waste or production of clean energy.

The formation of a biofilm follows well-identified steps: at the beginning, micro-organisms colonise the surface, then the colony grows and organises within an extra-cellular matrix and finally some of the organisms are dispersed in order to colonise another location of the surface. The extra-cellular matrix plays a particular role in this organisation, since it acts as a barrier and Bastien.Polizzi@unice.fr † olivier.bernard@inria.fr ‡ Corresponding author : magali.ribot@univ-orleans.fr 1 provides a protection for the whole colony. However, the physical structure of a biofilm, which has a gel-like structure, is still unclear: in the same biofilm, some parts have a viscous behavior, whereas other parts have a visco-elastic structure [START_REF] Galy | Mapping of bacterial biofilm local mechanics by magnetic microparticle actuation[END_REF]. Moreover, even if the biofilm contains liquid inside, it is mainly in a solid phase.

During the last decades, the study of micro-algae biofilm has become an important field of research. Although micro-algae biofilms are almost always detrimental in industrial fields [START_REF] Kjellerup | Biocorrosion and biofilm formation in a nutrient limited heating system subjected to alternating microaerophilic conditions[END_REF], it is a credible alternative for the production of clean energy, which is one of the biggest challenge of our century.

Indeed, there exist a wide variety of micro-algae species that are able to develop as a biofilm and to produce chemical components for green oil production. The main interest of a biofilm culture, in comparison with usual cultures of micro-algae in suspension, is the reduction of the harvesting cost. Indeed, it is estimated that the biomass obtained from biofilm contains between 80% and 90% of water instead of more than 99% for cultures in suspension [START_REF] Gross | Development of a rotating algal biofilm growth system for attached microalgae growth with in situ biomass harvest[END_REF]. Moreover, for a classical micro-algae culture in suspension, several steps are needed to remove the water and to concentrate the products, which is both time consuming and expensive. On the opposite, for a biofilm system, micro-algae are stuck on a surface and are harvested by scraping.

Since the comprehension of the development of phototrophic biofilms is an important topic of research, several types of mathematical models have already been proposed. A first kind of models are based on ordinary differential equations (ODE) [START_REF] Wolf | Kinetic modeling of phototrophic biofilms: The PHOBIA model[END_REF] and take into account the kinetic reactions involved in the biofilm growth. In these models, the biological mechanisms are well detailed by considering many components and using a detailed modelling of the mass transfers between components. However, there is no description of the spatial behavior of the biofilm. Another type of models are multidimensional and multi-species cell centered models [START_REF] Van Loosdrecht | Mathematical modelling of biofilm structures[END_REF]. The solutions of these discrete models are obtained by individual based approaches or cellular automata. In these models also, many biological features are included thanks to the precise description of the mechanisms. However, fluid mechanics effects are difficult to take into account and, despite the recent advances in parallel computing that enable to simulate the evolution of a few hundreds of cells, this approach remains quite expensive from a computational perspective. Thus, hybrid models, coupling some cell-centered approaches with partial differential equations, have been proposed to improve this point [START_REF] Alpkvist | Three-dimensional biofilm model with individual cells and continuum EPS matrix[END_REF].

At the same time, fully continuous models have also been developed. In such a model, the biofilm is described as a viscoelastic material that expands in response to the pressure induced by mass exchanges between the biofilm and the surrounding liquid. Two classes of partial differential equations (PDE) models can be distinguished: in the first kind of PDE models, biofilm and liquid are separated by a physical interface, the evolution of which is computed by moving front techniques [START_REF] Alpkvist | A multidimensional multispecies continuum model for heterogeneous biofilm development[END_REF]. On the opposite, other models are based on mixture theory and do not treat the evolution of an interface. Mixture theory was introduced in the 1960s by a series of articles by Trusdell and provides continuous models based on PDEs for multi-component fluids by assuming that several components may be present locally; for more details about mixture theory, see [START_REF] Truesdell | An Introduction to the Mechanics of Fluids[END_REF][START_REF] Rajagopal | Mechanics of mixtures[END_REF]. Mixture models have been successfully used later on for the description of several biological systems and, in particular, for modelling biofilms [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF][START_REF] Clarelli | A fluid dynamics multidimensional model of biofilm growth: stability, influence of environment and sensitivity[END_REF][START_REF] Alpkvist | A multidimensional multispecies continuum model for heterogeneous biofilm development[END_REF][START_REF] Dockery | Finger formation in biofilm layers[END_REF][START_REF] Zhang | Phase field models for biofilms. i. theory and onedimensional simulations[END_REF][START_REF] Zhang | Phase-field models for biofilms ii. 2-D numerical simulations of biofilm-flow interaction[END_REF].

In this article, we focus on the mathematical modelling of micro-algae biofilms through mixture models with a special attention to the biological mechanisms involved in its growth. To do so, we improve and adapt the previous mixture model for biofilms of cyanobacteria in a fountain, proposed by Clarelli et al. in [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF][START_REF] Clarelli | A fluid dynamics multidimensional model of biofilm growth: stability, influence of environment and sensitivity[END_REF]. Note that, from the analytical perspective, the existence of solutions for such mixture models has been established very recently and only for some simplified models [START_REF] Bianchini | Global existence and asymptotic stability of smooth solutions to a fluid dynamics model of biofilms in one space dimension[END_REF].

In our model, the components under consideration are micro-algae, extra-cellular matrix and liquid, on the one hand, and some dissolved components, such as carbon dioxide, oxygen and substrate, on the other hand. In line with quota models usually used to describe growth limited by a nutrient, we consider separately the pool of carbon storage, mainly made of sugars and lipids, and the functional part of micro-algae, composed of proteins, RNA, etc... Therefore, the behavior of each of the seven components considered in the model is described by two conservation laws: mass balance equation, which accounts for mass transfers between the different components and gives the evolution of the component density, and force balance equation, which gives the evolution of momentum and contains physical properties of the phase. The model is supplemented by an incompressibility constraint, obtained from the total mass conservation. This constraint is expressed on the mean hydrodynamic velocity and leads to an hydrostatic pressure term in the force balance equations of the model. As in [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF][START_REF] Clarelli | A fluid dynamics multidimensional model of biofilm growth: stability, influence of environment and sensitivity[END_REF], we keep the inertial terms in momentum equations to guarantee the hyperbolicity of the system and the finite speed of propagation of the front.

However, the adaptation of the model in [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF] to the present context imposes several improvements: first, we consider a larger number of components and therefore we need to work with more velocities. Moreover, we take into account some components dissolved in the water, for which we add a diffusion term in the mass balance equations, leading to reaction-advection-diffusion equations. Finally, the main improvement is the fine description of the mass transfers between components, focusing on the biological mechanisms that drive the development of the biofilm, namely photosynthesis, respiration , functional biomass synthesis, extra-cellular matrix excretion and micro-algae death. To do so, all these mechanisms are expressed under the form of chemical reactions and the kinetic of these reactions is carefully studied. To express the reaction rates, we need to identify clearly the influence of each of the components of the model on the reaction under consideration. Some inhibiting factors or some threshold effects are described by kinetic laws, such as Monod's law [START_REF] Monod | The growth of bacterial cultures[END_REF] or Haldane's law [START_REF] Andrews | A mathematical model for the continuous culture of microorganisms utilizing inhibitory substrates[END_REF]. In particular, we use Droop's model [START_REF] Droop | Vitamin B12 and Marine Ecology. iv. the Kinetics of Uptake, Growth and Inhibition in Monochrysis Lutheri[END_REF] in order to take into account an activation threshold and a saturation effect of an internal ratio for the photosynthesis reaction rate. Finally, combining all this information, we are able to define in an appropriate way the source terms of the mass balance equations. Note that this fine description brings us to consider an incompressibility constraint with a non zero right-hand side, unlike in [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF].

Once the model is established, we propose a numerical scheme in order to perform some accurate numerical simulations, using some finite differences techniques. This numerical scheme is presented in details in the appendices. The basis of the scheme relies on what has been done in [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF], that is to say we compute the hyperbolic part of the problem with a robust Riemann solver based on a relaxation scheme. The same difficulties as in [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF] have to be faced: first, in some regions, some phases may vanish. Unlike gas theory, where we can work in absence of vacuum, the situation when one phase or several phases vanish is relevant in a biological context. In such a case, the computation of the velocities has to be done carefully, using an implicit discretisation in time for the momentum equations. Another common difficulty between the two works is the computation of the pressure, linked to the treatment of the incompressibility constraint. As in [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF], we use a fractional step approach, based on Chorin-Temam projection method [START_REF] Temam | Une méthode d'approximation de la solution des équations de Navier-Stokes[END_REF][START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF], adapted to the particular incompressibility constraint we deal with. However, another difficulty arises in our case: the source terms appearing in the mass balance equations, built from the kinetic laws, may be stiff and lead to numerical instabilities. We therefore use a partially implicit time discretisation for some of the components in the source terms, which ensures that all the densities remain positive.

This numerical scheme enables us to present some numerical simulations in the one-dimensional case in order to study the influence of various parameters of the model, such as the value of the substrate supply or the parameters involved in the light modelling. More precisely, some tests are presented at the end of the article to show the effect of the intensity on the upper surface of the water, the absorption coefficient of light by the biofilm and the day-night variation of light. Some numerical simulations in the two-dimensional case have also been performed, but will be presented in a forthcoming article with the aim of comparing them with experimental results.

The paper is organised as follows: in the first section, we present the components considered in the mathematical model and we detail the biological mechanisms, expressing them under the form of chemical reactions. In the second section, we describe the equations that compose the model, namely the mass balance equations and the force balance equations. A subsection is devoted to the expression of the source terms for the mass balance equations, based on the previously mentioned chemical reactions. To do so, we need to define precisely in a third subsection the form of the reaction rates, taking into account all the factors that may favor or inhibit the reactions. Finally, at the end of this section, we give an estimate of all the parameters included in the equations of the model. We display, in the last section, some numerical results in the one-dimensional case, investigating the influence of several parameters on the model. Finally, conclusion and perspectives are given. Two appendices are dedicated respectively to the description of the numerical scheme and to supplementary numerical simulations.

Description of the biological system

Since the scope of the model is to describe precisely the production of lipids by the micro-algae within the framework of biofilms, we need to consider several mechanisms. In particular, we will study in details the reactions that drive the growth of the biofilm and the production of extra-cellular matrix (ECM). This section is dedicated first to the presentation of the components that are described in our model. Then, we will detail the mechanisms and summarise them through chemical reactions.

Considered components in the mathematical model

In this section, we aim at giving the list of the components quantified by the model, that is to say micro-algae, extra-cellular matrix, liquid and constituents dissolved in liquid phase.

In a biofilm, micro algae are embedded within a self produced extra-cellular matrix. Micro-algae and extra-cellular matrix are both organic compounds but they have some very different roles in the biofilm formation: photosynthesis takes place in the micro-algae, whereas extra-cellular matrix is an inert component. Indeed, extra-cellular matrix is made of dead cells and extra-polymeric substance released by micro-algae such as polysaccharides and small quantities of proteins.

Photosynthesis supports growth of micro-algae and takes place in organelles called chloroplasts, which are part of functional biomass. In order to have a better description of photosynthesis, we track the evolution of the functional biomass present in micro-algae. Therefore, instead of considering micro-algae as a whole, we will consider separately the functional biomass and the pool of carbon storage, following [START_REF] Bernard | A 2D model for hydrodynamics and biology coupling applied to algae growth simulations[END_REF][START_REF] Mairet | Modelling neutral lipid production by the microalga Isochrysis aff. galbana under nitrogen limitation[END_REF]. On the one hand, the functional biomass gathers biosynthetic apparatus such as proteins, nucleic acids and structural material. On the other hand, the pool of carbon storage is defined as the micro-algae without the functional part and is mainly made of sugar, lipids and carbohydrates.

Consequently, the mathematical model under consideration contains four main constituents, that is to say: pool of carbon storage denoted by A, functional biomass denoted by N, extra-cellular matrix denoted by E and liquid phase denoted by L. Note that the two elements A and N form micro-algae, that may be denoted by M when necessary. We also consider three main components dissolved in liquid phase which play an active role in biofilm growth and extra-cellular matrix production: carbon dioxide denoted by C, oxygen denoted by O and substrate (mainly nitrogen) denoted by S.

Before going further, let us detail the definition of A, N and E. Micro-algae are mostly made of water: we consider that 90% of their mass is water. So components A and N represent the pool of carbon storage and the functional biomass within the associated liquid present inside the cellular membrane of micro-algae. For example, the pool of carbon storage gathers sugars, lipids and carbohydrates with a part of the internal water of micro-algae. It is assumed that water inside micro-algae is proportionally divided into A and N according to their mass. Similarly the extra-cellular matrix is a tangle of polymers soaked in water so E denotes the ECM with internal water. Observe that liquid phase L represents the water outside of the biofilm and should not be confused with water trapped into biofilm.

Finally, let us notice that the effect of the light intensity is of course of great importance in order to model photosynthesis; in the following, light intensity is denoted by I. The leading mechanism responsible for the biofilm development is photosynthesis, that is to say the absorption of the carbon dioxide dissolved in water using the light. More precisely, during the photosynthesis process, micro-algae consume carbon dioxide and water, release oxygen and produce carbohydrate in order to gain energy. Inorganic carbon uptake depends on many factors and we summarise in what follows the overriding phenomena that occur during photosynthesis. First, the photosynthesis rate increases with the amount of reactants (carbon dioxide and water), although several studies highlight that this rate saturates when the carbon dioxide concentration exceeds a certain value. Then, the photosynthesis rate also increases with light intensity up to a threshold, above which it decreases due to photo-inhibition [START_REF] Eilers | Dynamic behaviour of a model for photosynthesis and photoinhibition[END_REF][START_REF] Bernard | Hurdles and challenges for modelling and control of microalgae for CO 2 mitigation and biofuel production[END_REF]; since biofilm develops during several days, it also benefits from a periodical enlightment, following the day and night cycle. Another inhibiting factor for photosynthesis is the oxygen concentration [START_REF] Costache | Comprehensive model of microalgae photosynthesis rate as a function of culture conditions in photobioreactors[END_REF]. Finally, as explained in previous subsection, photosynthesis happens within the functional biomass and, therefore, the quota of functional biomass present in the micro-algae will affect the photosynthesis rate.

Main biological mechanisms considered

Moreover, when the micro-algae are not illuminated, they use another source of energy provided by the respiration. Basically, it is the reverse of photosynthesis: the cells use the oxygen dissolved in water to oxidise the carbohydrate contained in the pool of carbon storage of micro-algae. During this process, carbon dioxide and water are released.

Finally, since we quantify the production of functional biomass, we consider that this process consumes the carbohydrate contained inside micro-algae and the nitrate dissolved in the water.

Extra cellular matrix formation

Let us describe the different mechanisms that lead to the formation of extra-cellular matrix, which is composed both of dead cells and extra-polymeric substance (EPS) released by micro-algae.

First, part of the extra-cellular matrix is composed of EPS excreted by the micro-algae. During the photosynthesis process, when micro-algae are illuminated in the presence of carbon dioxide, they store carbon under the form of carbohydrate. Then, in the presence of nitrogen, they produce functional biomass and divide. However, if no nitrogen is available and the pool of carbon storage is full, micro-algae start to release polysaccharides [START_REF] Staats | Exopolysaccharide production by the epipelic diatom Cylindrotheca closterium: effects of nutrient conditions[END_REF], but also some parts of the functional biomass, composed of proteins.

The other mechanism which produces extra-cellular matrix is the death of micro-algae which occurs in absence of light and nutrients. We assume that micro-algae death rate depends exclusively on the dissolved oxygen concentration: the micro-algae start to die when the oxygen concentration becomes smaller than a threshold. On the opposite, the dissolved oxygen concentration can exceed more than three times the oxygen saturation in liquid at equilibrium, due to photosynthesis. In such a case, the medium is unbalanced which also causes the death of micro-algae. In practice, when the micro-algae die, the functional biomass and the pool of carbon storage contained in the cells are released and become entangled within the surrounding extra-cellular matrix.

Representations of the biological mechanisms under the form of chemical reactions

The aim of this subsection is to give a schematic representation of five of the mechanisms considered in Subsections 2.2.1 and 2.2.2 under the form of chemical reactions, involving the components listed in Sec. 2.1. More precisely, we restrict ourselves to the description of photosynthesis, respiration, functional biomass synthesis, EPS excretion and micro-algae death. All these chemical reactions depend on the stoichiometry and the reaction rates.

The stoichiometric coefficient of a chemical reaction is equal to the ratio between the number of units of reactants consumed and the number of units of synthesised products. However, some components of the model stand for several chemical elements and some reactions are not perfectly balanced. As a consequence, we would rather use pseudo-stoichiometric coefficients to take into account losses that may occur during the reactions. In the following, pseudo-stoichiometric coefficients are denoted by

η φ P (resp. η φ R , η φ N , η φ E and η φ D )
for Photosynthesis (resp. Respiration, functional biomass synthesis, EPS excretion and micro-algae Death) where φ stands for one of the components φ ∈ {A, N , E, L, S, C, O}. Note that the values of all these parameters are fixed later on in Table 1.

Moreover, the reaction rate describes the number of units of the main synthesised product by unit of time. This function depends on several parameters, in particular the concentrations of reactants. In the following, we denote by ϕ P (resp. ϕ R , ϕ N , ϕ E , ϕ D ) the reaction rate for photosynthesis (resp. respiration, functional biomass synthesis, EPS excretion and micro-algae death) and, in the following subsections, we will make precise for each reaction with which component we normalise the reaction rate. All these reaction rates will be modeled in details as source terms of the mass balance equations later on, in Subsection 3.2.2.

Photosynthesis

As described in Sec. 2.2.1, during photosynthesis, micro-algae assimilate the carbon dioxide contained in liquid phase thanks to light. More precisely in the presence of light, micro-algae consume carbon dioxide and liquid to produce sugar and reject oxygen. In terms of the components of Sec. 2.1, photosynthesis consumes C and L to produce A and O. Therefore, the schematic representation of photosynthesis is

η C P C + η L P L ϕ P ----------→ A + η O P O, (2.1)
where ϕ P is the reaction rate of photosynthesis normalised with respect to the pool of carbon storage A and η φ P for φ = C, L, O are the pseudo-stoichiometric coefficients of photosynthesis.

Respiration

In absence of light, respiration provides energy to the micro-algae cells by breaking large molecules, such as carbohydrate, contained in the pool of carbon storage into smaller ones. The respiration process consumes carbohydrate contained into micro-algae and oxygen dissolved in liquid to produce energy; this mechanism releases carbon dioxide and water. It can therefore be written under the following form:

A + η O R O ϕ R ----------→ η C R C + η L R L, (2.2)
where ϕ R is the reaction rate of respiration normalised with respect to A and η φ R for φ = C, L, O are the pseudo-stoichiometric coefficients of respiration.

Functional biomass synthesis

Functional biomass is composed of nitrogen embedded into long carbon chains. Carbon chains are provided by the pool of carbon storage A, whereas nitrogen comes from the substrate S. Therefore, the reaction is represented by

η A N A + η S N S ϕ N ----------→ N, (2.3)
where ϕ N is the reaction rate, normalised by the functional biomass N and η A N and η S N are the pseudo-stoichiometric coefficients for A and S of the synthesis of functional biomass.

Extra-polymeric substance (EPS) excretion

As explained in Sec. 2.1, extra-cellular matrix E constitutes a fundamental component of the biofilm and is mainly made of polysaccharides and of a few proteins. The polysaccharides released by micro-algae come from the pool of carbon storage A, whereas proteins are produced by the functional biomass N. Therefore EPS (extra-polymeric substrate) excretion is divided into two reactions:

A ϕ A E ----------→ E, (2.4a
)

N ϕ N E ----------→ E, (2.4b)
where ϕ A E (resp. ϕ N E ) is the reaction rate for the excretion of polysaccharides (resp. of functional biomass).

Micro-algae death

Due to environmental factors, micro-algae may die and the content of the cells goes into the extra-cellular matrix. Consequently, the death can be represented by two reactions, one for the pool of carbon storage and the other one for the functional biomass, namely:

A ϕ A D ----------→ E, (2.5a)
N ϕ N D ----------→ E, (2.5b)
where ϕ A D and ϕ N D are the two reaction rates.

Description of the mathematical model

In this section, we introduce a mathematical model describing the formation and the growth of a biofilm of micro-algae together with the biofilm composition. All the phenomena mentioned in previous section will be integrated into the model. To avoid the mathematical complexity associated with the introduction of a physical interface between the biofilm and the external medium, we use the mixture theory framework, following [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF]. This model is composed of a mass balance equation for each component coupled with force balance equations in which the physical properties of each phase are described. Since the mass balance equations are written in terms of volume fractions, the model includes a volume constraint, which can be expressed equivalently as an incompressibility constraint, in the spirit of Navier-Stokes incompressibility constraint. Notice that, since we are dealing with bio-chemical reactions, the expressions of the source terms for the mass balance equations are complex nonlinear functions which will be detailed in section 3.2.

Mass balance equations

As mentioned in Section 2.1, we split the components considered in this model in two parts: on the one hand, the ones that constitute biofilm and liquid, and on the other hand, the ones dissolved in liquid.

Notations and conventions

In the following, we denote by t 0 the time variable and x ∈ Ω the position with Ω the spatial domain under consideration which may be a subset of R, R 2 or R 3 . Note that, although this article focuses on the one-dimensional case, our aim is to perform numerical simulations and comparisons with experiments in a 2D or 3D setting. Therefore, we present the following mathematical model in a general context.

For a component, we denote by φ its volume fraction, by ρ φ the volumetric mass density, by v φ its velocity and by Γ φ the source term of the mass balance equation describing the evolution of its volume fraction. In particular, we denote by A (resp. N , E, L) the volume fraction of pool of carbon storage A (resp. of functional biomass N, extra-cellular matrix E and liquid L), namely the volume occupied by the component in an elementary volume V divided by the volume of V.

Remark that pool of carbon storage A and functional biomass N are both parts of micro-algae; therefore, we use for these two elements the same density, denoted by ρ M and the same velocity called v M , where M stands for the whole micro-algae.

Mass balance equations for biofilm and liquid

By a balance of mass in the elementary volume V, we find that the evolution with respect to time of the mass of one of the components verifies

d dt V ρ φ φ dV = - ∂V ρ φ φv φ • ndσ + V Γ φ dV ,
where ∂V is the boundary of V and n its unitary external normal. We can thus deduce the following mass balance equation:

∂ t ρ φ φ + ∇ X • ρ φ φv φ = Γ φ .
Assuming that this component is incompressible, that is to say that ρ φ is constant in time and space, we finally obtain the following equality:

∂ t φ + ∇ X • φv φ = Γ φ ρ φ .
Thus, for φ = A, N , E and L, the corresponding mass balance equations write as follows:

∂ t A + ∇ X • (Av M ) = Γ A ρ M (3.1a) ∂ t N + ∇ X • (N v M ) = Γ N ρ M (3.1b) ∂ t E + ∇ X • (Ev E ) = Γ E ρ E (3.1c) ∂ t L + ∇ X • (Lv L ) = Γ L ρ L .
(3.1d)

Incompressibility constraint

Assuming that the volume of the dissolved components is negligible, for any time t 0 and position X ∈ Ω, the volume fractions A, N , E and L satisfy the volume condition:

A(t, X) + N (t, X) + E(t, X) + L(t, X) = 1. (3.2)
Consequently, by summing the four equations (3.1) and using condition (3.2), we obtain an incompressibility constraint:

∇ X • ((A + N )v M + Ev E + Lv L ) = Γ A + Γ N ρ M + Γ E ρ E + Γ L ρ L . (3.3)
Let us notice that we can work with the two equations (3.2) and (3.3) in the following, omitting therefore the mass balance equation for liquid (3.1d).

Mass balance equations for components dissolved in liquid

Let us now consider the components dissolved in liquid phase, that is to say substrate S, carbon dioxide C and oxygen O. Since their volume fraction is negligible, we consider their mass fraction, denoted by S, C and O, that is to say the mass contained in liquid divided by the mass of liquid.

Their volumetric mass density is the same as the liquid density ρ L and their velocity is driven by the velocity of liquid phase v L . Moreover, in addition to the convection phenomenon, we take into account the diffusion and we add a diffusive term in the mass balance equations. Following [START_REF] Ambrosi | On the closure of mass balance models for tumor growth[END_REF] we assume that the diffusion follows the Fick's law, namely the diffusive flux is assumed to be proportional to the concentration gradient of dissolved component in liquid. The concentration of dissolved components in liquid phase contained in the elementary volume V is φρ L for φ = S, C and O. So the integral form of the diffusion terms writes as

∂V Lδ φ ∇ X (φρ L ) dσ = V ∇ X • Lδ φ ∇ X (φρ L ) dV
where δ φ are the diffusion coefficients. Therefore the local formulation of the mass balance equation for the components dissolved in liquid phase writes as

∂ t (SL) + ∇ X • (SLv L ) -∇ X • (δ S L∇ X S) = Γ S ρ L (3.4a) ∂ t (CL) + ∇ X • (CLv L ) -∇ X • (δ C L∇ X C) = Γ C ρ L (3.4b) ∂ t (OL) + ∇ X • (OLv L ) -∇ X • (δ O L∇ X O) = Γ O ρ L . (3.4c)

Source terms and reaction rates

In this section, we will make precise the form of the source terms that appear on the right-hand side of equations (3.1) and (3.4) and that describe the mechanisms presented in Section 2.3. We will detail this expression for each component in the following subsection and then we will give the expressions of the reaction rates that may depend on several factors.

Source terms

Pool of carbon storage A is a product of photosynthesis, see Sec. By combining all these terms, we obtain that the source term for A can be written as:

Γ A = ϕ P -ϕ R -η A N ϕ N -ϕ A E -ϕ A D . (3.5a)
Now, considering Sec. 2.3.3, functional biomass N is produced with rate ϕ N . However, the loss of functional biomass is due to EPS (extra-polymeric substance) excretion with rate ϕ N E and to mortality with rate ϕ N D , see equations (2.4b) and (2.5b). Therefore, we can define the source term for N as:

Γ N = ϕ N -ϕ N E -ϕ N D . (3.5b)
Extra-cellular matrix E is only produced and comes from two mechanisms: EPS excretion and mortality, which both can use either the pool of carbon storage or the functional biomass as a reactant. In all cases, according to Sec. 2.4a, 2.4b, 2.5a and 2.5b, pseudo stoichiometric coefficients are equal to 1 and we obtain the following source term:

Γ E = ϕ A E + ϕ N E + ϕ A D + ϕ N D . (3.5c)
Water which is the main component of liquid phase L is a reactant in photosynthesis with rate η L P ϕ P and a product in respiration with rate η L R ϕ R , see Sec. 2.1 and 2.2, so:

Γ L = η L R ϕ R -η L P ϕ P . (3.5d)
Now, substrate S is only consumed by functional biomass synthesis (2.3), with rate η S N ϕ N ; thus, we define the source term for S as:

Γ S = -η S N ϕ N . (3.5e)
Carbon dioxide C, as liquid, is a reactant for photosynthesis (2.1) and a product of the respiration process (2.2); therefore we set:

Γ C = η C R ϕ R -η C P ϕ P . (3.5f)
Finally, oxygen O is released by the photosynthesis mechanism (2.1) with rate η O P ϕ P and consumed by the respiration mechanism with rate η O R ϕ R . Therefore, the source term for C can be written as:

Γ O = η O P ϕ P -η O R ϕ R .
(3.5g)

Expressions of reaction rates

Now, in this section, let us give the expressions of the reaction rates presented in Section 2.3 and used in previous subsection to model mass exchanges between species. To do so, we will use classical kinetics models, like Michaëlis-Menten's law for enzyme catalyzed reactions or Haldane's kinetics law when inhibition is considered. These laws are defined as dimensionless elementary functions with values in the interval [0, 1] and each reaction rate is built as a product of elementary functions multiplied by the maximal rate per unit of time. In what follows, this maximal rate is called µ α where α represents the considered mechanism; the estimated values of these rates are given in Table 1, as well as the other constants of the model. Note that some mechanisms may depend on the intra-cellular quota Q defined by

Q = N N + A ,
which represents the amount of functional biomass per unit of micro algae biomass. In particular, this quota will be used to express the photosynthesis rate. Photosynthesis rate. Photosynthesis 2.1, which is the most complex mechanism since its reaction rate is driven by all components except substrate, is also obviously influenced by the light intensity, denoted by I. For the reader's convenience, we consider separately the effects of each component and we denote by f φ , with φ = Q, N , L, C, O, I, the corresponding elementary function which depends only on φ and takes values in [0, 1], except f N . The expression for the photosynthesis rate will finally be the product

ϕ P = µ P f N f Q f L f C f O f I , (3.6)
where µ P is the maximal photosynthesis rate by unit of micro-algae and unit of time.

First, we assume that the growth of the pool of carbon storage depends linearly on the functional biomass amount, so f N is equal to ρ M N up to a constant k, to be determined.

Then, µ P is defined as the maximal photosynthesis rate by unit of mass of micro-algae and unit of time, that is to say as the maximal value for 1 ρ M ϕ P N +A . But the maximal value for 1

ρ M ϕ P
N +A is also equal to µ P k N N +A max = µ P kQ max according to eq. (3.6).

Therefore k = 1 Q max
and we set

f N = ρ M N Q max (3.7a)
with Q max a threshold above which the functional biomass production is stopped, that is to say the quota Q does not have any impact on the photosynthesis rate if Q ≥ Q max . Moreover, as mentioned in Sec. 2.2.1, photosynthesis is performed by chloroplasts, contained in the functional biomass. Their quantity increases with the intra-cellular quota and photosynthesis occurs when Q exceeds a given threshold Q min , which can be described by Droop's model [START_REF] Bernard | Transient behavior of biological loop models with application to the Droop model[END_REF][START_REF] Droop | Vitamin B12 and Marine Ecology. iv. the Kinetics of Uptake, Growth and Inhibition in Monochrysis Lutheri[END_REF]. We obtain therefore:

f Q = Q max Q max -Q min max 0, 1 - Q min min {Q, Q max } , (3.7b) which implies that f Q = 0 for Q < Q min and f Q = 1 when Q > Q max .
Now, to take into account the impact of reactants, namely carbon dioxide and liquid, we use Michaëlis-Menten's law which reflects an almost linear increase for a low amount of reactants and a threshold effect when reactants are in excess. The photosynthesis rate depends consequently on

f C = C K C + C , f L = L K L + L , ( 3.7c)
with K C and K L the half saturation constants of carbon dioxide and liquid. So we have

f φ ∼ φ K φ
when φ goes to 0 and f φ ∼ 1 when φ goes to infinity, for φ = C and L. Moreover, according to [START_REF] Costache | Comprehensive model of microalgae photosynthesis rate as a function of culture conditions in photobioreactors[END_REF], photosynthesis is inhibited when the oxygen concentration is too high. We model this effect with the following sigmoid function:

f O = 1 1 + O K O α (3.7d)
with K O the half saturation constant and α > 0 a strictly positive parameter that determines the stiffness of the sigmoid. Finally, we need to consider a last effect, which is the impact of light intensity. Since the reaction is inhibited when light intensity is too high, we model the effect of light on the photosynthesis rate by Haldane's law [START_REF] Eilers | Dynamic behaviour of a model for photosynthesis and photoinhibition[END_REF][START_REF] Thébault | Comparison between two mathematical formulations of the phytoplankton specific growth rate as a function of light and temperature, in two simulation models[END_REF], that is to say:

f I = 2 (1 + K I ) Î Î2 + 2K I Î + 1 , where Î = I I opt , ( 3.7e)
with I the received light intensity, I opt the optimal light intensity and K I a parameter. This function behaves almost like a linear function of Î for Î < 1 and like 1 Î for Î greater than 1.

The received light intensity I is a space-dependent parameter, since light is absorbed according to the various layers it encounters. So local light intensity I(t, x, y, z) depends on the incident light intensity on the upper surface of water I 0 (t, x, y) and on an attenuation coefficient, which evolves with the depth and the composition of the above medium. Precisely, the attenuation coefficient at position X = (x, y, z) depends on the depth and on the medium composition at points (x, y, ξ) for z ≤ ξ ≤ H, with H the height of the domain.

Following [START_REF] Stomp | Colorful niches of phototrophic microorganisms shaped by vibrations of the water molecule[END_REF], we can describe the received light intensity by:

I(t, X) = I 0 (t, x, y) exp - H z τ L L (t, x, y, ξ) + τ M (A (t, x, y, ξ) + N (t, x, y, ξ) + E (t, x, y, ξ)) dξ , (3.8)
where τ L and τ M are the absorption coefficients of liquid and biofilm, estimated in Table 1 Finally, by combining equations (3.7a) to (3.7e), we deduce the form of the photosynthesis rate, that is to say:

ϕ P = µ P ρ M Q max -Q min N C K C + C L K L + L max 0, 1 - Q min min {Q, Q max } 1 1 + O K O α 2 (1 + K I ) Î Î2 + 2K I Î + 1 . (3.9a)
Respiration rate. For the respiration, we only take into account the effects of the reactants, namely the oxygen and the pool of carbon storage, see Sec. 2.3.2. We assume that ϕ R depends linearly on the quantity of the pool of carbon storage ρ M A and that the Michaëlis-Menten's law applies for the dependence on oxygen; we denote by K R the half saturation constant for the oxygen and by µ R the maximal respiration rate. Thus, we obtain the following form for the respiration rate:

ϕ R = µ R ρ M A O K R + O (3.9b)
which induces a quasi linear increase for small values of O and a threshold effect for O K R .

Functional biomass synthesis rate. The functional biomass is synthesised thanks to a pool of enzymes and other molecules containing nitrogen such as DNA, RNA. Thus ϕ N (see 2.3.3) is assumed to depend linearly on the quantity of functional biomass ρ M N . It also depends on the availability of the two reactants, that is to say the substrate and the pool of carbon storage, and on the need for functional biomass. The influence of the substrate S can be modelled by Michaëlis-Menten's kinetic with the half saturation constant K S . In line with Droop's model, we assume that the growth process depends on the fraction of nitrogen per unit of biomass. The availability in carbohydrate contained in the pool of carbon storage is related to the intra-cellular quota: a low value of Q means high disposability whereas large values of Q indicate a low availability of A. Similarly, since low values for the intra-cellular quota indicate a lack of functional biomass, the production of functional biomass decreases with the intra-cellular quota. As a consequence the availability in carbohydrate and the need for functional biomass are modelled together by a linear decreasing function of the intra-cellular quota which vanishes for Q Q max . Therefore, we set

ϕ N = µ N ρ M N S K S + S max 0, Q max -max {Q, Q min } Q max -Q min , (3.9c)
with µ N the maximal rate of functional biomass synthesis. Note that ϕ N = 1 for Q smaller than Q min and ϕ N = 0 when Q is greater than Q max . EPS excretion rates. The EPS is mainly excreted from the pool of carbon storage and, in smaller amounts, from the functional biomass. Therefore, the EPS excretion rate ϕ A E (resp. ϕ N E ) is proportional to the quantity of the pool of carbon storage (resp. functional biomass) contained in micro-algae. In both cases, the mechanism is the same but the kinetic of excretion is slower for the EPS coming from functional biomass. So, the maximal excretion rate µ N E for the EPS coming from N in the reaction (2.4b) is smaller than the maximal excretion rate µ A E for EPS coming from A in reaction (2.4a) -their values are estimated in Table 1. Moreover, we assume that the EPS excretion rate depends on the intra-cellular quota Q similarly as for the functional biomass synthesis rate. Combining all these effects, we obtain:

ϕ A E = µ A E ρ M A max 0, Q max -max {Q, Q min } Q max -Q min , (3.9d) ϕ N E = µ N E ρ M N max 0, Q max -max {Q, Q min } Q max -Q min . (3.9e)
Micro-algae death rates. Finally, the micro-algae death rates ϕ A D and ϕ N D introduced in Sec. 2.3.5 are two parts of the same mechanism. Therefore, their expressions are identical except that ϕ A D is proportional to the quantity of the pool of carbon storage and ϕ N D to the quantity of functional biomass. The variations of the death rates are driven essentially by the dissolved oxygen concentration: if the oxygen concentration is below the reference concentration K D , the mortality of micro-algae increases. On the opposite, due to photosynthesis, the oxygen concentration can exceed the reference concentration , which also increases the death rate. Consequently the death rate depends on the ratio between the oxygen concentration and the optimal oxygen concentration:

Ô = O K D
. In order to model a quasi linear decrease of the death for O < K D we use a modified Michaëlis-Menten's law:

f D Ô = 1 - β Ô Ôβ + β -1 (3.9f)
with β > 1 a parameter estimated in Table 1, that settles the stiffness of the increase of the death rate for O > K D . Observe that we put 1β in the denominator in order that f D remains in [0, 1] for Ô 0 and f D (1) = 0. Finally we set:

ϕ A D = µ D ρ M A 1 - β Ô Ôβ + β -1 , (3.9g) ϕ N D = µ D ρ M N 1 - β Ô Ôβ + β -1 , (3.9h)
where µ D is the maximal mortality rate.

Force balance equations

In this section, we establish some equations to compute the velocities used in the mass balance equations. To do so, we write some evolution equations for the momentum thanks to some force balance considerations. Following what we have done in Sec. 3.1, we consider the three following velocities: velocity of micro-algae v M , velocity of extra-cellular matrix v E and liquid velocity v L . Keeping the notations of Sec. 3.1, for φ = M, E or L, the evolution in time of momentum in an elementary volume V satisfies

∂ t V ρ φ φv φ dV = F vol + F surf , (3.10)
where F vol are the forces inside the volume V and F surf are the forces acting through the surface of V. Forces inside the volume can be decomposed as body forces, contact forces between species and source of momentum resulting from mass exchanges, leading to the following expression:

F vol = V ρ φ b φ dV + V mφ dV + V Γ φ v φ dV . (3.11)
Moreover, forces acting through the boundary come from the flux of momentum and from contact forces within the component, that is to say:

F surf = - ∂V ρ φ φv φ v φ • n dσ + ∂V Tφ • ndσ . (3.12)
Finally, we obtain the following formulation for the force balance equation:

∂ t ρ φ φv φ + ∇ X • ρ φ φv φ ⊗ v φ = ∇ X • Tφ + ρ φ b φ + mφ + Γ φ v φ . (3.13)
Since the mixture is saturated, i.e. without vacuum, following [START_REF] Ambrosi | On the closure of mass balance models for tumor growth[END_REF][START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF], we can decompose Tφ = -φP I + φT φ and mφ = P ∇ x φ + m φ , where P is the hydrostatic pressure, T φ is the excess stress tensor, m φ is the excess interaction force and I is the identity matrix. In line with previous works [START_REF] Clarelli | Mathematical models for biofilms on the surface of monuments[END_REF][START_REF] Zhang | Phase field models for biofilms. i. theory and onedimensional simulations[END_REF][START_REF] Zhang | Phase-field models for biofilms ii. 2-D numerical simulations of biofilm-flow interaction[END_REF], we assume that the body forces are negligible. So, equation (3.13) can be written as:

∂ t ρ φ φv φ + ∇ X • ρ φ φv φ ⊗ v φ = m φ -φ∇ x P + ∇ X • φT φ + Γ φ v φ . (3.14)
In a close mixture system, the sum of the momentum supply due to the interaction forces and of the momentum transfers caused by the mass exchanges is equal to zero, that is to say φ=A,N ,E,L mφ + Γ φ v φ = 0. Then, using the definition of mφ and condition (3.2), we deduce that we have the following relation for liquid phase:

m L + Γ L v L = -(m M + (Γ A + Γ N )v M + m E + Γ E v E ) . ( 3.15)
Moreover, following [START_REF] Ambrosi | On the closure of mass balance models for tumor growth[END_REF][START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF], the general form of the excess stress tensor for an elastic fluid is given by

φT φ = -Σ φ (φ) -φκ φ (φ)∇ X • v φ I + φµ φ (φ) ∇ x v φ + ∇ x v φ T ,
with Σ φ (φ) the elastic interaction and µ φ (φ) the viscosity coefficient. We assume that the stress is constant and we neglect all the shear stress effects for the biofilm components A, N and E, namely κ φ (φ) = µ φ (φ) = 0 for φ = A, N , E. Furthermore, we consider that Σ φ depends linearly on φ: Σ φ (φ) = γ φ φ with γ φ > 0 such that for φ = A, N , E, equation (3.14) becomes

∂ t ρ φ φv φ + ∇ X • ρ φ φv φ ⊗ v φ = -φ∇ x P -∇ X • γ φ φ + m φ + Γ φ v φ . (3.16)
Regarding liquid phase, we make the classical assumption that T L = 0 meaning that the excess stress tensor is present only in solid components [START_REF] Rajagopal | Mechanics of mixtures[END_REF]. Thus, in liquid phase, the only external force is the hydrostatic pressure in order to have liquid at rest in absence of biofilm. Consequently, using eq. (3.15) to simplify eq. (3.14) for φ = L, we deduce that the force balance equation for liquid phase can be written as:

∂ t (ρ L Lv L ) + ∇ X • (ρ L Lv L ⊗ v L ) = -L∇ x P -m M -m E -Γ M v M -Γ E v E .
(3.17)

Finally, let us detail the interaction forces. Following [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF], we assume that the contact forces obey Darcy's law and are proportional to the relative difference between the velocities of the components. Thus, we set

m M = -m ML (v M -v L ) -m ME (v M -v E )
and

m E = -m EL (v E -v L ) -m EM (v E -v M ) ,
where the friction coefficients m ML , m ME , m EL and m EM are experimental parameters. In addition, we assume that the friction force of extra-cellular matrix on biofilm and the friction force of biofilm on extra-cellular matrix have the same magnitude, which leads to m EM = m ME . These parameters are estimated in Table 1.

To sum up, from equations (3.16) and (3.17), we deduce that force balance equations are equal to:

∂ t ((A + N )v M ) + ∇ X • ((A + N )v M ⊗ v M ) = 1 ρ M -(A + N )∇ x P -∇ X • γ M (A + N ) (3.18a) -m ML (v M -v L ) -m ME (v M -v E ) + Γ M v M ∂ t (Ev E ) + ∇ X • (Ev E ⊗ v E ) = 1 ρ E (-E∇ x P -∇ X • (γ E E) -m EL (v E -v L ) + m ME (v M -v E ) + Γ E v E ) (3.18b) ∂ t (Lv L ) + ∇ X • (Lv L ⊗ v L ) = 1 ρ L (-L∇ x P + m ML (v M -v L ) + m EL (v E -v L ) -Γ M v M -Γ E v E ) . (3.18c)

Boundary conditions

We complement this system with some boundary conditions as follows: we impose Neumann boundary conditions for the components constituting biofilm and for the liquid. For components dissolved in liquid phase, we also impose Neumann boundary conditions on all boundaries except on the boundary on top, where we impose some non homogeneous Dirichlet boundary conditions. These boundary conditions stand for external supply of substrate coming from the top of the domain.

The definition of the boundary conditions for the velocities requires more attention. Indeed by integrating the incompressibility constraint (3.3) on the whole domain Ω and using Stokes' theorem, we get the following compatibility condition

∂Ω ((A + N )v M + Ev E + Lv L ) • n dν = Ω Γ A + Γ N ρ M + Γ E ρ E + Γ L ρ L dω.
In order to satisfy this compatibility condition, we impose for all the velocities some no flux boundary conditions on each boundary except on the top of the domain where the flux is imposed so that the compatibility condition holds.

Complete system of PDEs in 1D

In this subsection, let us summarise the full set of equations we derived in previous subsections.

On the one dimensional domain Ω = [0, L x ] the full system consists in mass balance equations (3.1) -(3.4) coupled with source terms (3.5) and reaction rates (3.9), force balance equations (3.18), condition (3.2) and the incompressibility constraint (3.3), complemented by boundary conditions:

Reaction rates

ϕ P = µ P ρ M Q max -Q min N C K C + C L K L + L 2 (1 + K I ) Î Î2 + 2K I Î + 1 max 0, 1 - Q min min {Q, Q max } 1 1 + O K O α , ϕ R = µ R ρ M A O K R + O , ϕ N = µ N ρ M N S K S + S max 0, Q max -max {Q, Q min } Q max -Q min , ϕ A E = µ A E ρ M A max 0, Q max -max {Q, Q min } Q max -Q min , ϕ N E = µ N E ρ M N max 0, Q max -max {Q, Q min } Q max -Q min , ϕ A D = µ D ρ M A 1 - β Ô Ôβ + β -1 , ϕ N D = µ D ρ M N 1 - β Ô Ôβ + β -1 , Q = N N + A , Ô = O K D , Î = I I opt , I(t, x) = I 0 (t) exp - L x x τ L L (t, ξ) + τ M A (t, ξ) + N (t, ξ) + E (t, ξ) dξ
Source terms

Γ A = ϕ P -ϕ R -η A N ϕ N -ϕ A E -ϕ A D , Γ N = ϕ N -ϕ N E -ϕ N D , Γ E = ϕ A E + ϕ N E + ϕ A D + ϕ N D , Γ L = η L R ϕ R -η L P ϕ P , Γ S = -η S N ϕ N , Γ C = η C R ϕ R -η C P ϕ P , Γ O = η O P ϕ P -η O R ϕ R ,

Mass balance equations

∂ t A + ∂ x (Av M ) = Γ A ρ M , ∂ t N + ∂ x (N v M ) = Γ N ρ M , ∂ t E + ∂ x (Ev E ) = Γ E ρ E ,

Volume constraint and incompressibility constraint

A + N + E + L = 1, ∂ x ((A + N )v M + Ev E + Lv L ) = Γ A + Γ N ρ M + Γ E ρ E + Γ L ρ L ,
Mass balance equations for components dissolved in liquid

∂ t (SL) + ∂ x (SLv L ) -∂ x (δ S L∂ x (S)) = Γ S ρ L , ∂ t (CL) + ∂ x (CLv L ) -∂ x (δ C L∂ x (C)) = Γ C ρ L , ∂ t (OL) + ∂ x (OLv L ) -∂ x (δ O L∂ x (O)) = Γ O ρ L ,
Force balance equations

∂ t ((A + N )v M ) + ∂ x (A + N )v M 2 = 1 ρ M -(A + N )∂ x P -∂ x γ M (A + N ) -m ML (v M -v L ) -m ME (v M -v E ) + (Γ A + Γ N )v M , ∂ t (Ev E ) + ∂ x Ev E 2 = 1 ρ E -E∂ x P -∂ x (γ E E) -m EL (v E -v L ) + m ME (v M -v E ) + Γ E v E , ∂ t (Lv L ) + ∂ x Lv L 2 = 1 ρ L -L∂ x P + m ML (v M -v L ) + m EL (v E -v L ) -(Γ A + Γ N )v M -Γ E v E .
Boundary conditions

∂ x φ (t, 0) = 0, ∂ x φ (t, L x ) = 0, for φ = A, N , E, L, v φ (t, 0) = 0, v φ (t, L x ) = L x 0 Γ A + Γ N ρ M + Γ E ρ E + Γ L ρ L dx, for φ = M, E, L, ∂ x (φL) (t, 0) = 0, ∂ x (φL) (t, L x ) = θ φ , for φ = S, C, O,
where θ S , θ C , θ O are the external supplies of substrate, carbon and oxygen.

Parameters estimates

The parameter values chosen for the simulations are presented in Table 1. Most of the parameters are derived from the literature and some parameters are derived from the state of the art on photosynthesis. The half saturation constants used in the reaction rates are generally expressed in terms of concentration in kg L -1 or in mol L -1 . However, in our case, the amount of each component is dimension free and the values given in Table 1 are therefore converted, using the volumetric mass density or the molar mass.

For the photosynthesis and the respiration processes, the pseudo-stoichiometric coefficients are estimated using the molar mass of the components involved in the reactions and considering that, during the reaction, around 90% of the exchanged matter is water.

Moreover, the values of K O and α are chosen to fit the curve proposed in [START_REF] Costache | Comprehensive model of microalgae photosynthesis rate as a function of culture conditions in photobioreactors[END_REF] as a sigmoid function. Finally, the values of Q min and Q max are deduced from Mairet et al. [START_REF] Mairet | Modelling neutral lipid production by the microalga Isochrysis aff. galbana under nitrogen limitation[END_REF]

m ML 20 kg m -3 d -1
Friction coefficient of M over L [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF][START_REF] Clarelli | A fluid dynamics multidimensional model of biofilm growth: stability, influence of environment and sensitivity[END_REF] m EL 20 kg m -3 d -1

Friction coefficient of E over L [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF][START_REF] Clarelli | A fluid dynamics multidimensional model of biofilm growth: stability, influence of environment and sensitivity[END_REF] m ME 20 kg m -3 d -1

Friction coefficient of M over E [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF][START_REF] Clarelli | A fluid dynamics multidimensional model of biofilm growth: stability, influence of environment and sensitivity[END_REF]

γ M 2.8 • 10 -7 kg m -1 d -2
Tensor coefficient for micro-algae M [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF][START_REF] Clarelli | A fluid dynamics multidimensional model of biofilm growth: stability, influence of environment and sensitivity[END_REF]]

γ E 2.8 • 10 -7 kg m -1 d -2
Tensor coefficient for extra-cellular matrix E [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF][START_REF] Clarelli | A fluid dynamics multidimensional model of biofilm growth: stability, influence of environment and sensitivity[END_REF] Table 1: Estimated values of the biological and physical parameters used in the system described in subsection 3.5

Numerical simulations

In this section, we present some numerical results in the one dimensional case. Although the one dimensional case is not realistic, it is very useful to make preliminary tests on the efficiency of the model and to analyse solutions of the system. For example, in our case, the following simulations help to understand effects of the limiting processes present in the expression of the reaction rates and the role of some key parameters. More precisely, after giving a first test case which is considered as a reference case, we study the influence of the following parameters : substrate supply θ S , intensity on the upper surface of the water I 0 , absorption coefficient τ M and variation of light. Note that other numerical simulations are presented in appendices to show the influence of the elastic tensor coefficients γ M and γ E and of the maximal photosynthesis rate µ P . The one dimensional setting has to be interpreted as the case when biofilms have a uniform horizontal distribution and the space variable under consideration is therefore the height. In all the simulations presented here, we take the following values: the length of the domain is equal to L x = 5 mm, the number of discretisation points is N x = 400 and the initial data are taken as: -4 ]. Boundary conditions are described in section 3.5; they are all homogeneous, except for the components dissolved in liquid on top of the domain (i.e. x = L x ), for which we take the same values θ S , θ C and θ O as for initial conditions above. Values of all the other parameters of the system can be found in Table 1. Finally, light intensity on the upper surface of water (x = L x ) is given by I 0 (t, L x ) = I opt and is independent of time.

A 0 = 5 • 10 -2 χ x 10 -4 , N 0 = 8.38 • 10 -3 χ x 10 -4 , E 0 = 0, L 0 = 1 -A 0 -N 0 -E 0 , S 0 = θ S , C 0 = θ C , O 0 = θ O ,

First tests

In Figure 1, we represent the volume fractions of the biofilm components, that is to say A, N , E and A + N + E, and the mass fractions of the components dissolved in liquid, that is to say S, C and O, for t = 10, 45 and 90 days. Note that the left boundary of the domain corresponds to the bottom of the tank and the right boundary to the top of the tank, where nutrients are brought to the domain. We can observe the evolution of a front, corresponding to the development of biofilm within water. The total volume fraction of biofilm increases with time at the location of the front (see the purple curve of A + N + E on Fig. 1a,1c,1e). Moreover we observe that the composition of the biofilm is not homogeneous in space: in the region of the front, the amount of micro-algae A (in green in Fig. 1) and of N (in red) is particularly high, whereas the region of the bottom is almost exclusively made of extra-cellular matrix (in brown).

Now, looking at the figures displaying the mass fractions of the components dissolved in liquid, we distinguish three regions. The interpretation of these regions is tightly linked to the limiting factors involved in the reaction rates and displayed in Fig. 4, that is to say, starting from the right side of the domain:

• On top of the tank, corresponding to x ∈ [2.4; 5] in Fig. 1f, the mass fractions of S and C are smaller than their values on the top boundary, which are fixed by boundary conditions, owing to their consumption by the biofilm at the front: x ∈ [2.2; 2.4], whereas the mass fraction of O is higher due to the significant release of oxygen by photosynthesis.

• In the area directly behind the front, corresponding to x ∈ [1.8; 2.2] in Fig. 1f, the mass fraction of C slowly increases due to release induced by respiration (see the red curve of f R in Fig. 4d) and due to absence of consumption because of the null photosynthesis rate (see the orange curve of f P in Fig. 4c). For oxygen, an opposite phenomenon occurs : since oxygen is no longer active, its mass fraction decreases due to respiration. In this area, substrate is slightly consumed and its mass fraction decreases slowly.

• Finally, on bottom of the tank, where biofilm is almost exclusively made of extracellular matrix, corresponding to x ∈ [0; 1.8] in Fig. 1f, mass fractions are nearly constant in space.

Assuming that biofilm is uniform in the two other directions, the daily production rate is estimated to 0.982 g m -2 d -1 of dry matter (considering that biofilm contains about 90% of water) at t = 10 days, 1.084 g m -2 d -1 at t = 45 days and 1.374 g m -2 d -1 at t = 90 days. On average, daily production rate is equal to 1.131 g m -2 d -1 . This value is comparable with biological experiments made by Schnurr et al. [START_REF] Schnurr | The effect of light direction and suspended cell concentrations on algal biofilm growth rates[END_REF][START_REF] Schnurr | Factors affecting algae biofilm growth and lipid production: A review[END_REF] or Gross et al. [START_REF] Gross | Development of a rotating algal biofilm growth system for attached microalgae growth with in situ biomass harvest[END_REF], in which daily productivity is estimated between 0.7 g m -2 d -1 and 6.5 g m -2 d -1 depending on conditions, species and cultivation methods.

Note that the mass fractions of components dissolved in liquid phase influence the reaction terms through the reaction rates. However, the amount of matter is given by the quantities SL, CL and OL that are represented in Figure 2. For convenience, we only describe Figure 2c but the two other figures in Figure 2 read similarly. On top of the tank (i.e. x ∈ [2.3; 5]) the volume fraction of liquid is equal to 1 so SL, CL and OL have the same behavior as their mass fractions S, C and O, that is to say the amount of S and of C decreases, whereas the mass fraction of O increases. In the biofilm front area (i.e. x ∈ [1.8; 2.3]), the volume fraction of liquid decreases suddenly, leading to a decrease of SL, CL and OL. On bottom of the tank, the amount of these quantities slowly increases along with the volume fraction of liquid.

Figure 3 represents the velocities of the different components at t=90 days. We observe that the velocities of micro-algae and of extra-cellular matrix behave similarly, that is to say they are positive around the biofilm front region, as expected from the evolution of biofilm front. Indeed, the front location, defined as the point where the gradient of A + N + E is the largest, moves with a speed equal to 26.1 µm d -1 , which is in the range of 23 -27 µm d -1 measured in [START_REF] Schnurr | The effect of light direction and suspended cell concentrations on algal biofilm growth rates[END_REF]. On the opposite, the velocity of liquid is negative, which can be explained by the fact that liquid is consumed by biofilm through photosynthesis. Figures 4a and4c represent functions f Q , f L , f C , f O and f I described at equations (3.7b), (3.7c), (3.7d), (3.7e) and involved in the expression of the photosynthesis reaction rate (3.9a), evaluated at t=45 and 90 days, respectively. At t = 45 days, we can observe that elementary functions f C and f L have values greater than 0.88, meaning that they are not limitant for photosynthesis. On the opposite, functions f I , f O and f Q take smaller values, meaning that photosynthesis is limited by the lack of light and functional biomass and inhibited by high concentration of oxygen. More precisely, 90% of light is absorbed in the first 240 µm after the biofilm front (i.e. x ∈ 0.82 • 10 -3 ; 1.06 • 10 -3 ) leading to f I (in yellow) equal to zero for x ∈ 0; 0.7 • 10 -3 . The function f Q (in red) is equal to 1 in the interval 0; 0.8 • 10 -3 and rapidly decreases to around 0.67 at x = 1.2 mm to stabilise around 0.57 in the last part of the domain. So, the quota of functional biomass is not limiting inside the biofilm but is partially limiting on the biofilm front. On the interval 0; 1.3 • 10 -3 , where oxygen concentration is high, see Fig. 1d, the function f O has values under 0.7, meaning that oxygen inhibits the photosynthesis process. Finally, the product

f Q f C f L f I f O f P (a)
f Q f C f L f I f O f P (c)
(Q) = max Q max -max{Q,Q min } Q max -Q min
f P = f Q f L f C f O f I ,
which is nearly equal to the photosynthesis rate is also plotted on the same figure (in orange) and it confirms that photosynthesis takes place on top of the biofilm, that is to say in the interval 1 • 10 -3 ; 1.2 • 10 -3 .

At t=90 days, we observe in Fig. 4c that f O is now above 0.82, meaning that oxygen concentration is hardly inhibiting. Indeed, according to Figures 1d and1f, the mass fraction of oxygen in biofilm area has decreased between 45 and 90 days. At the same time, the total volume fraction of biofilm has increased (see Fig. 1c and1e), leading the volume fraction of water to decrease. f L takes values around 0.82 in the area of the biofilm front, so the water becomes slightly limiting. For the light intensity and the effect of the quota of the functional biomass, we observe the same effects as at t=45 days: 90% of the light is absorbed in the first 250 µm after the biofilm front (see f I in yellow) and Q is not limiting inside the biofilm but is partially limiting in the biofilm front. Finally, we also plot f P , that confirms that photosynthesis still takes place on top of the biofilm (in the interval 2.2 • 10 -3 ; 2.4 • 10 -3 ), even if the order of influence of the elementary functions has changed. Now, in Figures 4b and4d, we show the other functions used to compute the reactions rates (3.9). First, we notice that Fig. 4d is basically a shift of Fig. 4b meaning that the influence of f R , f S , f N and f D , unlike the elementary functions that compose ϕ P , does not change inside the biofilm over time. We can also observe that the graph of f R , defined by f R (O) = O K O +O and used in respiration rate (3.9b), remains close to 0.8 in the biofilm area, so oxygen is almost not limitant for respiration. Then, concerning the functional biomass synthesis rate (3.9c), function f S (S) = S K S +S is nearly equal to 1, whereas

f N (Q) = max Q max -max{Q,Q min } Q max -Q min
vanishes in the biofilm region. Therefore, functional biomass synthesis happens mostly in the region of the biofilm front. Finally, f D (O) = 1 -

β Ô Ôβ +β-1 , with Ô = O O ref
, present in death rates (3.9g), takes higher values in the biofilm, since the presence of oxygen increases these rates.

Influence of the value of substrate supply θ S

We study in this subsection the influence of the value of substrate supply θ S .

Figure 5a represents the average daily production rate (estimated on 90 days) of dry biomass for the different constituents of the biofilm with respect to the parameter θ S . According to this graph, the daily production rate of dry biofilm (A + N + E) increases until θ S = 20 µgS/gL, then stabilises around 1.13 g m -2 d -1 . However, differences can be observed between the production rate of each component. Indeed the production rates of carbon pool (A) and functional biomass (N) increase linearly until θ S = 35 µgS/gL; then they stabilise at 1.63 g m -2 d -1 for A and 0.68 g m -2 d -1 for N. Behavior of ECM is different: its production rates increases until 9.77 g m -2 d -1 for θ S = 17.5 µgS/gL, then decreases and stabilises at 8.8 g m -2 d -1 for θ S > 35 µgS/gL. These production rate changes impact the biofilm composition, as it can be observed in Figure 5b. Indeed, from θ S = 0 to θ S = 35 µgS/gL, the percent of EPS composing the biofilm decreases from 93.1 to 77.2%, whereas the percentage of carbon pool (resp. functional biomass) increases from 6.1 to 16.7% (resp. 0.8 to 6.1%). Then, for θ S > 35 µgS/gL, the biofilm composition remains stable. The front velocity behaves as the total production rate, namely it increases until θ S = 15 µgS/gL and then stabilizes at 26.1 µm d -1 .

In order to better understand the influence of substrate supply, we present the results of numerical simulations performed with θ S = 1.65 • 10 of Section 4.1; the other parameters remain unchanged. Figure 6 represents volume and mass fractions at t = 90 days and in Fig. 7, we plot the elementary functions that compose the reaction rates.

We remark that these parameters lead to a substrate deficiency, since mass fraction S (in red in Fig. 6b) is equal to zero in the biofilm region. Substrate deficiency promotes extra-cellular matrix excretion: the volume fraction of the whole biofilm is almost the same as previously, with a difference of less than 5%, but the shape and the composition change significantly. The pool of carbon storage and the functional biomass volume fractions are 38.5% and 52.9% smaller, while the ECM volume fraction is 6.9% larger, see Figures 6a and1e. Moreover, the peak on the biofilm front is smaller. Regarding the elementary functions that compose the reaction rates, we notice that the photo-synthesis rate is the same as in previous simulations, see Fig. 7a. However, the functional biomass synthesis rate (3.9c) is influenced by function f S (S) = S K S +S , which is nearly equal to zero in Fig. 7a inside the biofilm region x ∈ [0; 2.1] mm. This effect is compensated by the functional biomass quota Q through the function f

N (Q) = max Q lim -max{Q,Q min } Q lim -Q min
, which is non zero above 2 mm (instead of 2.1 mm in Fig. 4d).

f Q f C f L f I f O f P (a)
f R (O) = O K O +O present in respiration rate (3.9b), f S (S) = S K S +S and f N (Q) = max Q lim -max{Q,Q min } Q lim -Q min

Influence of light intensity on the upper surface of water: I 0 (t, L x )

In this subsection, we investigate the influence of light intensity on the average daily production rate and on the velocity of the biofilm front.

Figure 8a represents the average daily production rate (estimated on 90 days) of dry biomass for the different constituents of the biofilm with respect to the light intensity on the upper surface I 0 (t, L x ). We can observe that the average production rate for A, N and E increases until I 0 (t, L x ) = 1.8I opt and then decreases. So the maximal productivity is reached for I 0 (t, L x ) = 1.8I opt , whereas I opt is the light intensity for which f I is maximal in the photosynthesis rate. However, the improvement of daily productivity for I 0 (t, L x ) = 1.8I opt compared to the productivity for I 0 (t, L x ) = I opt is very low: ∼ 3%. We can also observe that the average production rate of A and N for I 0 (t, L x ) > I opt is constant, meaning that the light intensity on the upper surface of water I 0 (t, L x ) does not change the productivity of these compounds. Now, Figure 8b shows that the velocity of the biofilm front is maximal for I 0 (t, L x ) = 0.6I opt , that is to say the maximal daily production rate and the maximal velocity of the front are not reached for the same values of I 0 (t, L x ).

As we can remark in Figure 9, for a light intensity on the upper surface I 0 (t, L x ) = 1.8I opt , shape of the biofilm components (A, N, E) and mass fraction of the substrate (S) are comparable to the case when I 0 (t, L x ) = I opt , represented in Fig. 1e and1f. But in the biofilm region, for a higher value of light intensity, mass fraction of inorganic carbon is smaller whereas mass fraction of oxygen is greater.

f Q f C f L f I f O f P (a)
f R f S f N f D (b)
Other elementary functions that compose the other reaction rates (3.9) at t=90 days.

Figure 10: Plots at t = 90 days of the elementary functions used to compute the reactions rates with respect to space. On graph 10a (on the left), we display the functions (3.7) used for the photosynthesis rate and their product f P ; on graph 10b (on the right), we plot

f R (O) = O K O +O present in the respiration rate (3.9b), f S (S) = S K S +S and f N (Q) = max Q max -max{Q,Q min } Q max -Q min present in the functional biomass synthesis rate (3.9c) and f D (O) = 1 - β Ô Ôβ +β-1 , with Ô = O O ref
, present in the death rates (3.9g) and (3.9h).

The curves that represent in Fig. 10b the elementary functions used in the reaction rates differ also partly, especially f I . Outside of the biofilm region, namely x ∈ 2.2 • 10 -3 ; 5 • 10 -3 , the value of f I is about 0.86, meaning that if micro-algae were present in this area, photosynthesis would be inhibited due to a too high light intensity. Then, gradually, as light penetrates into the biofilm and is absorbed, its intensity decreases until it reaches the optimal light intensity, here for x = 2.1 mm. Finally, as light intensity continues to be absorbed and thus decreases, f I decreases until it reaches 0 for x ∼ 1.8 mm. We can also observe that f O takes smaller values in the area of the biofilm front, when the light intensity on the upper surface is higher.

Influence of light absorption coefficient for the biofilm τ M

In this subsection, we investigate the influence of the light absorption coefficient of micro-algae τ M , see equation (3.8), which accounts directly for the light penetration in the biofilm. In Figure 11 we represent the volume and mass fractions for the constituents considered in the model with the light absorption parameter of micro-algae set to τ M = 1.25 • 10 4 m -1 , i.e. twice smaller than the simulations of Fig. 4.1. First of all, although the biofilm shape (in purple) remains consistent with the results presented in the first test case, see Fig. 1e, we can observe that the biofilm composition differs. Indeed after 90 days, the volume fraction of A (in green) is 21.9% larger and the volume fraction of ECM (in brown) is 2.7% larger, whereas the volume fraction of functional biomass (in red) is 5.7% smaller. These discrepancies lead to an increase of 5.4% for the total mass of the biofilm. In Figure 11a, we can also notice that the biofilm front has only reached x = 2 mm whereas in the first test case in Fig. 1e, it has reached x = 2.4 mm, meaning that the velocity of the biofilm front is smaller. Indeed, the front velocity t is estimated to 21.4 µm d -1 , instead of 26.1 µm d -1 in Fig. 4.1, so the light absorption parameter has an indirect impact on the front velocity.

As regards the components dissolved in liquid represented in Fig. 11b, we can observe that, in the biofilm area, mass fractions for carbon dioxyde and for substrate are lower while oxygen mass fraction is higher compared to Fig. 1f. In order to explain these differences, let us take a look at the elementary functions used to build the reaction rates and represented in Figure 12. Although light intensity through function f I remains the main limitant parameter, the decrease as it penetrates into the biofilm is slower. Moreover we can observe that f O (in cyan) takes values below 0.6 in the area of the biofilm front, so photosynthesis is here inhibited by the excess of oxygen. We can also notice that the function f P (in orange), which is nearly the photosynthesis rate, decreases more smoothly, so photosynthesis occurs deeper into the biofilm. Regarding the other elementary functions used to build the reaction rates and represented in Fig. 12b, only f N has a different shape: like photosynthesis, assimilation of substrate continues to occur deeper into the biofilm.

Influence of variation of light

Finally, we consider that light intensity on top of the domain depends on time as follows:

I 0 (t, L x ) = 3I opt max {0, sin (2πt)} ,
such that the biofilm is periodically enlightened with an intensity between 0 and 3I opt . Results can be observed in Figure 13 for the volume fractions and in Figures 14 and15 for the mass fractions and the elementary functions that compose the reaction rates. Average daily production rate over the first 150 days is equal to 0.495 g m -2 d -1 and the biofilm front moves with a speed equal to 15.2 µm d -1 . Volume fractions and composition of the biofilm are unchanged. However the mass fractions of the components dissolved into liquid phase (substrate, inorganic carbon and oxygen) represented in Fig. 14 and 15 evolve all day long with the enlightenment of the biofilm. Indeed, assimilation of inorganic carbon and release of oxygen are tightly linked to photosynthesis and so to light intensity. Regarding the elementary functions that compose the photosynthesis rate and displayed in the middle column in Fig. 14 and15, we can observe that the limiting factors are light intensity f I and oxygen excess f O defined in equations (3.7e) and (3.7d) respectively. As long as light intensity increases and remains smaller than I opt , the photosynthesis rate on the biofilm front increases; then, when light intensity becomes greater than I opt , it starts to inhibit photosynthesis on the front but enables micro-algae to make photosynthesis behind the front (i.

f Q f C f L f I f O f P (a)
f R f S f N f D (b)
Other elementary functions that compose the other reaction rates (3.9) at t=90 days.

Figure 12: Plots at t = 90 days of the elementary functions used to compute the reactions rates with respect to space, in the case a smaller light absorption coefficient for the biofilm, i.e. τ M = 1.25 • 10 4 m -1 . On the left, we display the functions (3.7) used for the photosynthesis rate and their product f P ; on the right, we plot

f R (O) = O K O +O present in the respiration rate (3.9b), f S (S) = S K S +S and f N (Q) = max Q lim -max{Q,Q min } Q lim -Q min
present in the functional biomass synthesis rate (3.9c) and

f D (O) = 1 - β Ô Ôβ +β-1 , with Ô = O O ref
, present in the death rates (3.9g) and (3.9h).

to an inhibition of photosynthesis; this phenomenon is confirmed by the curve of f O in Fig. 15b and 15f that takes small values in the area of the biofilm front. (3.9h). In this simulation, substrate supply is in excess which is confirmed by the values of f S that remain close to 1. During night (first row in Fig. 14 and15), in the biofilm area, f D takes values around 0.2 meaning there is death induced by lack of oxygen. The deficiency of oxygen is caused by the respiration process which is represented by f R , which is equal to about 0.2 in the biofilm area. Moreover, oxygen concentration increases with the light intensity which induces f R to increase: in the biofilm area, f R goes from 0.2 for I 0 (t, L x ) = 0 to 0.8 when I 0 (t, L x ) = 3I opt . Finally, let us take a look at the death process (see the curve of f D ) which increases with the oxygen concentration, so it is minimal for I 0 (t, L x ) = I opt /2 and it increases up to 0.2 for I 0 (t, L x ) = 3I opt .

Conclusion

In this paper, we have proposed a new mixture model for micro-algae biofilm producing lipids. The main originality of this model is the precise description of the biological mechanisms involved in its formation and development, combined with the use of partial differential equations to take into account the spatial variations of the biofilm. In particular, we use Droop's theory to describe the assimilation of carbon through photosynthesis and we model finely functional biomass growth, extra-cellular matrix excretion process and variation of light intensity. The general form of our model is the same as the model in [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF], but it differs on several points, among which the source terms of the mass balance equations and therefore the incompressibility constraint, the consideration of a third velocity in force balance equations and the addition of a diffusion term for the components dissolved in liquid.

In addition, we have proposed a well adapted numerical scheme, able to deal with the stiffness of the reaction terms. This numerical scheme can be straightforwardly extended in the two and three dimensional cases. However, the model contains a great number of parameters: the values of many of them are given in the literature, but some of them remain unknown and have to be calibrated. Consequently, in a further work, a sensibility analysis, as well as biological experiments, are scheduled in order to set up more precisely the parameters. A careful confrontation of the model in 2D or 3D with experimental results is foreseen. Some first tests are performed in the one-dimensional case. This case is not realistic, but it remains relevant for a uniform horizontal biofilm, the spatial variable being the height. Moreover, it enables to obtain some first results on the model. Comparison of the numerical simulations presented here with experimental results that can be found in the literature shows that our model predicts realistic daily production rates and front speeds, which gives a first validation. In these simulations, we can observe that the biofilm is composed of three distinct regions: the front, made of micro-algae performing photosynthesis, the middle of the biofilm composed of extra-cellular matrix and micro-algae and finally the bottom of the tank where only extra-cellular matrix remains. We expect this composition of the biofilm to hold true also in the two-dimensional and the threedimensional cases and to observe in these cases some mushroom-shaped structure, as described in [START_REF] Zhang | Phase field models for biofilms. i. theory and onedimensional simulations[END_REF][START_REF] Zhang | Phase-field models for biofilms ii. 2-D numerical simulations of biofilm-flow interaction[END_REF].

This first model can be extended to more complicated domains than a simple tank: several biological teams experiment the idea of cultivating micro-algae biofilms on rotative systems, where the biofilm is periodically immersed and enlightened. It would be of great interest to adapt the present model to this case and to study the dependence of such a model on factors, such as light frequency, light intensity or rotating frequency.

A Numerical scheme

The aim of this section is to present the numerical scheme we use to discretize the full set of equations described in Subsection 3.5. This scheme follows, for the spatial discretization, the strategy of [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF] derived for a similar model, that is to say some finite-differences scheme derived from relaxation techniques. In this article, the authors dealt with the two following issues: the computation of the velocities in the case of a vanishing phase and the computation of the pressure term. However, we have to face another difficulty here: the mass exchanges between components x (in mm) mass fractions are non linear and stiff, so we use some explicit-implicit treatment of the source terms of the mass balance equations in order to preserve the non-negativity of the solutions.

f Q f C f L f I f O f P (b)
f Q f C f L f I f O f P (b)
We consider the one dimensional domain Ω = [0, L x ] of width L x . The extension of the following scheme to the 2D and 3D cases is straightforward, but in this paper we focus our attention on the one dimensional case. We denote by ∆x the space step in the x direction and we use a regular grid on Ω, namely the discretization points x i = i∆x, 1 ≤ i ≤ N x . The k th time step is denoted by ∆t k and the time discretizations are therefore equal to t n = n k=1 ∆t k . We begin with rewriting the full system of Subsection 3.5 under the following form:

∂ t U + ∂ x F (U , W ) = Γ(U ) + ∂ x M δ L∂ x U L , (A.1a) ∂ t W + ∂ x G (U , W ) = G I (U , W ) + G T (U , W ) + G P (U ), (A.1b) L = 1 -A -N -E, (A.1c) ∂ x ((A + N )v M + Ev E + Lv L ) = Γ A + Γ N ρ M + Γ E ρ E + Γ L ρ L , (A.1d)
that is to say one equation for the mass balances, one equation for the force balances, the volume condition and the incompressibility constraint. Here, U is a vector containing the mass fractions, W is vector containing the three velocity, F (resp. G) is the flux of vector U (resp. of vector W ), M δ is the diffusion matrix for U , Γ(U ) is the mass exchanges source term in the mass balance equations. The source terms in the force balance equations are split into three parts: the pressure part G P (U ), the interactions term G I (U , W ) and the mass exchanges term G T (U , W ), that is to say:

U =                       A N E SL CL OL                       , F (U , W ) =                       Av M N v M Ev E SLv L CLv L OLv L                       , Γ (U ) =                       Γ A • ρ -1 M Γ N • ρ -1 M Γ E • ρ -1 E Γ S • ρ -1 L Γ C • ρ -1 L Γ O • ρ -1 L                       , M δ =                           0 0 0 0 δ S δ C δ O                           , W =          (A + N )v M Ev E Lv L          , G (U , W ) =            (A + N ) v M 2 + γ M E v E 2 + γ E Lv L 2            , G P (U ) =          -(A + N ) ∂ x P -E∂ x P -L∂ x P          , G I (U , W ) =           -m ML ρ M (v M -v L ) -m ME ρ M (v M -v E ) -m EL ρ E (v E -v L ) + m ME ρ E (v M -v E ) m ML ρ L (v M -v L ) + m EL ρ L (v E -v L )           , G T (U , W ) =             Γ A +Γ N ρ M v M Γ E ρ E v E -1 ρ L ((Γ A + Γ N )v M + Γ E v E )             .
We denote by U n i and W n i the discrete approximations of U (t n , x i ) and W (t n , x i ).

A.1 Numerical approximation of mass balance equations

The first step consists in solving equation (A.1a), using an explicit discretization based on a relaxation technique, presented in [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF] and used in [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF], for the flux term and a mixed strategy for the diffusion term and the source term. More precisely for the diffusion, U is treated implicitly to reduce the CFL condition, whereas L is explicitly. The spatial discretization is done with classical finite differences method which involves numerical approximations of L on the edges of the mesh cells; they are interpolated as the mean values of L on the cells sharing this edge, namely we use the approximation L i+ 1 2 = 1 2 (L i+1 + L i ). The strategy for the source term is explained further on. We obtain therefore a scheme of the form:

U n+1 i -∆tΓ U n i , U n+1 i - ∆t ∆x 2       L n i+ 1 2       U n+1 i+1 L n i+1 - U n+1 i L n i       -L n i-1 2       U n+1 i L n i - U n+1 i-1 L n i-1             = U n i - ∆t 2∆x F W n i+1 -F W n i-1 + λ ∆t 4∆x U n i+1 -2U n i + U n i-1 (A.2)
where the numerical velocity λ is common to the mass balance equations and to the force balance equations and is equal to the maximum of the eigenvalues of the Jacobian matrix of the fluxes F G , that is to say:

λ = max |2v L | , |v M | + γ M ρ M , |v E | + γ E ρ E .
For stability reasons, the time step ∆t is computed at each time step such that the stability condition λ ∆t ∆x ≤ 1 is satisfied. Now, let us explain how to compute the source terms with the explicit -implicit strategy, in order to guarantee that each component of U remains in [0, 1]. From now on in this subsection, we will drop the spatial indices, since no confusion is possible. We therefore discretise, at each point of the domain, the photosynthesis rate, the respiration rate and the functional biomass synthesis rate using implicit expressions for SL, CL and OL at the numerator and explicit expressions for all the other terms, which leads to the following discretizations:

ϕ P (U n , U n+1 ) = µ A ρ M N n 1 - Q min Q n (CL) n+1 (K L + L n ) (K C + C n ) , ϕ R (U n , U n+1 ) = µ R ρ M A n (OL) n+1 L n K R + (OL) n , ϕ N (U n , U n+1 ) = µ N ρ M N n (SL) n+1 L n K S + (SL) n Q max -Q n Q max -Q min .
Note that the two other reaction rates, namely the EPS excretion rate and the micro-algae death rate, are treated explicitly. The source terms are consequently approximated at each point of the domain by:

Γ A U n , U n+1 = 1 ρ M ϕ P (U n , U n+1 ) -ϕ R (U n , U n+1 ) -η A N ϕ N (U n , U n+1 ) -ϕ A E (U n ) -ϕ A D (U n ) , Γ N U n , U n+1 = 1 ρ M ϕ N (U n , U n+1 ) -ϕ N E (U n ) -ϕ N D (U n ) , Γ E (U n ) = 1 ρ E ϕ A E (U n ) + ϕ N E (U n ) + ϕ A D (U n ) + ϕ N D (U n ) , Γ S U n , U n+1 = - η S N ρ L ϕ N (U n , U n+1 ), Γ C U n , U n+1 = 1 ρ L -η C P ϕ P (U n , U n+1 ) + η C R ϕ R (U n , U n+1 ) , Γ O U n , U n+1 = 1 ρ L η O P ϕ P (U n , U n+1 ) -η O R ϕ R (U n , U n+1 ) .
In practice, since we do not treat implicitly any variable in the denominator of the source terms, each component of U n+1 can be easily computed through the resolution of linear systems. The overall strategy at time step t n+1 is the following: we first obtain the volume fractions (SL) n+1 , (CL) n+1 and (OL) n+1 as the solutions of linear systems, which allows to compute the reaction rates ϕ P (U n , U n+1 ), ϕ R (U n , U n+1 ) and ϕ N (U n , U n+1 ). It is then straightforward to compute the remaining components of U n+1 , that is to say A n+1 , N n+1 and E n+1 . Finally, an approximation of the volume fraction of liquid is given thanks to condition (A.1c).

A.2 Numerical approximation of force balance equations

Dealing with the computation of force balance equations, we face two difficulties: vanishing phases and computation of the pressure term.

Firstly, force balance equations give the evolution of the momentum of each component, whereas the friction forces depend on the velocities. As a consequence, when one of the phases is vanishing, it is not clear how to define its velocity, that is needed to approximate the friction forces at the following time step. Note that, in a biological context, situations where for example L = 1 and A = N = E = 0 are relevant and we cannot claim to be far from vacuum, as can be done in a physical context. Therefore, in order to compute the velocities, we use an implicit-explicit time discretization strategy for the momentum equations, where the interaction forces term G I (U , W ) is treated implicitly, see [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF].

Secondly, to compute the velocities, we need to know the gradient of the hydrostatic pressure, which is another unknown of the system. A natural approach consists in finding an equation verified by P , by taking the divergence of the sum of the momentum equations and by using eq. (A.1d). However, this method, which is known to be inefficient, leads to an elliptic equation for P with a non-unique solution. To overcome this difficulty we use a splitting approach which is basically an adaptation of the Chorin-Temam projection method [START_REF] Temam | Une méthode d'approximation de la solution des équations de Navier-Stokes[END_REF][START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF], see again [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF]. This method uses a projection-correction approach: first, we compute an approximation of the velocities using the force balance equations without the pressure terms; then, using the predicted velocities and the average incompressibility constraint, we compute the pressure as the solution of an elliptic equation. Finally, the velocities are corrected thanks to the value of ∂ x P .

Let us recall that the vector U n+1 is computed at a previous step. We now give some details on the scheme for the force balance equations computation.

First, the projection step enables us to estimate some predicted values V n+ 1 2 for the velocities and this first step of the scheme writes as:

M n+1 i V n+ 1 2 i = W n i - ∆t 2∆x (G (U , W ) n i+1 -G (U , W ) n i-1 ) + λ ∆t 4∆x W n i+1 -2W n i + W n i-1 + ∆tG T U n+1 i , W n i (A.3) where V =          v M v E v L          , G T U n+1 , W n =                       Γ M n+1 ρ M v M n Γ E n+1 ρ E v E n - 1 ρ L Γ M n+1 v M n + Γ E n+1 v E n                       , M n+1 i =                   
) A + N ρ E ρ L + E ρ M ρ L + L ρ M ρ E ,
which is a polynomial of degree 2 in ∆t with positive coefficients. Since all the friction coefficients are strictly positive and since the constraint A + N + E + L = 1, with A, N , E, L positive, is satisfied, the leading coefficient of det M is non zero. So, for all ∆t strictly positive, det M does not vanish and linear system (A.3) has a unique solution.

Then in the second step of splitting we have to solve

∂ t φv φ = - φ ρ φ ∂ x P (A.4)
for φ = M, E, L in the interval [t, t + ∆t n ] with initial data φ n+1 v n+ 1 2 φ . The discrete approximation of these equations for φ = M, E, L is given by

φ n+1 v n+1 φ -φ n+1 v n+ 1 2 φ = -∆t φ n+1 ρ φ ∂ x P n+1 .
Then taking the divergence of the sum of these equations over φ we get

∂ x A n+1 + N n+1 v n+1 M + E n+1 v n+1 E + L n+1 v n+1 L -∂ x A n+1 + N n+1 v n+ 1 2 M + E n+1 v n+ 1 2 E + L n+1 v n+ 1 2 L = -∆t∂ x A n+1 + N n+1 ρ M + E n+1 ρ E + L n+1 ρ L ∂ x P n+1 .
Now using the discrete approximation of the incompressibility constraint (A.1d), we deduce that P is solution of an elliptic equation with non constant coefficients:

∆t∂ x                 A n+1 + N n+1 ρ M + E n+1 ρ E + L n+1 ρ L         ∂ x P n+1         (A.5) = ∂ x A n+1 + N n+1 v n+ 1 2 M + E n+1 v n+ 1 2 E + L n+1 v n+ 1 2 L - Γ n+1 A + Γ n+1 N ρ M - Γ n+1 E ρ E - Γ n+1 L ρ L (A.6)
This equation is completed with Neumann boundary conditions on P . This boundary condition can be directly deduced from boundary conditions on the velocities and equation (A.4) on the boundary defined by x = 0. For x = L x , we need to note that, in this step, v T does not change since none of the volume or mass fractions are modified, so using equation (A.4) we also get Neumann boundary condition for the pressure on x = L x . However, with Neumann boundary conditions on P , equation (A.5) does not have a unique solution. In practice, we choose the solution for which the average value of P is null, that is to say the solution satisfying Finally in order to get a symmetric system we rather solve the minimisation problem

∆t∂ x A n+1 + N n+1 ρ M + E n+1 ρ E + L n+1 ρ L ∂ x P n+1 + λ = ∂ x A n+1 + N n+1 v n+ 1 2 M + E n+1 v n+ 1 2 E + L n+1 v n+ 1 2 L - Γ n+1 A + Γ n+1 N ρ M - Γ n+1 E ρ E - Γ n+1
L ρ L , L x 0 P n+1 dy = 0, ∂ x P n+1 (0) = ∂ x P n+1 (L x ) = 0, with λ the Lagrange multiplier associated to the constraint L x 0 P n+1 dy = 0. This system is discretized using classical centered finite differences method and leads to the resolution of a linear symmetric system.

Finally since the volume fractions do not change in this step, we update velocities as follows:

v n+1 M = v n+ 1 2 M - ∆t ρ M (∂ x P ) n+1 , v n+1 E = v n+ 1 2 E - ∆t ρ E (∂ x P ) n+1 , v n+1 L = v n+ 1 2 L - ∆t ρ L (∂ x P ) n+1 .
B Numerical simulations B.1 Influence of the values of elastic interaction tensor coefficients γ M and γ E

According to [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF], the value of the elastic tensor coefficient γ drives the speed of the biofilm front. In our model, unlike in [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF], we consider two distinct elastic tensor coefficients: γ M for micro-algae and γ E for extra-cellular matrix. In . The values of the volume fractions and the composition of the biofilm are the same as in Fig. 4.1, with a difference of less than 3.2%. However, the velocity of extra-cellular matrix is lower than the micro-algae velocity in this case. Note that when γ E ≥ γ M , the conclusion is similar: the composition of the biofilm does not change and the velocity of extra-cellular matrix is greater than the micro-algae velocity. In practice, the influence of γ E appears to be low, unlike γ M that drives the velocity of the biofilm front.

As in [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF], we can observe in Figure 18a that the front velocity depends linearly on √ γ M . Indeed, the front velocity determined with numerical experiments (red marks) fits the linear regression (in blue) given by y = 4.4311 • 10 -2 x + 2.6488 • 10 -6 with correlation coefficient equal to 0.99996. In Fig. 18b we represent the average daily production rate (defined as the difference of dry biomass between two consecutive days) of dry biomass over 90 days. Again we observe that numerical experiments (red marks) fit the linear regression (in blue) given by y = 517.04x + 0.85866 with correlation coefficient equal to 0.9997.

B.2 Influence of maximal photosynthesis rate µ P

In this subsection, we investigate the influence of maximal photosynthesis rate µ P . Figure 19a represents the average daily production rate (estimated on 90 days) of dry biomass for the different constituents of the biofilm with respect to maximal photosynthesis rate µ P . As it can be expected, the production rate of pool of carbon storage increases with the value of the maximal photosynthesis rate. Since all other biofilm components (N and E) are synthesised from carbohydrate of the pool of carbon storage, their production rates also increase with the maximal photosynthesis. So the growth of the whole biofilm (A, N and E) increases with µ P . Consequently, the velocity of the front, which is represented on figure 19b, also rises with the maximal photosynthesis rate. However, the production rates and the velocity of the biofilm front increase more slowly when µ P increases, meaning that biofilm growth is limited. Indeed, as µ P increases, the photosynthesis rate rises and also leads to an increase of oxygen released by this mechanism. So, the oxygen mass fraction near the biofilm front becomes very high and inhibits increasingly the photosynthesis rate.

2. 2 . 1

 21 Micro-algae growth Now, let us describe what are the main mechanisms that occur during the formation of micro-algae biofilms.

with θ S = 5 •

 5 10 -5 kgS/kgL θ C = 100 • 10 -5 kgC/kgL θ O = 7.2 • 10 -6 kgO/kgL, and where χ x 10 -4 denotes the characteristic function of interval [0, 10

Figure 1 :Figure 2 :

 12 Figure 1: On Figures 1a, 1c and 1e (on the left), volume fractions of A (in green), N (in red), E (in brown) and A + N + E (in purple) with respect to space at t = 10 days (on top), t = 45 days (in the middle) and t = 90 days (on bottom) respectively; on Figures 1b, 1d and 1f (on the right), mass fractions of S (in red) , C (in green) and O (in blue) with respect to space at t = 10 days, t = 45 days and t = 90 days.

Figure 3 :

 3 Figure 3: Velocities of micro-algae (on the left), extra-cellular matrix (in the middle) and liquid (on the right) with respect to space at t = 90 days.

  Elementary functions (3.7) that compose the photosynthesis rate (3.9a) at t=45 days. Other elementary functions that compose the other reaction rates (3.9) at t=45 days.

  Elementary functions (3.7) that compose the photosynthesis rate (3.9a) at t=90 days. Other elementary functions that compose the other reaction rates (3.9) at t=90 days.

Figure 4 :

 4 Figure 4: Plots at t = 45 days (on top) and t = 90 days (on bottom) of the elementary functions used to compute the reactions rates, with respect to space. On the left, we display the functions (3.7) used for the photosynthesis rate and their product, denoted by f P ; on the right, we plot f R (O) = O K O +O

  synthesis rate (3.9c) and f D (O) = 1 -

Figure 5 :

 5 Figure5: Average production rate of dry constituent in g m -2 d -1 for 90 days with respect to substrate supply θ S (on the left) and biofilm composition in percent with respect to substrate supply θ S (on the right).

Figure 6 :

 6 Figure 6: On graph 6a (on the left), volume fractions of A (in green), N (in red), E (in brown) and A+N +E (in purple) and, on graph 6b (on the right), mass fractions of S (in red) , C (in green) and O (in blue) with respect to space at t = 90 days in the case of substrate deficiency, i.e. θ S = 1.65 • 10 -5 .

  Elementary functions (3.7) that compose the photosynthesis rate (3.9a) at t=90 days. Other elementary functions that compose the other reaction rates (3.9) at t=90 days.

Figure 7 :

 7 Figure 7: Plots at t = 90 days of elementary functions used to compute the reactions rates with respect to space, in the case of substrate deficiency θ S = 1.65 • 10 -5 . On graph 7a (on the left), we display the functions (3.7) used for photosynthesis rate and their product f P ; on graph 7b (on the right), we plot f R (O) = O K O +O present in respiration rate (3.9b), f S (S) = S K S +S and f N (Q) = max Q lim -max{Q,Q min }

  present in functional biomass synthesis rate (3.9c) and f D (O) = 1death rates (3.9g) and (3.9h).

5 Front velocity m d - 1 (Figure 8 :

 518 Figure8: Average production rate of dry constituent in g m -2 d -1 for 90 days with respect to light intensity on the upper surface of water in graph 8a (on the left) and average velocity of the biofilm front with respect to light intensity in graph 8b (on the right).

Figure 9 :

 9 Figure9: On graph 9a (on the left), volume fractions of A (in green), N (in red), E (in brown) and A + N + E (in purple) and, on graph 9b (on the right), mass fractions of S (in red) , C (in green) and O (in blue) with respect to space at t = 90 days, when I 0 (t, L x ) = 1.8I opt .

  Elementary functions (3.7) that compose the photosynthesis rate (3.9a) at t=90 days.

Figure 11 :

 11 Figure 11: On the left, volume fractions of A (in green), N (in red), E (in brown) and A + N + E (in purple) and, on the right, mass fractions of S (in red) , C (in green) and O (in blue) with respect to space at t = 60 days in the case of a smaller light absorption coefficient for the biofilm, i.e. τ M = 1.25 • 10 4 m -1 .

  e.: x ∈ [2.6; 3]). At the same time, oxygen saturation increases due to the released induced photosynthesis and leads 29

  Elementary functions (3.7) that compose the photosynthesis rate (3.9a) at t=90 days.

Figure 13 :

 13 Figure13: Volume fractions of A (in green), N (in red), E (in brown) and A + N + E (in purple) with respect to space at t = 150 days, when the biofilm is periodically enlightened.

2 0 2 0 2 0

 222 Elementary functions (3.7) that compose the photosynthesis rate (3.9a) for I 0 (t, L x ) Disolved components for I 0 (t, L x ) = I opt /photosynthesis rate (3.9a) for I 0 (t, L x ) = I opt /Other elementary functions that compose the other reaction rates (3.9) for I 0 (t, L x ) = I opt /

  photosynthesis rate (3.9a) for I 0 (t, L x ) = I opt Other elementary functions that compose the other reaction rates (3.9) for I 0 (t, L x ) = I opt

Figure 14 :

 14 Figure14: Plots of the mass fractions on the components dissolved in liquid (on the left) with respect to space for different times of the 150 th day: night (on top), I 0 (150.027, L x ) = I opt /2 (in the middle), I 0 (150.0545, L x ) = I opt (on bottom). In the middle and on the right, the elementary functions used to compute the reactions rates with respect to space; more precisely: in the middle, the elementary functions (3.7) used for the photosynthesis rate and their product f P and, on the right, we plotf R (O) = O K O +O present in the respiration rate (3.9b), f S (S) = S K S +S and f N (Q) = max Q max -max{Q,Q min } Q max -Q min

Figure 15 :

 15 Figure 15: Plots of the mass fractions of the components dissolved in liquid (on the left) with respect to space for different times of the 150 th day: I 0 (150.1165, L x ) = 2I opt (on top) and I 0 (150.25, L x ) = 3I opt (on bottom). In the middle and on the right, the elementary functions used to compute the reactions rates with respect to space; more precisely: in the middle the elementary functions (3.7) used for the photosynthesis rate and their product f P and on the right we plot f R (O) = O K O +O

L x 0 P

 0 (t, y) dy = 0.

Figure 4 . 1 ,

 41 simulations were performed with γ M = γ E ; however, in Figures 16 and 17, we show the volume and mass fractions and the velocities obtained for γ M = 2.8 • 10 -7 kg m -1 d -2 and γ E = 5 • 10 -8 kg m -1 d -2

Figure 16 :

 16 Figure 16: On the left, volume fractions of A (in green), N (in red), E (in brown) and A + N + E (in purple) and, on the right, mass fractions of S (in red) , C (in green) and O (in blue) with respect to space at t = 90 days, when γ E = 3 • 10 -10 m 2 d -2 and γ M = 3 • 10 -9 m 2 d -2 .

Figure 17 :Figure 18 :

 1718 Figure 17: Velocities of micro-algae, extra-cellular matrix and liquid at t=90 when γ E = 3 • 10 -10 m 2 d -2 and γ M = 3 • 10 -9 m 2 d -2 .

  Velocity of the front with respect to µ P

Figure 19 :

 19 Figure19: Average production rate of dry constituent in g m -2 d -1 for 90 days with respect to the maximal photosynthesis rate in graph 19a (on the left) and average velocity of the biofilm front with respect to the maximal photosynthesis rate in graph 19b (on the right).

  2.3.1, with reaction rate ϕ P . Moreover, according to Sec. 2.3.2, 2.3.3, 2.3.4 and 2.3.5, it is a reactant for all the other mechanisms: A is consumed by respiration with rate ϕ R (see Sec. 2.2), by synthesis of functional biomass (2.

	3) E and by micro-algae death (2.5a) with N ϕ N , by the EPS excretion (2.4a) with rate ϕ A with rate η A rate ϕ A D .

  Note that the matrix M contains also the coefficients of the interaction terms, which are treated implicitly, and that this matrix can be computed thanks to the knowledge of U n+1 . Moreover, we can compute the determinant of M, namely det M = (A + N )EL EL m ML + m ME m ML + m ME m EL

	  (m + ∆t A n+1 i + N n+1 i + ∆t ρ M (m ML + m ME ) -∆t ρ M m ME -∆t ρ M m ML -∆t ρ E m ME E n+1 i + ∆t ρ E (m EL + m ME ) -∆t ρ E m EL -∆t ρ L m ML -∆t ρ L m EL L n+1 i + ∆t ρ L m ML ρ L (A + N )E + m EL ρ L (A + N )E + m EL ρ E (A + N )L + m ME ρ E (A + N )L + m ML ρ M EL + + ∆t 2 (m	m ME ρ M	EL

ML + m EL )                      .