# Numerical schemes for mixture theory models with filling constraint: application to biofilm ecosystems

Olivier Bernard, Mickael Bestard, Thierry Goudon, Leo Meyer, Sebastian

Minjeaud, Florent Noisette, Bastien Polizzi

► To cite this version:

Olivier Bernard, Mickael Bestard, Thierry Goudon, Leo Meyer, Sebastian Minjeaud, et al.. Numerical schemes for mixture theory models with filling constraint: application to biofilm ecosystems. ESAIM: Proceedings and Surveys, 2024, CEMRACS 2022 - Transport in physics, biology and urban traffic, 77, pp.249 - 266. 10.1051/proc/202477249 . hal-04907848

HAL Id: hal-04907848

https://hal.science/hal-04907848v1

Submitted on 23 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
NUMERICAL SCHEMES FOR MIXTURE THEORY MODELS WITH FILLING
CONSTRAINT: APPLICATION TO BIOFILM ECOSYSTEMS

OLIVIER BERNARD<fnref n="1" />, MICKAEL BESTARD<fnref n="2" />, THIERRY GOUDON<fnref n="3" />, LEO MEYER<fnref n="4" />, SEBASTIAN MINJEAUD<fnref n="5" />, FLORENT NOISETTE<fnref n="6" /> AND BASTIEN POLIZZI<fnref n="7" />

**Abstract.** The mixture theory framework is a powerful way to describe multi-phasic systems at an intermediary scale between microscopic and macroscopic scales. In particular, mixture theory reveals a powerful approach to represent microbial biofilms where a consortium of cells is embedded in a polymeric structure. To simulate a model of microalgal biofilm, we propose an upgraded numerical scheme, consolidating the one proposed by Berthelin et al. (2016) to enforce the volume-filling constraint in mixture models including mass exchanges. The strategy consists in deducing the discrete version of the incompressibility constraint from the discretized mass balance equations. Numerical simulations show that this method constrains the total volume filling constraint, even at the discrete level. Moreover, we add viscous terms in the biofilm model to properly represent biofilms interactions with its fluidic environment. It turns out that a well-balanced numerical scheme becomes of outmost importance to capture the biofilm dynamic when including the viscosity. This modelling upgrade also involves recalibrating model parameters. In particular, the elastic tensors to recover realistic front features. With the new parameters, the numerical set-up becomes more demanding to reach convergence.

INTRODUCTION

There are many physical cases of flows composed of different gas or liquids interacting together. For example, tissue bodies and tumors can be described as a set of interacting viscoelastic materials. Powder-snow avalanches can be described as a mixture of fluid phases. Similarly, the rheology of the gut microbiota and its interactions with chyme (a mixture of partially digested food and water) and the host can be modeled using mixture theory [17]. Complex flows can also be found in many engineering applications involving multiphase systems such as boiling water in nuclear reactors. Therefore, the framework of mixture theory is a common tool to model and study complex flows.

Mathematical models based on mixture theory take the form of systems of partial derivative equations, coupled with algebraic constraints. The theoretical analysis of such systems and the characterization of the qualitative properties of the solutions are extremely complicated [1, 15, 16, 24]. Thus, it is important to develop efficient numerical methods able to accurately capture the solutions [5, 9, 10, 13].

In this article, we are interested in applying mixture models to describe biofilm dynamics. Indeed, mixture theory revealed a powerful approach to represent microbial biofilms where a consortium of cells is embedded in a polymeric structure [10, 19, 20].

<footnote n="1">¹ BIOCORE, Inria, Sorbonne Université, Nice Sophia-Antipolis</footnote>

<footnote n="2">² IRMA, Université de Strasbourg, CNRS UMR-7501, Strasbourg, France</footnote>

<footnote n="3">³ Université Côte d'Azur, Inria, CNRS, LJAD</footnote>

<footnote n="4">⁴ IDP, Université d'Orléans, CNRS, CNRS UMR-7013, Orléans, France</footnote>

<footnote n="5">⁵ Université Côte d'Azur, Inria, CNRS, LJAD</footnote>

<footnote n="6">⁶ IMB, Université de Bordeaux, UMR CNRS 5251, Bordeaux</footnote>

<footnote n="7">⁷ LmB, Université de Franche-Comté, CNRS UMR-6623, Besançon, France</footnote>
In mixture theory, the unknowns of the model are requested to satisfy certain constraints. As far as the continuous equations are considered, several equivalent formulations of these constraints can be derived and used to bring out the properties of the model. However, the preservation of these constraints by a numerical scheme is a challenge and, once a discretization setup has been adopted, it is not clear that all the formulations of the constraints remain equivalent. This issue can induce a loss of stability and accuracy, and eventually a dramatic loss of key physical properties of the simulated flows. Thus, we adapt and extend the numerical scheme proposed in [2] in order to preserve these constraints. The numerical scheme will be tested and illustrated with a multiphasic model representing the development of a photosynthetic biofilm, with the application for biofuel, protein, or drug production.

The paper is organized as follows. The first section is dedicated to the mixture theory framework with a presentation of the simplified model used to test our numerical scheme. The second section details the numerical scheme and its properties. The third section presents the results and comparison with standard numerical schemes.

# 1. MIXTURE THEORY FRAMEWORK: APPLICATION TO BIOFILMS

## 1.1. Mixture theory framework

The mixture theory framework [22], also known as mixture mechanics or continuum mechanics for fluid dynamics, enables describing multi-phasic systems at the mesoscopic scale which is an intermediary scale between microscopic and macroscopic scales. It was introduced in the 1960s by Truesdell [26–28] and generalizes Navier-Stokes equations to multi-phasic systems. The mixture theory framework assumes that each component of the mixture might be present at every point in space and at any time. Moreover, the system’s physical properties (ex. viscosity, incompressibility, ...) are naturally included.

Consider a mixture of $k$ components indexed by $\alpha \in [1, k]$. Each component is locally described by its volumetric mass density $\rho_{\alpha}$, its volume fraction $\phi_{\alpha}$, and its local velocity $v_{\alpha}$. The volume fraction represents the relative volume occupied by a component in an elementary normalized piece of volume. Thus, assuming that there is no vacuum they satisfy the algebraic constraint

$$ \sum_{\alpha=1}^{k} \phi_{\alpha} = 1. \qquad (1) $$

The mixture dynamic depends on mass transfers which are modeled through mass balance equations (2a) and the local forces applied to the system which are accounted for through momentum balance equations (2b). Thus, for each component the state variables satisfy the equations:

$$ \partial_t (\rho_{\alpha} \phi_{\alpha}) + \nabla_x \cdot (\rho_{\alpha} \phi_{\alpha} v_{\alpha}) = \Gamma_{\alpha}, \qquad (2a) $$

$$ \partial_t (\rho_{\alpha} \phi_{\alpha} v_{\alpha}) + \nabla_x \cdot (\rho_{\alpha} \phi_{\alpha} v_{\alpha} \otimes v_{\alpha}) + \nabla_x \pi_{\alpha} + \phi_{\alpha} \nabla_x P = \nabla_x \cdot (\phi_{\alpha} \tau_{\alpha}) + F_{\alpha} + \phi_{\alpha} \rho_{\alpha} g + \Gamma_{\alpha} v_{\alpha}, \qquad (2b) $$

where $\Gamma_i$ is the mass exchange term, $\pi_{\alpha}$ is the elastic tensor, $P$ is the common pressure, $\tau_{\alpha}$ the viscous stress tensor, $F_{\alpha}$ the friction forces, and $g$ the gravity force. Depending on the considered application some forces can be neglected and some others might be added.

Depending on the targeted application one can add for each component an extra equation for the evolution of the density $\rho_{\alpha}$. Nevertheless, liquids are weakly compressible, especially when pressure variations are small. Therefore, in most cases, for liquids the component densities $\rho_{\alpha}$ can be assumed constant. When all the component volumetric mass densities are assumed constant, the mass balance equations (2a) are equivalent to $\partial_t (\phi_{\alpha}) + \nabla_x \cdot (\phi_{\alpha} v_{\alpha}) = \Gamma_{\alpha}/\rho_{\alpha}$. Then summing these equations for each phase leads to the pseudo incompressibility constraint:

$$ \nabla_x \cdot \left( \sum_{\alpha} \phi_{\alpha} v_{\alpha} \right) = \sum_{\alpha} \frac{\Gamma_{\alpha}}{\rho_{\alpha}}. \qquad (3) $$

This means that the local divergence of the averaged mixture velocity is equal to the local volume variation induced by mass exchanges.
The elastic tensor $\pi_{\alpha}$ can be interpreted as the internal pressure of the component. There are several ways to model this term depending on the nature of the component. When the component $\alpha$ represents particles, as in [2], there is a close-packing limit. This property can be enforced by using an appropriate expression for $\pi_{\alpha}$ as

$$ \pi_{\alpha} = \gamma_{\alpha} \frac{\phi_{\alpha}^{\beta_{\alpha}}}{\phi_{\alpha}^{\star} - \phi_{\alpha}}, \quad \text{with} \quad \gamma_{\alpha} > 0, \quad \text{and} \quad \beta_{\alpha} > 1, \qquad (4) $$

where $0 < \phi_{\alpha}^{\star} < 1$ is the so-called close-packing volume fraction limit. When the component $\alpha$ represents softer material like living tissues it can take the form of standard pressure law:

$$ \pi_{\alpha} = \gamma_{\alpha} \left( \frac{\phi_{\alpha}}{\phi_{\alpha}^{\star}} \right)^{\beta_{\alpha}}, \quad \text{with} \quad \gamma_{\alpha} > 0, \quad \text{and} \quad \beta_{\alpha} \ge 1, \qquad (5) $$

where $0 < \phi_{\alpha}^{\star} < 1$ is a threshold, see [10, 21]. More complex laws, based on the Flory–Huggins theory:

$$ \pi_{\alpha} = -\gamma_{\alpha}(\ln(1 - \phi_{\alpha}) + \phi_{\alpha} + \phi_{\alpha}^{2}), \quad \text{with} \quad \gamma_{\alpha} > 0, \qquad (6) $$

enable accounting for colligative properties at low concentrations, see [11].

The viscous stress tensor $\tau_{\alpha}$ is defined by

$$ \tau_{\alpha} = \mu_{\alpha} \phi_{\alpha} \left( \nabla v_{\alpha} + {}^{t}\nabla v_{\alpha} - \frac{2}{3} (\nabla \cdot v_{\alpha}) \mathrm{Id} \right), \qquad (7) $$

where the constant $\mu > 0$ stands for the component dynamic viscosity and ${}^{t}\nabla v_{\alpha}$ stands for the transpose of the velocity differential matrix.

The friction force $F_{\alpha}$ is induced by the difference in the relative speed of the mixture components:

$$ F_{\alpha} = \sum_{\alpha' \neq \alpha} f_{\alpha, \alpha'} (v_{\alpha'} - v_{\alpha}) \qquad (8) $$

with $f_{\alpha, \alpha'}$ the friction force law between the components pair $\alpha$ and $\alpha'$. As a first approximation, it can be assumed that $f_{\alpha, \alpha'}$ is a strictly positive constant. However, the friction between two components should vanish when one of them disappears. Thus, a more realistic alternative is to consider that friction depends on the local composition and use instead $f_{\alpha, \alpha'} (\phi_{\alpha} \phi_{\alpha'})^{\alpha, \alpha'}$. Nevertheless, the total momentum conservation principle enforces that

$$ \sum_{\alpha} F_{\alpha} = 0. $$

Dissolved components, like substrate, can be included. A dissolved component $p$ within a phase $\alpha$ is described through its concentration $\theta_{p}$. In addition to the transport by the phase, it can also diffuse within the phase at a rate $D_{p}$. Thus, the mass balance equations for a dissolved component within the phase $\alpha$ writes:

$$ \partial_t (\rho_{\alpha} \phi_{\alpha} \theta_{p}) + \nabla_x \cdot (\rho_{\alpha} \phi_{\alpha} \theta_{p} v_{\alpha}) - \nabla_x \cdot (\rho_{\alpha} \phi_{\alpha} D_{p} \nabla_x \theta_{p}) = \Gamma_{p}. \qquad (9) $$

where again the source term $\Gamma_{p}$ represents the mass exchange associated to component $p$.

## 1.2. Mixture model for biofilm

We focus on a simplified 1D model for biofilms. Biofilms are made of microorganisms $\mathcal{A}$ (microalgae, bacteria, or a consortium of both) and an extra-cellular matrix $\mathcal{E}$. The biofilm is usually immersed in water $\mathcal{L}$. Therefore, according to mixture theory framework, see section 1.1, each component $\alpha \in \{\mathcal{A}, \mathcal{E}, \mathcal{L}\}$ is described through three macroscopic variables: the mass density $\rho_{\alpha}$, the volume fraction $\phi_{\alpha}$, and the velocity $v_{\alpha}$. By definition, the volume fractions satisfy at any time the algebraic volume-filling constraint (1) which reads in this case: $\phi_{\mathcal{A}} + \phi_{\mathcal{E}} + \phi_{\mathcal{L}} = 1$. In the one-dimensional case, the mass balance equations (2a) writes:

$$ \partial_t (\rho_{\alpha} \phi_{\alpha}) + \partial_x (\rho_{\alpha} \phi_{\alpha} v_{\alpha}) = \Gamma_{\alpha}, \quad \alpha \in \{\mathcal{A}, \mathcal{E}, \mathcal{L}\}. \qquad (10) $$
In this context, the volumetric mass densities $\rho_\alpha$ can be assumed to be constant. Thus, the mixture averaged velocity satisfies the pseudo incompressibility constraint (11) which writes here:

$$ \partial_x (\phi_{\mathcal{A}} v_{\mathcal{A}} + \phi_{\mathcal{E}} v_{\mathcal{E}} + \phi_{\mathcal{L}} v_{\mathcal{L}}) = \frac{\Gamma_{\mathcal{A}}}{\rho_{\mathcal{A}}} + \frac{\Gamma_{\mathcal{E}}}{\rho_{\mathcal{E}}} + \frac{\Gamma_{\mathcal{L}}}{\rho_{\mathcal{L}}}. \quad (11) $$

For biofilms, there are various biological processes to be taken into account. The main processes are growth, extra-cellular matrix excretion, and death. These reactions are schematically represented in Table 1. The parameters $\eta_\alpha$ are pseudo-stoichiometric coefficients that quantify how much a reactant (ex. liquid, algae, substrate, ...) or a product (ex. algae, extra-cellular matrix, ...) is consumed or produced when a reaction occurs. The functions $\psi_i$ are the reaction rates. They describe the speed at which reactions take place as a function of the local composition of the mixture. The source terms read as follows:

$$ \Gamma_{\mathcal{A}} = \psi_g - \psi_e - \psi_d, \quad \Gamma_{\mathcal{E}} = \psi_e + \eta_{\mathcal{E}} \psi_d, \quad \Gamma_{\mathcal{L}} = (1 - \eta_{\mathcal{E}}) \psi_d - \eta_{\mathcal{L}} \psi_g. $$

<table><caption>TABLE 1. Schematic representation of the biochemical reactions considered in the model.</caption><thead><tr><th rowspan="2">Name</th><th colspan="3">Biological reaction representation</th></tr><tr><th>Reactant(s)</th><th>Rate</th><th>Product(s)</th></tr></thead><tbody><tr><td>Growth</td><td>η<sub>L</sub> L + η<sub>S</sub> S</td><td>ψ<sub>g</sub> →</td><td>Α</td></tr><tr><td>Excretion</td><td>Α</td><td>ψ<sub>e</sub> →</td><td>Ε</td></tr><tr><td>Death</td><td>Α</td><td>ψ<sub>d</sub> →</td><td>η<sub>Ε</sub> Ε + (1 - η<sub>Ε</sub>) L</td></tr></tbody></table>

The growth is mainly induced by substrate (S) assimilation and liquid (L) absorption. However, as a first approximation, we assume that the substrate is in excess. Thus, the growth rate $\psi_g$ takes the form $\psi_g = \mu_g \rho_{\mathcal{A}} \phi_{\mathcal{A}} \phi_{\mathcal{L}}$, where $\mu_g$ is the maximal growth rate. The extra-cellular matrix excretion $\psi_e$ and the death rate $\psi_d$ are assumed to be proportional to the quantity of microalgae, thus $\psi_e = \mu_e \rho_{\mathcal{A}} \phi_{\mathcal{A}}$ and $\psi_d = \mu_d \rho_{\mathcal{A}} \phi_{\mathcal{A}}$ respectively. Nevertheless, biofilms are very complex ecosystems and the biological processes are very simplified here. Thus, a model extension accounting for substrate and oxygen is presented in section 4.3.

In the one-dimensional case and neglecting the gravity, for $\alpha \in \{\mathcal{A}, \mathcal{E}, \mathcal{L}\}$ the momentum balance equations simplify into:

$$ \partial_t (\rho_{\alpha} \phi_{\alpha} v_{\alpha}) + \partial_x (\rho_{\alpha} \phi_{\alpha} v_{\alpha}^2) + \partial_x \pi_{\alpha} = -\phi_{\alpha} \partial_x P + \frac{4}{3} \partial_x (\mu_{\alpha} \phi_{\alpha} \partial_x v_{\alpha}) + F_{\alpha} + \Gamma_{\alpha} v_{\alpha}. $$

To keep the model as simple as possible, let us assume that the elastic tensor takes the form of a pressure law, see equation (5), for the tissues (ie. algae and extra-cellular matrix). Since the liquid phase is not elastic this term is null for the liquid, namely $\pi_{\mathcal{L}} = 0$. Similarly, let us assume that the friction forces are constant and symmetric. Thus, in the expression (8) for $F_\alpha$, the term $f_{\alpha,\alpha'}$ for $(\alpha, \alpha') \in \{\mathcal{A}, \mathcal{E}, \mathcal{L}\}^2$ and $\alpha \neq \alpha'$ are constant and such that $f_{\alpha,\alpha'} = f_{\alpha',\alpha}$.

The model is supplemented by boundary conditions. Let $\Omega = [0, L]$ be the domain and $\partial\Omega$ its boundary. In 1D, the domain should correspond to a biofilm core drilling in the orthogonal axis of the support where the biofilm develops. The velocities at the bottom of the domain, which corresponds to the surface on which the biofilm develops, vanish $v_\alpha(0) = 0$, $\alpha \in \{\mathcal{A}, \mathcal{E}, \mathcal{L}\}$. However, the velocity on the top must satisfy a constraint induced by the incompressibility constraint (3). Indeed, the integration over the whole domain of equation (3) combined with the null velocity at the bottom leads to

$$ (\phi_{\mathcal{A}} v_{\mathcal{A}} + \phi_{\mathcal{E}} v_{\mathcal{E}} + \phi_{\mathcal{L}} v_{\mathcal{L}})(x=L) = \int_{0}^{L} \left( \frac{\Gamma_{\mathcal{A}}}{\rho_{\mathcal{A}}} + \frac{\Gamma_{\mathcal{E}}}{\rho_{\mathcal{E}}} + \frac{\Gamma_{\mathcal{L}}}{\rho_{\mathcal{L}}} \right) dx $$

To enforce this condition, let assume that on the top, the velocities are given by $v_{x=L} = \int_{0}^{L} \left( \frac{\Gamma_{\mathcal{A}}}{\rho_{\mathcal{A}}} + \frac{\Gamma_{\mathcal{E}}}{\rho_{\mathcal{E}}} + \frac{\Gamma_{\mathcal{L}}}{\rho_{\mathcal{L}}} \right) dx$.
Remark 1.1. Although there is no biophysical reason to impose the equality between the top velocities, this assumption remains acceptable in this context. Indeed, our focus concerns the biofilm development and the final time considered prevents the biofilm to reach the top of the domain. Therefore, in our context, the hypothesis that all top velocity are equals should not affect the dynamics of the biofilm growth.

### 1.3. Synthesis of model equations

According to the previous section the PDE system under consideration writes:

$$ \phi_{\mathcal{A}} + \phi_{\varepsilon} + \phi_{\mathcal{L}} = 1, \qquad (12a) $$

$$ \partial_t \phi_{\mathcal{A}} + \partial_x (\phi_{\mathcal{A}} v_{\mathcal{A}}) = \frac{\Gamma_{\mathcal{A}}}{\rho_{\mathcal{A}}}, \qquad (12b) $$

$$ \partial_t \phi_{\varepsilon} + \partial_x (\phi_{\varepsilon} v_{\varepsilon}) = \frac{\Gamma_{\varepsilon}}{\rho_{\varepsilon}}, \qquad (12c) $$

$$ \partial_t \phi_{\mathcal{L}} + \partial_x (\phi_{\mathcal{L}} v_{\mathcal{L}}) = \frac{\Gamma_{\mathcal{L}}}{\rho_{\mathcal{L}}}, \qquad (12d) $$

$$ \partial_t (\rho_{\mathcal{A}} \phi_{\mathcal{A}} v_{\mathcal{A}}) + \partial_x (\rho_{\mathcal{A}} \phi_{\mathcal{A}} v_{\mathcal{A}}^2) + \partial_x \pi_{\mathcal{A}} = -\phi_{\mathcal{A}} \partial_x P + \frac{4}{3} \partial_x (\mu_{\mathcal{A}} \phi_{\mathcal{A}} \partial_x v_{\mathcal{A}}) + F_{\mathcal{A}} + \Gamma_{\mathcal{A}} v_{\mathcal{A}}, \qquad (12e) $$

$$ \partial_t (\rho_{\varepsilon} \phi_{\varepsilon} v_{\varepsilon}) + \partial_x (\rho_{\varepsilon} \phi_{\varepsilon} v_{\varepsilon}^2) + \partial_x \pi_{\varepsilon} = -\phi_{\varepsilon} \partial_x P + \frac{4}{3} \partial_x (\mu_{\varepsilon} \phi_{\varepsilon} \partial_x v_{\varepsilon}) + F_{\varepsilon} + \Gamma_{\varepsilon} v_{\varepsilon}, \qquad (12f) $$

$$ \partial_t (\rho_{\mathcal{L}} \phi_{\mathcal{L}} v_{\mathcal{L}}) + \partial_x (\rho_{\mathcal{L}} \phi_{\mathcal{L}} v_{\mathcal{L}}^2) = -\phi_{\mathcal{L}} \partial_x P + \frac{4}{3} \partial_x (\mu_{\mathcal{L}} \phi_{\mathcal{L}} \partial_x v_{\mathcal{L}}) + F_{\mathcal{L}} + \Gamma_{\mathcal{L}} v_{\mathcal{L}}, \qquad (12g) $$

where the sources terms $(\Gamma_{\alpha})_{\alpha}$, the elastic tensors $(\pi_{\alpha})_{\alpha}$ and the drag forces $(F_{\alpha})_{\alpha}$ are given by:

$$ \Gamma_{\mathcal{A}} = \psi_g - \psi_e - \psi_d, \qquad \Gamma_{\varepsilon} = \psi_e + \eta_{\varepsilon} \psi_d, \qquad \Gamma_{\mathcal{L}} = (1 - \eta_{\varepsilon}) \psi_d - \eta_{\mathcal{L}} \psi_g, \qquad (13a) $$

$$ \psi_g = \mu_g \rho_{\mathcal{A}} \phi_{\mathcal{A}} \phi_{\mathcal{L}}, \qquad \psi_e = \mu_e \rho_{\mathcal{A}} \phi_{\mathcal{A}}, \qquad \psi_d = \mu_d \rho_{\mathcal{A}} \phi_{\mathcal{A}}, \qquad (13b) $$

$$ \pi_{\alpha} = \gamma_{\alpha} \left( \frac{\phi_{\alpha}}{\phi_{\alpha}^{*}} \right)^{\beta_{\alpha}}, \qquad \alpha \in \{\mathcal{A}, \mathcal{E}\}, \qquad (13c) $$

$$ F_{\alpha} = \sum_{\alpha' \neq \alpha} f_{\alpha, \alpha'} (v_{\alpha'} - v_{\alpha}), \qquad f_{\alpha, \alpha'} = f_{\alpha', \alpha}, \qquad \alpha \in \{\mathcal{A}, \mathcal{E}, \mathcal{L}\}. \qquad (13d) $$

The system (12) is supplemented with the boundary conditions

$$ v_{\alpha}(x = 0) = 0, \quad \text{and} \quad v_{\alpha}(x = L) = \int_{0}^{L} \left( \frac{\Gamma_{\mathcal{A}}}{\rho_{\mathcal{A}}} + \frac{\Gamma_{\varepsilon}}{\rho_{\varepsilon}} + \frac{\Gamma_{\mathcal{L}}}{\rho_{\mathcal{L}}} \right) dx, $$

for all $\alpha \in \{\mathcal{A}, \mathcal{E}, \mathcal{L}\}$.

The initial data for the volume faction can be chosen arbitrarily provided they are biologically relevant. However, to enforce the algebraic constraint on the sum over all the volume fractions (1), the velocities have to satisfy the incompressibility constraint (11) at all times and therefore the initial velocities must verify this constraint as well. Thus, the initial velocities are defined through a pressure $P$ computed using the incompressibility constraint, see section 2.3, by $v_{\alpha}^0 = \tilde{v}_{\alpha}^0 - \frac{\partial_x P}{\rho_{\alpha}}$ where $\tilde{v}_{\alpha}^0$ is the initial desired velocity. Here, the system is assumed to be initially at rest so $\tilde{v}_{\alpha}^0 = 0$ for all the phases.

Most of the parameters come from [2] or [19, 20]. The viscosity coefficient for microalgae and the extra-cellular-matrix are taken from [18]. All the parameter values are gathered in table 2.

## 2. NUMERICAL SCHEME

In this section, we are interested in the numerical approximation of the PDE system (12). Nevertheless, the general principles and in particular the treatment of the pseudo incompressibility constraint remain valid in a more general context. In such PDE systems, the pressure is defined through the volume filling constraint (1),
<table><caption>TABLE 2. Model parameters. The parameters come from [2, 18, 19].</caption><thead><tr><th>Symbol</th><th>Name</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>μ<sub>g</sub></td><td>Microalgae maximal growth rate</td><td>2</td><td>1/day</td></tr><tr><td>μ<sub>e</sub></td><td>Microalgae maximal ECM excretion rate</td><td>0.4</td><td>1/day</td></tr><tr><td>μ<sub>d</sub></td><td>Microalgae maximal death rate rate</td><td>0.2</td><td>1/day</td></tr><tr><td>μ<sub>r</sub></td><td>Microalgae maximal respiration rate rate</td><td>0.2</td><td>1/day</td></tr><tr><td>η<sub>A</sub></td><td>Microalgae pseudo-stoichiometric coefficient</td><td>1.0</td><td>Ø</td></tr><tr><td>η<sub>L</sub></td><td>Liquid pseudo-stoichiometric coefficient</td><td>0.96</td><td>Ø</td></tr><tr><td>η<sub>S</sub></td><td>Substrate pseudo-stoichiometric coefficient</td><td>8.67 · 10<sup>-2</sup></td><td>Ø</td></tr><tr><td>η<sub>C</sub></td><td>Inorganic carbon pseudo-stoichiometric coefficient</td><td>0.146</td><td>Ø</td></tr><tr><td>η<sub>O</sub></td><td>Oxygen pseudo-stoichiometric coefficient</td><td>0.106</td><td>Ø</td></tr><tr><td>η<sub>E</sub></td><td>Liquid pseudo-stoichiometric coefficient for death</td><td>0.90</td><td>Ø</td></tr><tr><td>K<sub>I</sub></td><td>Light parameter</td><td>0.1</td><td>Ø</td></tr><tr><td>τ</td><td>Light absorption coefficient for the biofilm</td><td>2.5 · 10<sup>4</sup></td><td>m<sup>-1</sup></td></tr><tr><td>I<sub>surf</sub></td><td>Light intensity at the surface</td><td>100</td><td>μmol m<sup>-2</sup>s<sup>-1</sup></td></tr><tr><td>I<sub>opt</sub></td><td>Optimal light intensity</td><td>100</td><td>μmol m<sup>-2</sup>s<sup>-1</sup></td></tr><tr><td>K<sub>S</sub></td><td>Substrate half saturation coefficient</td><td>6.2 · 10<sup>-8</sup></td><td>kg/L</td></tr><tr><td>K<sub>C</sub></td><td>Inorganic carbon half saturation coefficient</td><td>4.4 · 10<sup>-6</sup></td><td>kg/L</td></tr><tr><td>K<sub>O</sub></td><td>Oxygen threshold for growth</td><td>3.2 · 10<sup>-5</sup></td><td>kg/L</td></tr><tr><td>n<sub>O</sub></td><td>Oxygen exponent for growth</td><td>14</td><td>Ø</td></tr><tr><td>K<sub>I</sub></td><td>Light coefficient for Haldane law</td><td>0.1</td><td>Ø</td></tr><tr><td>K<sub>r</sub></td><td>Oxygen half saturation coefficient</td><td>1.0 · 10<sup>-6</sup></td><td>kg/L</td></tr><tr><td>θ<sub>in,S</sub></td><td>Input concentration for substrate</td><td>4 · 10<sup>-5</sup></td><td>kg/L</td></tr><tr><td>θ<sub>in,C</sub></td><td>Input concentration for inorganic carbon</td><td>10 · 10<sup>-5</sup></td><td>kg/L</td></tr><tr><td>θ<sub>in,O</sub></td><td>Input concentration for oxygen</td><td>7.2 · 10<sup>-6</sup></td><td>kg/L</td></tr><tr><td>D<sub>S</sub></td><td>Diffusion coefficient for substrate</td><td>1.47 · 10<sup>-4</sup></td><td>m<sup>2</sup>/day</td></tr><tr><td>D<sub>C</sub></td><td>Diffusion coefficient for inorganic carbon</td><td>1.80 · 10<sup>-4</sup></td><td>m<sup>2</sup>/day</td></tr><tr><td>D<sub>O</sub></td><td>Diffusion coefficient for oxygen</td><td>1.98 · 10<sup>-4</sup></td><td>m<sup>2</sup>/day</td></tr><tr><td>ρ<sub>A</sub></td><td>Microalgae volumetric mass density</td><td>1050</td><td>kg/m<sup>3</sup></td></tr><tr><td>ρ<sub>C</sub></td><td>Extra-cellular matrix volumetric mass density</td><td>1050</td><td>kg/m<sup>3</sup></td></tr><tr><td>ρ<sub>L</sub></td><td>Liquid volumetric mass density</td><td>1025</td><td>kg/m<sup>3</sup></td></tr><tr><td>φ<sup>*</sup><sub>A</sub></td><td>Microalgae close packing threshold</td><td>0.75</td><td>Ø</td></tr><tr><td>γ<sub>A</sub></td><td>Microalgae viscoelastic tensor coefficient</td><td>1.2 · 10<sup>-9</sup></td><td>kg m<sup>-1</sup>day<sup>-1</sup></td></tr><tr><td>β<sub>A</sub></td><td>Microalgae viscoelastic tensor exponent</td><td>1</td><td>Ø</td></tr><tr><td>φ<sup>*</sup><sub>E</sub></td><td>Extra-cellular matrix close packing threshold</td><td>0.75</td><td>Ø</td></tr><tr><td>γ<sub>E</sub></td><td>Extra-cellular matrix viscoelastic tensor coefficient</td><td>1.2 · 10<sup>-9</sup></td><td>kg m<sup>-1</sup>day<sup>-1</sup></td></tr><tr><td>β<sub>E</sub></td><td>Extra-cellular matrix viscoelastic tensor exponent</td><td>1</td><td>Ø</td></tr><tr><td>μ<sub>L</sub></td><td>Liquid viscosity</td><td>10<sup>-3</sup></td><td>Pa s</td></tr><tr><td>μ<sub>A</sub></td><td>Microalgae viscosity</td><td>0.25</td><td>Pa s</td></tr><tr><td>μ<sub>E</sub></td><td>Extra-cellular matrix viscosity</td><td>0.75</td><td>Pa s</td></tr><tr><td>f<sub>A,E</sub></td><td>Friction coefficient between A and E</td><td>20</td><td>kg m<sup>-3</sup>day<sup>-1</sup></td></tr><tr><td>f<sub>A,L</sub></td><td>Friction coefficient between A and L</td><td>20</td><td>kg m<sup>-3</sup>day<sup>-1</sup></td></tr><tr><td>f<sub>E,L</sub></td><td>Friction coefficient between E and L</td><td>20</td><td>kg m<sup>-3</sup>day<sup>-1</sup></td></tr></tbody></table>
namely $\phi_A + \phi_E + \phi_L = 1$ for the considered model. The treatment of this constraint and thus the definition of the pressure is always an issue and requires specific treatment. To this end, the momentum equations are treated using a projection correction method inspired by the numerical method introduced by Chorin [6–8] and Temam [25] for incompressible viscous flows. In a nutshell, the momentum equation is decomposed using a time splitting to separate the contribution of the pressure as follows:

$$ \partial_t (\rho_\alpha \phi_\alpha v_\alpha) + \partial_x (\rho_\alpha \phi_\alpha v_\alpha^2) + \partial_x \pi_\alpha = \frac{4}{3} \partial_x (\mu_\alpha \phi_\alpha \partial_x v_\alpha) + F_\alpha + \Gamma_\alpha v_\alpha, \quad (14a) $$

$$ \partial_t (\rho_\alpha \phi_\alpha v_\alpha) + \phi_\alpha \partial_x P = 0. \quad (14b) $$

for $\alpha \in \{\mathcal{A}, \mathcal{E}, \mathcal{L}\}$.

## 2.1. Projection correction method

Let us start with the presentation of the time discretization. Let $T \in \mathbb{R}^+$ be the final time and $(t_n)_{n \ge 0}$ a subdivision of $[0, T]$ such that $t_n = \sum_{k=0}^n \Delta t_k$. Consider $\alpha \in \{\mathcal{A}, \mathcal{E}, \mathcal{L}\}$ a phase and its associated volume fraction $\phi_\alpha$ and velocity $v_\alpha$. Then, $\phi_\alpha^n(x)$ and $v_\alpha^n(x)$ denote, respectively, their approximation at time $t_n$. To shorten the notations, let us drop the space variable $x$ and denote $\delta t = \Delta t_{n+1}$. Assuming that all the quantities are known at time $t_n$, the approximated solution at time $t_{n+1} = t_n + \delta t$ is computed using the following steps:

(1) Update the volume fractions according to the mass balance equations (12b)-(12d):

$$ \phi_{\alpha}^{n+1} = \phi_{\alpha}^{n} - \frac{\delta t}{\rho_{\alpha}} \partial_{x} (\phi_{\alpha}^{n} v_{\alpha}^{n}) + \frac{\delta t}{\rho_{\alpha}} \Gamma_{\alpha}^{n}. \quad (15) $$

(2) Update the momentum equations without the contribution of the pressure term by solving the following system:

$$ \phi_{\alpha}^{n+1} v_{\alpha}^{n+\frac{1}{2}} - \phi_{\alpha}^{n} v_{\alpha}^{n} = \frac{\delta t}{\rho_{\alpha}} \left( -\partial_x (\phi_{\alpha}^{n} (v_{\alpha}^{n})^2) - \partial_x \pi_{\alpha}^{n} + \frac{4}{3} \partial_x \left( \mu_{\alpha} \phi_{\alpha}^{n+1} \partial_x v_{\alpha}^{n+\frac{1}{2}} \right) + F_{\alpha} (\phi^{n+1}, v^{n+\frac{1}{2}}) + \Gamma_{\alpha}^{n} v_{\alpha}^{n} \right). \quad (16) $$

(3) Compute the pressure using the incompressibility constraint (3). This step is detailed in subsection 2.3.

(4) Update the velocity using the pressure with:

$$ v_{\alpha}^{n+1} = v_{\alpha}^{n+\frac{1}{2}} - \frac{\delta t}{\rho_{\alpha}} \partial_x P^{n+1}. $$

## 2.2. 1D space discretization

Following [2], the space is discretized using staggered grids. This enables avoidance of any odd/even decoupling in the stencil of the discrete version of the system. Moreover, the use of staggered grids also allows to have or deduce naturally the quantity of interest (e. g. deduce the pressure gradient on the velocity mesh grid). Let $(x_i)_{i \in \llbracket 0, I \rrbracket}$ be a regular subdivision of the domain $\Omega$ such that $x_i = i\Delta x$ with $\Delta x = L/I$ the mesh step. Let also define the mesh cell centers: $x_{i+\frac{1}{2}} = (i + \frac{1}{2})\Delta x$ for $i \in \llbracket 0, I - 1 \rrbracket$. The model variables are located:

* at the mesh cell centers for the volume fraction and the pressure: $\phi_{\alpha, i+\frac{1}{2}}$, $P_{i+\frac{1}{2}}$ for $0 \le i \le I-1$

* at the mesh cell edges for the velocities: $v_{\alpha, i}$ for $0 \le i \le I$.

Figure 1 gives an example of the staggered grids with the localization of model variables.

Model unknowns are discretized using a finite volume scheme. The transport terms in the mass balance equations (12b)-(12d) are written:

$$ \phi_{\alpha, i+\frac{1}{2}}^{n+1} = \phi_{\alpha, i+\frac{1}{2}}^{n} - \frac{\delta t}{\rho_{\alpha} \Delta x} \left( \mathcal{F}_{i+1}(\phi_{\alpha}^{n}, v_{\alpha}^{n}) - \mathcal{F}_{i}(\phi_{\alpha}^{n}, v_{\alpha}^{n}) \right) + \frac{\delta t}{\rho_{\alpha}} \Gamma_{\alpha, i+\frac{1}{2}} \quad (17) $$

where $\mathcal{F}_i$ represents the numerical mass flux at the interface $x_i$, which is a function of the neighboring cells. There are multiple relevant choices for the definition of the numerical flux. For the sake of simplicity, to ensure
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Regular staggered grid in one dimension for 5 mesh cells with the volume fractions and the velocities locations. The pressure $P$ and the phase volume fractions $(\phi_{\alpha})_{\alpha}$ are located at the mesh cell centers $(x_{i+\frac{1}{2}})_{0\le i\le 4}$.</figcaption></figure>

stability and since it is well adapted to staggered grids, it is convenient to use upwind numerical flux. Thus, the discrete mass flux is defined by $\mathcal{F}_i(\phi, v) = \mathcal{F}^+(\phi_{i-\frac{1}{2}}, v_i) + \mathcal{F}_i^-(\phi_{i+\frac{1}{2}}, v_i)$ with

$$ \mathcal{F}^{+}(\phi, v) = \begin{cases} 0 & \text{if } v \le 0, \\ \phi v & \text{if } v > 0, \end{cases} \quad \text{and} \quad \mathcal{F}^{-}(\phi, v) = \begin{cases} \phi v & \text{if } v < 0, \\ 0 & \text{if } v \ge 0. \end{cases} $$

All the volume fractions are updated using equation (17). Thus, the volume-filling constraint enforcement is not guaranteed and depends on the strategy used to compute the pressure, see sections 2.3 and 3.2.

**Remark 2.1.** To update the volume fractions and ensure volume-filling constraint enforcement another strategy consists to use equation (17) for all the components except one (usually the liquid) which is computed using the algebraic volume-filling constraint (1): $\phi_{\alpha'} = 1 - \sum_{\alpha \neq \alpha'} \phi_{\alpha}$ as done in [9, 10, 19, 20].

For the momentum balance equation, following [2], the transport term is also discretized using an upwind strategy based on the material velocity $v$, that is the momentum flux is defined by

$$ \mathcal{G}_{i+\frac{1}{2}} = \frac{v_{\alpha,i}^n}{2} \left( \mathcal{F}^+(\phi_{i-\frac{1}{2}}, v_i) + \mathcal{F}^+(\phi_{i+\frac{1}{2}}, v_{i+1}) \right) + \frac{v_{\alpha,i+1}^n}{2} \left( \mathcal{F}^-(\phi_{i+\frac{1}{2}}, v_i) + \mathcal{F}^-(\phi_{i+\frac{1}{2}}, v_{i+1}) \right). $$

The other terms of equation (16) are discretized using standard approximations. Remark that interpolation on the dual mesh is required only for the zeroth order terms like the momentum supply induced by mass exchanges or friction forces. For these terms, the approximation of the volume fraction on the dual mesh is obtained by approximating the volume fractions using the values in the neighboring cells: $\phi_i = \frac{1}{2}(\phi_{i-\frac{1}{2}} + \phi_{i+\frac{1}{2}})$. Therefore, dropping the $\alpha$ for readability, equation (16) is discretized as follows:

$$ \phi_i^{n+1} v_i^{n+\frac{1}{2}} - \frac{4\delta t}{3\Delta x^2} \mu \left( \phi_{i+\frac{1}{2}}^{n+1} v_{i+1}^{n+\frac{1}{2}} - \left( \phi_{i+\frac{1}{2}}^{n+1} + \phi_{i-\frac{1}{2}}^{n+1} \right) v_i^{n+\frac{1}{2}} + \phi_{i-\frac{1}{2}}^{n+1} v_{i-1}^{n+\frac{1}{2}} \right) - \delta t F \left( \phi_i^{n+1}, v_i^{n+\frac{1}{2}} \right) \\ = \phi_i^n v_i^n - \frac{\delta t}{\Delta x \rho} \left( \mathcal{G}_{i+\frac{1}{2}} + \pi \left( \phi_{i+\frac{1}{2}}^n \right) - \mathcal{G}_{i-\frac{1}{2}} - \pi \left( \phi_{i-\frac{1}{2}}^n \right) \right) + \delta t \Gamma(\phi_i^n) v_i^n. $$

**Remark 2.2.** In this projection step, the viscosity and the friction are treated implicitly. For the viscosity, this treatment enables the relaxation of the CFL constraint and avoids numerical instabilities.

**Remark 2.3.** Like in [10, 19, 20] the computation of friction forces requires a specific treatment. Indeed, the friction forces depend on the difference between the phase velocities, and when a phase vanishes the velocity can not be deduced from the momentum (ie. $\phi v$). In the considered applications, areas of pure liquid or biofilm are important so the adaptation of the initial data to avoid phase vanishing is irrelevant. To overcome this difficulty, a strategy consists to treat these terms implicitly so the velocity can be directly computed using the above equation. However, this is costly because it imposes to solve at each time step a linear system of size: number of phases×mesh grid size.

Finally, the space discretization of the correction step is given by: $v_{\alpha}^{n+1} = v_{\alpha}^{n+\frac{1}{2}} - \frac{\delta t}{\rho_{\alpha} \Delta x} \left( P_{i+\frac{1}{2}}^{n+1} - P_{i-\frac{1}{2}}^{n+1} \right)$
## 2.3. Pressure approximation

Let us detail the third step of the projection correction method. This is the key step to enforce the algebraic constraint on the sum over all the volume fractions (1). The standard strategy consists in plugging the time discrete version of equation (14b): $\phi_{\alpha}^{n+1} v_{\alpha}^{n+1} = \phi_{\alpha}^{n+1} v_{\alpha}^{n+\frac{1}{2}} - \frac{\delta t}{\rho_{\alpha}} \phi_{\alpha}^{n+1} \partial_x P^{n+1}$ into the incompressibility constraint (11) to obtain the following equation on the pressure:

$$ \partial_x \left( \sum_{\alpha} \phi_{\alpha}^{n+1} v_{\alpha}^{n+\frac{1}{2}} - \frac{\delta t}{\rho_{\alpha}} \phi_{\alpha}^{n+1} \partial_x P^{n+1} \right) = \sum_{\alpha} \frac{\Gamma_{\alpha}}{\rho_{\alpha}} \qquad (18) $$

Thus the pressure can be obtained by solving a non-linear and inhomogeneous Poisson equation. As mentioned above, this strategy relies on the use of the continuous version of the incompressibility constraint. Therefore, there is no guarantee that the algebraic volume-filling constraint will be fulfilled at the discrete level.

To enforce the algebraic volume filling constraint, we adapt the strategy proposed in [2], which consists in using the fully discretized mass balance equations to deduce the appropriate discrete incompressibility constraint. To this end, let us assume that the constraint $\sum_{\alpha} \phi_{\alpha}^{n} = 1$ is satisfied for all times $(t_n)_{n \ge 0}$ and in all the grid mesh cells. Thus, the sum of the equations (17) over the phases leads to

$$ \sum_{\alpha} \frac{1}{\rho_{\alpha}} (\mathcal{F}_{i+1}(\phi_{\alpha}^{n}, v_{\alpha}^{n}) - \mathcal{F}_{i}(\phi_{\alpha}^{n}, v_{\alpha}^{n})) = \Delta x \sum_{\alpha} \frac{\Gamma_{\alpha, i+\frac{1}{2}}}{\rho_{\alpha}}. \qquad (19) $$

Then, as in the standard strategy, an equation on the pressure or its gradient can be deduced by using the time discrete version of equation (14b). Since in the correction step, the volume fractions remain unchanged, the time discrete version of equation (14b) simplifies into $v_{\alpha}^{n} = v_{\alpha}^{n-\frac{1}{2}} - \frac{\delta t}{\rho_{\alpha}} \partial_x P^n$. Injecting this relation into equation (19) gives:

$$ \sum_{\alpha} \frac{1}{\rho_{\alpha}} \left( \mathcal{F}_{i+1} \left( \phi_{\alpha}^{n}, v_{\alpha}^{n-\frac{1}{2}} - \frac{\delta t}{\rho_{\alpha}} \partial_x P^n \right) - \mathcal{F}_{i} \left( \phi_{\alpha}^{n}, v_{\alpha}^{n-\frac{1}{2}} - \frac{\delta t}{\rho_{\alpha}} \partial_x P^n \right) \right) = \Delta x \sum_{\alpha} \frac{\Gamma_{\alpha, i+\frac{1}{2}}}{\rho_{\alpha}}. \qquad (20) $$

Consequently, defining the pressure as a solution of the non-linear equation (20) ensures that (19) holds. Next, going back to the mass balance equations (17), it guarantees that $\sum_{\alpha} \phi_{\alpha, i+\frac{1}{2}}^{n+1} = \sum_{\alpha} \phi_{\alpha, i+\frac{1}{2}}^{n}$, so that the algebraic volume-filling constraint is met at time $t_{n+1}$ if it holds at time $t_n$. The existence of such a pressure, solution of (20), can be obtained through topological degree arguments as in [2, Annex B]. The solution can be approximated using Newton's methods. In practice, although this method is more expensive than the standard approach its cost remains reasonable. Indeed, the Jacobian matrix is explicitly known and the solution at the previous time step reveals to be a good initial guess so only very few iterations are necessary to converge. Both strategies are compared in subsection 3.2.

## 3. NUMERICAL RESULTS

The paper aims to present and test a numerical method able to simulate mixture models for biofilms by guaranteeing the preservation of the algebraic volume filling constraint. Another challenge when one wants to go towards the applications, relies on the difficulty of calibrating the parameters of the model. Many parameters are, up to our knowledge, not available in the current literature and very difficult to extrapolate from experimental data. For example, in [9, 10, 19, 20] the elastic tensor (ie. $\pi_{\alpha}$) settings are calibrated so that the biofilm front velocity matches observations. Consequently, any modification of the model requires recalibration. To avoid such difficulties, subsection 3.1 presents numerical simulations based on the numerical scheme presented in section 2, but assumes that the viscosity can be neglected, which enables reusing parameters from [9, 10, 19, 20] for the elastic tensors. Secondly, subsection 3.2 presents comparisons between the two strategies to approximate the pressure, still neglecting the viscosity. Finally, subsection 4.1 presents the dynamic of the full model including viscosity and recalibration of the elastic tensors.
Initially, the mixture is only made of microalgae and liquid and the volume fractions are set by

$$ \phi_{\mathcal{A}}^{0}(x) = \max\{0, 0.05(x - 0.1)(x + 0.1)\}, \quad \phi_{\mathcal{E}}^{0} = 0, \quad \text{and} \quad \phi_{\mathcal{L}}^{0} = 1 - \phi_{\mathcal{A}}. \qquad (21) $$

As mentioned in subsection 1.3, the system is assumed to be at rest. Thus, the initial velocities are defined by $v_{\alpha}^{0} = -\frac{1}{\rho_{\alpha}}\partial_{x}P$ where the pressure $P$ is determined according to the strategy presented in subsection 2.3 to enforce the algebraic constraint on the sum over all the volume fractions (1).

Unless otherwise mentioned, all the simulations presented hereafter are made on a domain of length $L = 10mm$ so $\Omega = [0, 10]$ and using 512 mesh cells, namely $I = 512$. The Courant–Friedrichs–Lewy (CFL) constant used for the hyperbolic part is $\frac{1}{4}$. Here, the velocities are very low thus the CFL condition is not binding. However, there are also stability conditions associated with the explicit treatment of the source terms for negative contributions. These negative contributions are proportional to volume fractions, thus a CFL condition could be derived from (15). In practice, the maximal time step is fixed to $10^{-1}$.

### 3.1. Biofilm dynamic without viscosity

Figure 2 presents the numerical results for different times of system (12) where the viscous terms are neglected, namely $\mu_{\alpha} = 0$ for $\alpha \in \{\mathcal{A}, \mathcal{E}, \mathcal{L}\}$. The presence of non-conservative terms in the equations is a common feature of many systems describing multiphase flows. This leads to the issue of deciding whether such terms are meaningful when the solutions present discontinuities, which is likely the case, especially when the viscosity is neglected. This issue has been addressed, to our knowledge mostly in the case of compressible flows, where the mass densities of the constituents are not constant. In such cases, the pressure terms are defined by state laws and a key point is to check whether or not the system keeps the hyperbolic property that guarantees some kind of stability of the model, which eventually is an ingredient to prove the existence of solutions [3]. Then, the definition of the non-conservative products comes from the analysis of the jump relations and the structure of the discontinuous solutions [12, 14]. Several models of the interface pressure have been introduced that meet these requirements. Such a discussion is beyond the scope of the present work, and we refer the interested reader to the review [4]. Here, we have performed several runs of simulations, by varying time and space steps as well as CFL number to assess the robustness of the presented results.

The simulation is made using the numerical scheme presented in section 2 and using the strategy based on the adaptation of [2] for the computation of the pressure, see subsection 2.3. In these figures, the left side corresponds to the surface where the biofilm sticks and develops and the right side corresponds to the side covered by the liquid, where nutrients are brought.

According to Figures 2a, 2b and 2c, there is front propagation corresponding to the biofilm (dashed orange curve) development within the liquid. As in [19], two areas can be distinguished within the biofilm. For example in Figure 2c, on the left side, namely for $x \in [0, \sim 6]\mu m$, the biofilm is mainly made of extra-cellular matrix (ie. $\mathcal{E}$), whereas on the biofilm front, namely for $x \in [\sim 6, 7.6]\mu m$, the biofilm is mainly made of microalgae. On the opposite, the right side, namely for $x > 7.6$, is made of pure water.

For the biofilm components, the velocities are positive near the front, which is expected and explains the biofilm expansion. Otherwise, the liquid velocity is negative in the biofilm region, which means that the liquid is drained into the biofilm due to its consumption for the biofilm growth.

### 3.2. Volume filling constraint validation

Let us compare the two strategies presented in subsection 2.3 to enforce the algebraic volume-filling constraint, that is computing the pressure $P$ either as the solution of discretization the linear equation (18) (standard strategy), or as the solution of the non-linear equation (20) (adapted strategy). To this aim, as in subsection 3.1, the system (12) without the viscous terms is simulated, but using the standard strategy to enforce the volume filling constraint. The results for the mixture components volume fractions at $t = 360h$ are presented in Figure 3a and can be compared to Figure 2c. According to these plots, the results are comparable. Similarly, the shape of the pressure gradients curves are also similar, see Figure 3b. Nevertheless, according to the dotted purple curve in Figure 3b, there is a significant discrepancy close to the biofilm front (ie. at $x = 7.6mm$) in the pressure gradients between the two correction strategies. Note that the front is the active part of the biofilm. Namely, it is in this
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Mixture components volume fractions (first row) and velocities (second row) for different times.</figcaption></figure>

area that the source terms are the largest and lead to significant changes in mixture composition. Therefore, it is expected that the effect of the pressure gradient is observable notably there.

Besides, the pressure can be interpreted as the Lagrange multiplier associated with the volume-filling constraint. Thus, it is important to compare how these strategies enable enforcing at the discrete level the volume filling constraint (1). To this end, Figure 3c represents the sum of the volume fractions within the domain a time $t = 360h$. According to this plot, the strategy adapted from [2] enables ensuring the volume filling constraint, whereas the standard strategy does not. Numerically, the maximal error on the volume-filling constraint for the standard strategy is $1.007 \cdot 10^{-3}$, whereas with the adapted strategy, it is $5.107 \cdot 10^{-15}$, namely the order of magnitude of the precision used in Newton's method. Moreover, with this adapted method, the error remains negligible throughout the simulation whereas, with the standard strategy, it varies over time, see Figure 7 in Appendix B.

*Remark 3.1.* The strategy adapted from [2] requires the resolution of a non-linear equation which is here solved using Newton's method. As already mentioned, although this method is more expensive than the standard approach its cost remains reasonable. Indeed, the Jacobian matrix is explicitly known and the solution at the previous time step reveals to be a good initial guess so only very few iterations are necessary to converge. Typically, two or three iterations are required thus the pressure computation is more or less twice more costly.

# 4. MODEL EXTENSIONS

Following insights coming from [19], this section presents various relevant extensions of the model and their numerical simulations.

## 4.1. Including the viscosity

Adding the viscous terms for the components requires recalibrating the model parameters. Indeed, the viscosity is a measure of the component's resistance to deformation. Therefore, when accounting for the viscosity, the parameters associated with the component's ability to deform must be adapted. In particular, the elastic tensors for the microalgae and the extra-cellular matrix must be recalibrated. Moreover, up to our knowledge, there is no
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Mixture components volume fractions (left), pressure gradients (center) and the sum of the volume fractions (right) at time $t = 360h$. In figures 3b and 3c the blue curve represents correspond to a simulation made using the standard strategy for the pressure gradient computation (ie. solving a Poisson equation) and the dashed orange curve represents the results obtained using the strategy adapted from [2] which require the resolution of a non-linear equation (ie. Newton's method).</figcaption></figure>

direct measurement of the parameters and they are calibrated, see [9, 10, 19], such that the biofilm front velocity matches experimental measurements, see [23]. However, such calibration is extremely complex because the biofilm front velocity depends also on many other parameters like the growth or death rate. Nevertheless, to get the right order of magnitude of the biofilm front velocity the elastic tensor coefficients must be significantly increased: multiplied by $9 \cdot 10^7$ so is set to $\gamma_A = \gamma_E = 4.5 \cdot 10^{-3} \text{kg/m/day}$.

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Mixture components volume fractions (first row) and velocities (second row) for different times. In this simulation, the viscosity is included and the elastic tensors for the biological phases are multiplied $9 \cdot 10^7$ so the biofilm front velocity matches experimental measurements. The simulation is made using 2048 mesh cells for the space grid.</figcaption></figure>

Figure 4 represents the time dynamic of mixture components when accounting for the contribution of viscosity. The global dynamic remains comparable to the dynamic observed in Figure 2. In particular, there is still a biofilm
traveling front. Again, there are two areas within the biofilm: the back which is mainly made of an extra-cellular matrix (for $x \in [0, 5]$ in subfigure 4c), and the front which is mainly made of microalgae (for $x \in [3.97, 6.7]\mu m$ in subfigure 6c). Nevertheless, a major discrepancy is that the microalgae remain more located at the front when including the viscosity. This is particularly visible at $t = 240h$ when comparing Figure 2b and Figure 4b. In addition, at $t = 240h$, the velocities order of magnitude close to the front is larger when including the viscosity. However, the interpretation of this observation is tricky. Indeed, the shift in the elastic tensors for the biological phases imposes the use of very refined mesh grids to properly capture the biofilm dynamic. Thus, it would be of particular interest to design and use a well-balanced numerical scheme able to preserve the biofilm front structure. For more details about the numerical convergence of the scheme, see appendix C.

## 4.2. Including light intensity

A microalga is a photosynthetic organism. Thus, microalgae require light to grow. When microalgae develop within a biofilm, the upper layers overshadow the lower layers. Following [9, 10, 19] to account for these mechanisms, the microalgae growth rate becomes $\psi_g = \mu_g \rho_A \phi_A \phi_L f_I$, where $f_I$ accounts for the effect of light on growth. This term depends on the rescaled received light intensity $I$ and takes the form of the Haldane law:

$$f_I = \frac{2(1 + K_I)I}{I^2 + 2K_I I + 1}. \qquad (22)$$

The rescaled light intensity is the ratio between the received light and the optimal light intensity $I_{opt}$, namely:

$$I(t, x) = \frac{I_{\text{surf}}}{I_{\text{opt}}} \exp\left(-\int_{x}^{L} \tau(1 - \phi_{\mathcal{L}}(t, y))dy\right), \qquad (23)$$

where $I_{surf}$ is the light intensity at the surface of the tank (ie. $x = L$) and $\tau$ the attenuation coefficient of the biofilm, assuming that microalgae and extra-cellular matrix have the same attenuation rate. The parameter values associated to the light are gathered in Table 2.

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Mixture components volume fractions for different times. In this simulation, the effective microalgae growth rate (ie. $\psi_g$) accounts for the contribution of light intensity through Haldane's law (22) and light attenuation induced by biofilm layers, see equation (23)</figcaption></figure>

Figure 5 represents the time dynamic of mixture components when accounting for the contribution of light. The global dynamic is comparable to the dynamic observed in Figure 2. In particular, the biofilm front position travels at the same speed, and, again, there are two areas within the biofilm: the back which is mainly made of an extra-cellular matrix (for $x \in [0, 6]$ in subfigure 5c) and the front which is mainly made of microalgae (for $x \in [6, 7.6]\mu m$ in subfigure 5c). However, as expected, the volume fraction of biofilm is lower. Indeed, taking into account the effect of light reduces the growth in the shadowed areas and thus the biomass of microalgae. The extra-cellular matrix is also reduced since it is made from microalgae excretion and dead organisms.
### 4.3. Including light intensity and solutes

Following [19], let us now include three dissolved components: the substrate (S), the inorganic carbon (C), and the oxygen (O). As mentioned in subsection 1.1, the dynamic for dissolved components is modeled using a convection-diffusion reaction equation (9).

In a nutshell, the substrate represents the nitrate which is a nutrient of primary importance for the growth of autotrophic organisms like microalgae. Besides, roughly speaking, photosynthesis is the assimilation of inorganic carbon using light energy by autotrophic organisms. Photosynthesis releases oxygen. Thus, including these components is of primary interest. Taking into account these compounds also allows us to include the process of respiration. Basically, respiration is the opposite mechanism of photosynthesis and its consideration allows us to better describe the dynamic of thick biofilms. Indeed, the process of respiration becomes non-negligible in the absence of light, namely in the biofilm's inner layers.

As for the light, the contributions of the dissolved components to the photosynthesis process are accounted for in the growth through the multiplication by functions $f_p, p \in \{S, C, O\}$ which represent how the growth is modified by the local concentration of these components. On the one hand, limited contribution in high concentration regimes of the substrate and the inorganic carbon is modeled using Monod's law: $f_p = \frac{\theta_p}{K_p+\theta_p}$. On the other hand, the inhibition induced by high oxygen concentration is modeled by the sigmoidal function $f_O = \frac{1}{1+\left(\frac{\theta_O}{K_O}\right)^{n_O}}$. Thus, including the contribution of the dissolved components and the light intensity, the algae growth rate becomes: $\psi_g = \mu_g \rho_A \phi_A \phi_L f_I f_S f_C f_O$.

The respiration process is modeled by $\psi_r = \mu_r \phi_A \frac{\theta_O}{K_r+\theta_O}$ where $\mu_r$ is the maximal respiration rate and $K_r$ the half-saturation constant for the oxygen.

The modification of the microalgae growth rate and the inclusion of the respiration process requires to adapt the source terms for the phases as follows:

$$ \Gamma_A = \psi_g - \psi_e - \psi_d - \psi_r, $$

$$ \Gamma_E = \psi_e + \eta_E \psi_d, $$

$$ \Gamma_L = (1 - \eta_E) \psi_d + \eta_L (\psi_r - \psi_g). $$

As for a phase, the source term for a dissolved component is the sum of the pseudo-stoichiometric coefficients multiplied by the reaction rates. Thus, for the dissolved components, the source terms are

$$ \Gamma_S = -\eta_S \psi_g, $$

$$ \Gamma_C = -\eta_g^C \psi_g + \eta_r^C \psi_r, $$

$$ \Gamma_O = \eta_g^O \psi_g - \eta_r^O \psi_r. $$

The external supply for the dissolved components is modeled through Dirichlet boundary conditions at the top of the bioreactor, namely at $x = L$. Otherwise, the no flux boundary condition at the bottom of the bioreactor is modeled using the Neumann boundary condition: $\partial_x \theta_p|_{x=0} = 0$ for $p \in \{S, C, O\}$.

The parameter values associated with the inclusion of the dissolved components are gathered in Table 2.

Numerically, the transport and reaction terms in the mass balance equations for the solutes are treated similarly to the other components. The diffusion terms are treated implicitly to ensure stability without constraining the CFL condition.

Figure 6 represents the time dynamic of mixture components when accounting for the contribution of light and solutes. The global dynamic is comparable to the dynamic observed in Figure 5. In particular, the biofilm front position travels at a comparable speed. Again there are two areas within the biofilm: the back which is mainly made of an extra-cellular matrix (for $x \in [0, 6.5]$ in subfigure 6c), and the front which is mainly made of microalgae (for $x \in [5.5, 6.6]\mu m$ in subfigure 6c). However, the biofilm front velocity is slightly slower here. This can be explained by the fact that the lack or excess of solutes in the active part of the biofilm slightly reduces its growth. Indeed, for example, at $t = 360h$, within the biofilm area, the concentration of substrate is reduced by 12.4% and the concentration of inorganic carbon is reduced by 5.1% relatively to the input values (ie. $\theta_{in}$). Besides, the concentration of oxygen is increased by 46.3% relatively to $\theta_{in,C}$. These discrepancies are larger at the beginning and tend to decrease over time, see Figure 9 in the supplementary material. These results are in good agreement with the results presented in [19, 20].
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Mixture components volume fractions (first row) and velocities (second row) for different times.</figcaption></figure>

## 5. CONCLUSIONS AND PERSPECTIVES

This article proposes an adaptation of the numerical scheme presented in [2] able to enforce the volume filling constraint in mixture models including mass exchanges. As in [2] the strategy consists in deducing the discrete version of the incompressibility constraint from the discretized mass balance equations. Numerical simulations show that this method enables the enforcement of the total volume filling constraint at the discrete level.

In addition, on the modeling side, previous models from the literature are enriched by the inclusion of viscous terms. These terms are essential to properly model biofilms in their fluidic environment especially when there is a mixing of the surrounding fluid. In this context, this work has allowed us to highlight the importance of designing well-balanced numerical scheme able to efficiently capture the biofilm dynamic when including the viscosity. Indeed, including the viscosity requires recalibrating model parameters; in particular the elastic tensors need to be strongly rescaled to recover realistic front features. However, with these parameters, the numerical set-up is more demanding to reach convergence. This difficulty leads to consider further the design of a specific well-balanced scheme for the problem. To this end, the use of well-balanced numerical schemes able to preserve the equilibrium at the biofilm front can be considered.

Finally, in further works, it would be interesting to include additional biological features. Among others, biofilms are generally multi-species. The framework of mixture theory is well adapted to incorporate different species and such extensions are affordable if the interaction between the species and their metabolisms is known. To make the model even more realistic and predictive its calibration on experimental data is also particularly interesting. In conclusion, real-life biofilms are 3D and therefore the extension and implementation of the numerical method in 2D and 3D should be considered.

## Acknowledgement

* This work was supported by Conseil Regional de Bourgogne Franche-Comté (France).
REFERENCES

[1] H. Beirão da Veiga. Diffusion on viscous fluids. existence and asymptotic properties of solutions. *Annali della Scuola Normale Superiore di Pisa-Classe di Scienze*, 10(2):341-355, 1983.

[2] F. Berthelin, T. Goudon, and S. Minjeaud. Multifluid flows: A kinetic approach. *Journal of Scientific Computing*, 66(2):792-824, 2016.

[3] D. Bresch, B. Desjardins, J.-M. Ghidaglia, and E. Grenier. Global weak solutions to a generic two-fluid model. *Arch. Rational Mech. Anal.*, 196:599-629, 2010.

[4] D. Bresch, B. Desjardins, M. Ghidaglia, J. E. Grenier, and M. Hillairet. Multi-fluid models including compressible fluids. In *Handbook of mathematical analysis in mechanics of viscous fluids*, pages 2927-2978. Springer, Cham, 2018.

[5] C. Calgaro, E. Creusé, and T. Goudon. Modeling and simulation of mixture flows: Application to powder-snow avalanches. *Computers & Fluids*, 107:100-122, 2015.

[6] A. J. Chorin. The numerical solution of the navier-stokes equations for an incompressible fluid. *Bulletin of the American Mathematical Society*, 73(6):928-931, 1967.

[7] A. J. Chorin. Numerical solution of the navier-stokes equations. *Mathematics of computation*, 22(104):745-762, 1968.

[8] A. J. Chorin. On the convergence of discrete approximations to the navier-stokes equations. *Mathematics of computation*, 23(106):341-353, 1969.

[9] F. Clarelli, C. Di Russo, R. Natalini, and M. Ribot. Mathematical models for biofilms on the surface of monuments. In *Applied and industrial mathematics in Italy III*, pages 220-231. World Scientific, 2010.

[10] F. Clarelli, C. Di Russo, R. Natalini, and M. Ribot. A fluid dynamics model of the growth of phototrophic biofilms. *Journal of mathematical biology*, 66(7):1387-1408, 2013.

[11] N. Cogan and J. P. Keener. The role of the biofilm matrix in structural development. *Mathematical Medicine and Biology*, 21(2):147-166, 2004.

[12] F. Coquel, T. Gallouët, J.-M. Hérard, and N. Seguin. Closure laws for a two-fluid two-pressure model. *C. R. Math. Acad. Sci. Paris*, 334(10):927-932, 2002.

[13] D. Dutykh, C. Acary-Robert, and D. Bresch. Mathematical modeling of powder-snow avalanche flows. *Studies in Applied Mathematics*, 127(1):38-66, 2011.

[14] T. Gallouët, P. Helluy, J.-M. Hérard, and J. Nussbaum. Hyperbolic relaxation models for granular flows. *M2AN Math. Model. Numer. Anal.*, 44(2):371-400, 2010.

[15] T. Goudon and A. Vasseur. On a model for mixture flows: Derivation, dissipation and stability properties. *Archive for Rational Mechanics and Analysis*, 220:1-35, 2016.

[16] A. V. Kazhikhov and S. Smagulov. The correctness of boundary-value problems in a diffusion model of an inhomogeneous liquid. In *Soviet Physics Doklady*, volume 22, page 249, 1977.

[17] S. Labarthe, B. Polizzi, T. Phan, T. Goudon, M. Ribot, and B. Laroche. A mathematical model to investigate the key drivers of the biogeography of the colon microbiota. *Journal of Theoretical Biology*, 462:552-581, 2019.

[18] B. W. Peterson, Y. He, Y. Ren, A. Zerdoum, M. R. Libera, P. K. Sharma, A.-J. van Winkelhoff, D. Neut, P. Stoodley, H. C. van der Mei, and H. J. Busscher. Viscoelasticity of biofilms and their recalcitrance to mechanical and chemical challenges. *FEMS microbiology reviews*, 39(2):234-245, 03 2015.

[19] B. Polizzi, O. Bernard, and M. Ribot. A time-space model for the growth of microalgae biofilms for biofuel production. *Journal of Theoretical Biology*, 432:55-79, 2017.

[20] B. Polizzi, A. Fanesi, F. Lopes, M. Ribot, and O. Bernard. Understanding photosynthetic biofilm productivity and structure through 2d simulation. *PLOS Computational Biology*, 18(4):1-20, 04 2022.

[21] L. Preziosi and A. Tosin. Multiphase modelling of tumour growth and extracellular matrix interaction: mathematical tools and applications. *Journal of Mathematical Biology*, 58(4):625, Oct 2008.

[22] K. Rajagopal and L. Tao. *Mechanics of mixtures*. Series on Advances in Mathematics for Applied Sciences. World Scientific, Singapore, 1995.

[23] P. J. Schnurr and D. G. Allen. Factors affecting algae biofilm growth and lipid production: A review. *Renewable and Sustainable Energy Reviews*, 52:418-429, 2015.

[24] P. Secchi. On the motion of viscous fluids in the presence of diffusion. *SIAM Journal on Mathematical Analysis*, 19(1):22-31, 1988.

[25] R. Temam. Sur l'approximation de la solution des équations de navier-stokes par la méthode des pas fractionnaires (i). *Archive for Rational Mechanics and Analysis*, 32:135-153, 1969.

[26] C. Truesdell. Sulle basi della termomeccanica. I, II. *Atti Accad. Naz. Lincei. Rend. Cl. Sci. Fis. Mat. Nat.* (8), 1957.

[27] C. Truesdell. *Rational thermodynamics*. McGraw-Hill Book Co., New York-London-Sydney, 1969. A course of lectures on selected topics, With an appendix on the symmetry of the heat-conduction tensor by C. C. Wang.

[28] C. Truesdell and R. Toupin. The classical field theories. In *Handbuch der Physik, Bd. III/1*. Springer, Berlin, 1960. With an appendix on tensor fields by J. L. Ericksen.
# Appendices

## A. STATIONARY STATES

The equilibrium states for the system (12) correspond to the state solution where the source terms of all phases vanish, namely: $\Gamma_i = 0$ for $i \in \{\mathcal{A}, \mathcal{E}, \mathcal{L}\}$. In particular, $\Gamma_{\mathcal{E}} = 0$ induces $\psi_e + \eta_{\mathcal{E}}\psi_d = 0$ which lead to $\phi_{\mathcal{A}} = 0$. Thus, the spatiotemporal stationary state are the null state for microalgae, namely $\phi_{\mathcal{A}} = 0$ and $\mathcal{E} + \mathcal{L} = 1$.

## B. TIME DYNAMIC OF THE VOLUME FRACTION SUM

<figure><img src="image_8.png" /><figcaption>FIGURE 7. Time evolution for the maximal error within the domain on the sum of volume fractions: $E = \max_x |\Sigma_\alpha \phi_\alpha - 1|$.</figcaption></figure>

## C. NUMERICAL CONVERGENCE ANALYSIS FOR THE MODEL INCLUDING VISCOSITY

Numerical experiments have shown that numerical parameters need to be significantly reduced to reach convergence when the viscosity is included, and re-estimating the elastic tensor accordingly. Indeed, as mentioned in section 4.1, when including the viscosity, the elastic tensor coefficients must be rescaled and multiplied by $9 \cdot 10^7$ to obtain realistic front velocities for the biofilm. Figure 8 shows the convergences of the numerical scheme in both cases: with and without the viscous term. As expected, numerical convergence is obtained in both cases. Nevertheless, as presented in Figure 8 the convergences rate is lower when the viscous term is included. This explains at least partially why the numerical parameters need to be significantly reduced to reach acceptable precision for the application considered when the viscosity is included.

(A) Numerical convergence for the model without the viscous term.

<figure><img src="image_9.png" /></figure>

(B) Numerical convergence for the model including the viscous term.

<figure><img src="image_10.png" /><figcaption>FIGURE 8. Numerical analysis of the convergence of the scheme. The left side corresponds to the case without the viscous term, namely $\mu_\alpha = 0$ for $\alpha \in \{\mathcal{A}, \mathcal{E}, \mathcal{L}\}$ and the right side corresponds to the case with the viscous term and using very large values for the elastic tensor coefficients.</figcaption></figure>
D. RELATIVE VARIATION OF SOLUTES CONCENTRATION

Figure 9 shows the relative variation of solutes concentration for different times associated to the simulation presented in subsection 4.3. In this figure, we observe that the variations relatively to the input concentration are larger at the beginning (ie. *t* = 120h) than at the end (ie. *t* = 360h).

<figure><img src="image_11.png" /><figcaption>FIGURE 9. Solutes relative concentration for different times</figcaption></figure>