A network model for the biofilm growth in porous media and its effects on permeability and porosity

Successful microbial enhanced oil recovery depends on several factors like reservoir characteristics and microbial activity. In this work, a pore network is used to study the hydrodynamic evolution over time as a result of the development of a bioﬁlm in the pores. A new microscopic model is proposed for bioﬁlm growth which takes into account that nutrients might not fully penetrate the bioﬁlm. An important novelty in this model is that acknowledges the continuous spreading of the bioﬁlm over the network. The results from the current study can be used to obtain a new relation between the porosity and permeability which might be used as an alternative to the Kozeny Carman relation.


Introduction
The production of oil from the reservoir is initially accomplished by the internal pressure of the reservoir.However, when the primary production declines some external forces have to be applied, hence waterflooding or gas injection techniques are implemented to extract oil from the reservoir.These injection schemes are called the secondary oil recovery production.Nevertheless, even after primary and secondary recovery two-thirds of the oil are still trapped in the ground [4].The tertiary oil recovery extraction aims to increase the mobility of the remaining oil.One of the tertiary (or enhanced) oil recovery techniques is the microbial enhanced oil recovery (MEOR) which uses the growth of bacteria and the resulting by-products in order to increase the oil production.Microbial growth may enhance oil displacement by increasing the efficiency of the waterflooding process, by reducing interfacial tension and by changing the rock wettability [1,14].It has been observed that interfacial tension reduction and the increase of waterflooding efficiency Communicated by Gabriel Wittum.
B Fred Vermolen F.J.Vermolen@tudelft.nl 1 Delft Institute of Applied Mathematics, Delft University of Technology, Van Mourik Broekmanweg 6, 2628 XE Delft, The Netherlands caused by selective plugging are the mechanisms that have the largest impact on oil recovery [24].
Since it is hard to quantify the relation between the successful application of MEOR and parameters like the individual reservoir characteristics and the microbial activity, the development of computational models is of vital importance.These models are used to predict the bacterial growth and the in-situ regeneration of bioproducts in order to develop a proper field strategy [24].The influence of biofilm growth on porous media characteristics such as permeability and porosity has been modeled in several studies [5,6,9,12,16,18,29].The mathematical description is based on a theoretical framework and phenomenological relations obtained from experimental results [5,9,16,18,22,29].Biofilm growth models include Darcy continuum models [27,33], bacterially-based models [16], Lattice Boltzmann based simulations [17,32] and Pore Network Models (PNM) [5,9,18,21,26,31].Usually, in biofilm growth models the porous medium consists of three components: the grains, the biofilm which grows on the walls of the solid grains and the liquid in the pore space.The grains are assumed to be impermeable to the liquid and the nutrients, therefore hydrodynamic model equations are written only for the liquid and biofilm [17].
Cunningham et al. [8] showed experimentally the effect of the accumulation of biofilm on the porosity, permeability and friction factor of the porous media.The porosity of the media decreased between 50 and 96% due to the accumula-123 tion of biofilm, while permeability decreased between 92 and 98%.Taylor and Jaff [28] obtained an analytic expression to describe changes in the porous media as a result of biofilm growth in the continuum scale.However, in Taylor and Jaff [28] it is assumed that biofilm growth proceeds uniformly through the network which is an oversimplification according to laboratory experiments [15].Clement et al. [6] model the biofilm growth using a macroscopic approach.This model does not assume any specific pattern for biofilm accumulation, instead it is based on macroscopic estimates of average biomass concentrations.Seki and Miyazaki [23] proposed a mathematical model for bioclogging that takes into account the nonuniform microbial distribution of colonies which ranges from micro-colonies to biofilm.However assuming uniform biofilm thickness in their model gives an overestimation of the bioclogging process [36].Therefore, pore network models and pore-scale models are needed to describe the growth of biomass and its effects on the macroscopic properties of the media properly [35].
In PNMs, the porous medium is modeled by cylindrical interconnected tubes in which water or any fluid can flow.The biofilm development is stimulated by the injection of nutrients into the network.Transport of nutrients takes place within an aqueous phase and is described by a convectiondiffusion-reaction equation.The reaction term models the consumption of nutrients caused by bacterial population growth.The bacterial population will determine the development of biofilm in the pores of the medium.This biofilm will grow and will change the radii of the pores, leading to porosity and permeability reduction and hence to a modification in the flow pattern dynamics of the fluid that carries the nutrients through the network [5,26,31].Kim and Fogler [12] studied the effects of biofilm growth on porosity under starvation conditions.They show a good agreement with experimental results and show the existence of a critical shear stress.Raoof and Hassanizadeh [19] used a pore network model to describe two-phase flow in a porous media.They took into account the influence of the nodes of the network on the effective resistance of the fluids.They used a coordination number distribution which allows a maximum coordination number of 26.Additionally, they assigned a variety of cross-sectional shapes including circular, rectangular and triangular.They claimed that the inclusion of the volume of the nodes of the network affects the relation between the relative permeability and the saturation of the fluids.Despite the relevance of their work, in their model, they did not include the development of biofilm in the porous medium.In this study,as an approximation, we disregard the volume of the nodes to avoid additional complications in the model.Arns et al. [2] studied the effect of topology in the relative permeability of the networks.They found that the relative permeability curves obtained with stochastic networks are in good agreement with the ones obtained from imaged rock networks.The bacte-ria and Extracellular Polymeric Substance (EPS) in porous media are often lumped together and are represented as a continuous uniform layer of biomass attached on the surface of the solid grains of the porous media [9,26,31].This uniform layer of biomass is referred to as biofilm.Furthermore, the biofilm growth rate is usually assumed to be proportional to the volume of biomass.Nevertheless, the nutrients might not be available in the entire volume of the biofilm.This phenomenon occurs if the consumption of nutrients is faster than the diffusion rate within the biofilm so that the (diffusion) penetration of the nutrients into the biofilm proceeds at a slower rate than the other processes [10,13,25].Hence, the hypothesis that the nutrients are distributed over the whole volume of biofilm is questionable.Therefore, we assume that biofilm growth occurs only in a limited volume where the concentration of nutrients is maximal.
Usually, in PNMs the microbial activity is assumed to exist only within the tubes and no spread of biomass between neighboring tubes is described [5,18,21,26,31].However, experiments show that the biomass or biofilm continuously grows, extending through the whole medium [15].To model the inter-pore transport, Ezeuko et al. [9] consider a spreading potential among neighboring tubes.The spreading of the biofilm is allowed once the biomass has completely saturated the host pore.Thullner et al. [31] modeled the colony growth by assuming that a tube in the network was completely full or empty.Hence a binary switch mechanism is used to describe the spreading of biomass.The switch to completely filled tubes is determined by the size of the tubes.However, they did not consider any exchange of biomass between neighboring tubes.In our model, we describe the continuous spreading of the biofilm between adjacent tubes by computing the spreading of biomass from one pore to its neighbors, if there is a difference of volume of biomass between neighboring tubes.
In this study, we present a new biofilm growth model which takes into account that nutrients cannot fully penetrate the biofilm since consumption of nutrients is faster than the diffusion rate through the biofilm.We take into account that the biofilm growth is limited within a thin penetration layer, in which bacteria are in direct contact with the nutrients.In our model, there are two types of biofilm development: growth in the interior of the tube and growth at the extremes of the tube.Biofilm growth in the extremes of the tube will lead to the spreading of the biofilm to the neighboring tubes and through the whole network.The currently proposed biofilm growth model approach has several advantages over other models.Firstly, we incorporate the likely non-homogeneous distribution of the nutrients within the biofilm.Secondly, since biofilm growth takes place mainly in the boundary between water and biofilm, the internal biofilm growth will naturally stop if the tube is full of biofilm.Finally, the biofilm growth in the extremes of the tubes leads to spreading of biomass through the whole network.In this model there is no need to seed initially all the tubes in the network to observe the clogging of the network.This paper is focused on the presentation of biofilm growth model in a pore network.Future research might be the used of these results to obtain an alternate relation between porosity and permeability.The up-scaling of these results is beyond the scope of this paper.

Mathematical model
We represent the porous medium as a 2D rectangular network composed of interconnected cylindrical tubes.The point where these tubes are connected is called a node of the network and is indexed as node n i .The tube between the node n i and n j is indexed as the tube t i j (see Fig. 1).We assume that the radius is the same for all the tubes (which differs from previous studies because we want to express the spreading of the biofilm in a simple way, the modeling of this phenomenon is explained later) and the same length l.The number of tubes connected in each node is four for interior nodes, three for boundary nodes and two for the nodes in the corners of the network.We assume the bacteria and the biofilm are lumped together and hence we refer to them as the single phase: biofilm.We assume that nutrients are injected through the network and transported within a fluid phase.For simplicity we chose water as the fluid in which the nutrients are transported.We define the thickness of the biofilm in the tube t i j by r b i j , the radius available for water by r w i j and the total radius of the tube by R (see Fig. 1).The volumetric flow of the water phase q i j in the tube t i j is described by a modified form of the Poiseuille equation [30], where Δp is the pressure drop between neighboring nodes, μ is the viscosity of water that flows in the bulk, l is the length of the tube and the dimensionless number β is the ratio between the viscosity of water flowing through the biofilm and the viscosity of water flowing through the bulk.We use β = 10 7 which according to Thullner and Baveye [30] is a good approximation for an impermeable biofilm.Mass conservation is imposed on each of the nodes.For the node n i we have where, and further q i j is the flux in the tubes connected to node n i .The balance of nutrients is described by an advectiondiffusion-reaction equation.Denoting the concentration of nutrients by C, this gives where u is the advection velocity related to the local flux q by u = q/A.Here A denotes the area of the cross-section of the tube and D is the diffusion coefficient of water.Further b + represents the concentration of biofilm that grows as a result of consumption of nutrients (no detachment term is taken into account in this equation).In general, the concentration of nutrients b is linked to the volume of biofilm, where ρ and V T , respectively, denote the mass density of biofilm and the total volume of the tube.We describe the overall growth rate of the biofilm in the following paragraphs.
In this model we assume that nutrients might not penetrate completely through the biofilm since the reaction is 123 faster than the diffusion rate within the biofilm.Therefore we propose that there exists a maximal distance (from the water biofilm interface) that the nutrients can travel within the biofilm.The maximal distance is called the penetration layer, Γ , and implicitly defines a maximal volume in which the nutrients can diffuse.This volume is called the penetration volume V p and it is assumed to be constant during the whole process of biofilm growth.If the volume of biofilm is smaller than the penetration volume, the nutrients can penetrate the whole biofilm volume and hence the biofilm growth rate is proportional to the volume of biofilm.However, if the biofilm volume is much larger than the penetration volume, the nutrients react with the biofilm only within this penetration volume, which is adjacent to the water-biofilm interface.In this case, the biofilm growth rate is proportional to the area of the water biofilm interface.Further, since in general there are two regions in the tube where the biofilm encounters the nutrients, we model two kinds of biofilm growth: internal biofilm growth and biofilm growth at the extremes of the tube (see Fig. 2).Firstly, we describe the internal growth.
The biofilm growth rate in the interior of the tube t i j is modeled as follows, Here V i b f i j denotes the volume of interior biofilm in the tube t i j .Further, f (V b f i j ) ≥ 0 is the positive part of a sigmoid function for V b f i j that depends on the penetration volume V p , k 1 is the specific biofilm growth rate, A i wb f is the internal interface water biofilm area, A i T is the external area of the tube, C i j is the concentration of nutrients within the tube and E s is a saturation constant.The positive part of the sigmoidlike function is defined as, The dependence of biofilm growth rate (Eq.6) on the concentration of nutrients is given by the Monod equation which models the limiting nutrient consumption by the biofilm.
Next we explain the reason why we use the positive part of a sigmoid-like function: if the volume of biofilm is small, i.e.
Therefore, the biofilm growth rate is proportional to the volume of biofilm, . If the volume of biofilm is much larger than the penetration volume, V b f V p , then f ∼ 1 and therefeore the biofilm growth is proportional to the area of the interface between water and biofilm, Water-biofilm interface in the interior of the tube Water-biofilm interface in the extremes of neighboring tubes The biofilm growth rate is zero when there is no biofilm in the tube or when the tube is filled with biofilm, consequently, biofilm growth in the interior of the tube stops if there is no more space in the tube.Note that our approach is phenomenological.Further, the area A i wb f can be written in terms of the total volume of the tube V T and the volume of biofilm V b f , therefore Eq. ( 6) for the biofilm which grows in the interior of the tubes, is expressed as follows, Note that the above relation for V i b f i j represents a continuous relation of biomass growth with the volume of biofilm.
Secondly, we describe the biofilm that grows in the extremes of the tube.Since the penetration volume in the extremes is very small compared to the whole volume of biofilm, the biofilm growth rate in the extremes of the tubes is proportional to the area of the interface between water and biofilm (see Fig. 2).We assume binary interactions with the neighboring tubes.The area of the interface between water and biofilm A e wb f between the tube t i j and the tube t jk can be written in terms of the difference between biofilm volumes of neighboring tubes.The biofilm grows in the extreme of the tube with a larger volume of biofilm and it is given to the neighboring tube which has a smaller volume of biofilm.
If we assume that the volume of biofilm V b f jk in the neighboring tube t jk (connected to the node n j ) is larger than the volume of biofilm V b f i j in the tube t i j , then the biofilm growth in the extreme of the neighboring tube t jk is given by, Here V e b f jk represents the volume of biofilm at the extreme of the tube, A e wb f is the external interfacial water biofilm area and A e T is the cross-sectional area in the extreme of the tube.The ratio between the external interfacial water biofilm area and the cross-sectional area of the tube A e T is a measure of the biofilm growth in the extremes of the tube.This ratio is zero if the volume of biofilm is the same in both 123 interacting tubes which means there is no biofilm growth in the extreme of the tube and hence no volume of biofilm is added to either of them.On the other hand, when this ratio is one, the biofilm grows at a maximal rate and the accumulated biofilm is added to the tube t i j .Note that there is no biomass exchange between neighboring tubes; the biomass is produced in the extreme of the tube and it is given to the neighboring one, hence no loss term for the biomass growth is necessary to describe this phenomenon.In this way, this model for the biofilm growth allows the spreading of the biofilm through the whole network, which is consistent with experimental observations.The area A e wb f between the tube t i j and the tube t jk can be written in terms of the volume of the biofilm of the tubes.Hence Eq. ( 9) for the biofilm growth at the extreme of the tube t jk changes into, We take into account all the neighboring tubes whose volumes of biofilm are larger than the volume of biofilm in the tube t i j .To this extent we introduce the following index set notation for the tube t i j which connects nodes n i and n j .Consider the node n j then we define the set of neighboring nodes of it, except n i by Λ ji (see Fig. 1).Therefore, taking into account all neighboring tubes, the equation for the biofilm growth in the tube t i j due to biofilm growth in the extremes of the neighboring tubes is written as, in which we use the notation we combine the internal growth of biofilm with the biofilm growth in the extremes of the neighboring tubes and including possible detachment of biofilm, which is proportional to the interfacial water-biofilm area, we obtain, where k 2 is the detachment rate coefficient.Further, H (V b f i j ) is defined as, We include the function H because detachment occurs only when there is biofilm within the tube.In case there is no biofilm in the tube, H = 0, which means the detachment rate is zero.In Eq. ( 12) the first term is the interior biofilm growth, the second and third term describes the biofilm which grows in the extremes of the neighboring tubes and the fourth term is a term for the detachment of the biofilm.
The reaction rate of the nutrients is given by, In summary, we solve the following coupled mathematical problem: Find p, subject to such that, where, Here L x is the size of the network in the x direction and L y the size in y direction.Next to this r w i j decreases as a result of deposition of biofilm, which grows under the presence of nutrients.The balance of nutrients is given by, subject to, C(x, y, t 0 ) = 0, The biofilm grows according to, Subject to the initial condition Our routine randomly chooses 4% of the tubes.Note that b i j = ρV bf i j V T .We have chosen b 0 = 10 −4 (kg/m 3 ) for the initial tubes that were seeded with biofilm.The consumption of nutrients is modelled by, links the concentration of nutrients in the nodes and the concentration in the tubes.

Numerical method
The numerical approach and the computational procedure used in this work are described in this section.When mass  1), is combined to Eq. (2) a linear system for the pressures at the nodes p i arises.After solving this system, the flux q i j in each of the tubes is computed.
The equation for the balance of nutrients is solved for the concentration C i at each node n i of the network (see Fig. 3).To discretize the equation for the balance of nutrients, we write it in the following form, where the first-order upwind scheme for the advection term gives, where Further, the diffusion term of the Eq. ( 4) is discretized using a time-implicit method for the concentration.The area is used from the previous time step.We use a finite difference scheme in space.Therefore the discretization for the diffusion part, reads as, where A τ w i j is the area of the cross section of the bulk water in the tube t i j and A i j is the total area of cross section of the tube t i j .
To write the reaction term in each node, we assume that at each node there is a perfect mixture of biofilm.Therewith we get, where, The solution of the concentration of nutrients obtained from the advection-diffusion-reaction is used for the approximation of the biofilm volume.
The biofilm growth takes place within the tubes of the network.Here we use an explicit Euler time integration method to arrive at, The computational procedure used in this work is as follows.Firstly, the pressure is imposed in the left and right boundary of the network.Subsequently, the pressure in each node is computed from the linear system resulting from the mass conservation in each node.For solving this system, we consider Dirichlet boundary conditions in the left and right boundaries and homogeneous Neumann boundary condition for the upper and lower boundary.The pressures in each node are used to compute the flux in each tube by means of Eq. ( 1).After this step, we proceed to solve the advection diffusion reaction equation for the nutrients and we compute the concentration of nutrients in each node as well as the volume of biofilm in the tubes.The thickness of the biofilm and the radius of the void space available for water is updated and the process starts again at the next time step (See Fig. 4).

Concentration of biomass
Fig. 4 Full model algorithm that combines the transportation of nutrients and biofilm growth

Simulation results
In this section we describe the numerical experiments and the results obtained for the biofilm growth in a pore network.Firstly, in order to validate the advection-diffusion part of our code, we compare our results with an analytic solution and with a Continuous Time Random Walk (CTRW) transport model [7].Secondly, we studied the biofilm growth effects on the outflux and porosity.For this study, the biofilm growth rate k 1 is fixed but three different detachment rates k 2 are used.Finally, we compare our results with the Kozeny-Carman relation and with two quasi-steady biofilm growth models.
Firstly, the evolution of the concentration of nutrients through the network is studied without the presence of biofilm.We disregard the reaction term in order to be able to compare the transport and diffusion of nutrients with an analytic solution in 1-D and with an existing model based on CTRW.The CTRW transport model can consider classical and non-classical Fickian dispersion.In this case we use Fickian diffusion for the CTRW since in our model we are not considering other kinds of diffusion.We use a MATLAB toolbox developed by Cortis and Berkowitz [7] to obtain the breakthrough curve with the CTRW model.The diffusion coefficient and the pore velocity used in the CTRW transport model are listed in Table 1.
We solve the advection-diffusion equation for the concentration of nutrients with our model using a mesh with 201×11 123 nodes, which means there are 201 nodes in x direction and 11 nodes in y direction.The number of tubes is determined implicitly by the number of nodes and by the topology of the network.Further, we assume that all the tubes in the network have the same radius.We use the volumetric flows through the pores from the network model for the solution of the concentration of nutrients.Under these conditions for the size of the mesh and the uniform size of the radii in all the tubes, we can compare the results with a model based on CTRW and with an analytic solution in one dimension [34].The analytic solution of the advection-diffusion equation (Eq. 4 without reaction term) in 1-D is given by: in which erfc is the complementary error function, v * is the velocity and D * is the diffusion coefficient used in this first simulation.
Figure 5 shows the results for the normalized concentration of nutrients C/C 0 in one of the tubes that is located adjacent to to the outlet of the network for our model, a model based on CTRW and the analytic solution given by Eq. ( 30).We observe a good agreement between, the CTRW model, the analytic solution and our model, which indicates that our scheme produces consistent results.However we observe a small shift between our model, the CTRW and the analytic solution.The shift is attributed to the following cause: our model contains a Neumann boundary condition at the outflow boundary, whereas the analytic solution is valid in a domain with infinite size.Therefore the concentration calculated by our model is a little higher than the one computed the use of the analytic solution.This can be proved in more rigour using smoothness of the solution and the maximum principle.The Yield coefficient Y 0.34 [3] Half saturation constant for biofilm Biofilm/bulk water viscosity ratio β 10 7 [30] complete set of parameters for this simulation is presented in Table 1.The next step is to quantify the effects of biofilm growth on the porosity and permeability of the porous medium.Therefore, we solve the biofilm growth and the transport of nutrients as a coupled problem.Initially 4% of the tubes are seeded with an initial concentration of biomass b 0 = 1 × 10 −4 (kg/m 3 ).We performed three sets of simulations in which the biofilm growth rate is fixed, however three different values for the detachment rate factor are chosen, k 2 = 10 −6 (1/s), k 2 = 10 −7 (1/s) and k 2 = 0 (1/s).For this set of simulations we used a network with 101 x 61 nodes and we considered a radius R = 1.1937 × 10 −5 (m) for all the tubes of the network.The complete set of parameters for this set of simulations is listed in Table 2.For each pair of biofilm growth k 1 and detachment rate factor k 2 , we performed ten simulations where we fixed all the parameters, except the initial distribution of tubes seeded with biofilm.The normalized flux Q n is defined as, Q n = Q Q 0 , where Q 0 is the initial flux in the network (i.e.before biofilm growth).We compute the average of the normalized flux and we observe that the 95% confidence interval is very close to the average value of the normalized flux, therefore the initial random biofilm distribution does not have a significant effect on the results.
The evolution of the average normalized flux through the network for the detachment rate k 2 = 10 −7 (1/s) and k 2 = 0 (1/s) is shown in Fig. 6.For detachment rates k 2 = 10 −7 (1/s) and k 2 = 0 (1/s), we observe a decrease of the normalized flux due to the accumulation of biomass in the network.However, for k 2 = 10 −6 (1/s) the detachment of biofilm dominates over biofilm growth and the initial distribution of biomass is removed during the first stage.Therefore, in this case no biofilm develops in the medium and there will be no changes in the permeability and porosity of the network.This implies that Q n = 1 at all times.If the biofilm detachment rate is smaller, the development of biofilm attached to the walls of the pores leads to a reduction in the radius available for the water flow and consequently biofilm growth leads to a reduction of the normalized flux of the network.We observe very similar behavior for the detachment rate k 2 = 10 −7 (1/s) and k 2 = 0 (1/s).
In Fig. 7 the average of the fraction of biofilm volume is presented for k 2 = 10 −7 (1/s) and k 2 = 0 (1/s).The fraction of volume of biofilm in the network is given by . The sum is taken over all the tubes in the network.
Since we neglect the volume of the nodes, the volume of the tubes corresponds to the volume of the pore space.We observe that during the first minutes the volume of biofilm in the network increases monotonically for the two cases.Further, after approximately 300 min the biomass growth reaches a steady state.We observe that approximately 32% of the void space of the network is occupied by volume of biomass at the steady state for both cases.
Finally, in addition to the full model which considers the transport of nutrients and the biofilm growth as two coupled phenomena, two quasi-steady state models of biofilm growth are also considered in this work.In these models we set an amount of volume of biofilm in the network, then we compute the effect of the volume of biofilm in the radius available for water and finally we compute the flux through the network.Note that the transport-diffusion equation is not solved in these models.
In the first model we consider that initially biofilm is present in all the tubes of the network and that the biofilm grows at the same rate in all the tubes.Therefore we refer to this model as uniform biofilm growth.
In the second model we hypothesize that each tube in the network could either be completely filled with biofilm or completely empty.We vary the percentage of tubes filled with biofilm from 1% of the tubes to 100% of the tubes.In each stage, the tubes filled with biofilm are chosen randomly.We refer to this model as random biofilm growth.We perform 10 simulations and we determine the average flux in the outlet of the network.We found that the variance of the result was very small.We compare the results of the full biofilm growth model with the uniform growth, with the random growth and with the Kozeny-Carman relation.The Kozeny-Carman is a well known equation that provides a relation between the porosity φ and the permeability K and it is given by the following equation, 123 in which C k is a parameter related to the specific internal surface area of the pores in a porous media.In order to be able to compare the full model with the uniform growth model, with the random growth and with the Kozeny-Carman equation, we have to express the volume of biofilm in terms of porosity and the normalized flux in terms of the permeability.The relation between the fraction of biomass and porosity is given by the following equation, in which φ 0 is the initial porosity.
The relation between the normalized flux Q n and the permeability is determined by the Darcy's Law, If the pressure drop ΔP, the cross-sectional area of the network A n , the length in the flux direction L and the viscosity μ are constant during the process of biofilm growth, we have that In which K 0 is the initial permeability.Then, using Eqs.(31) and (34) we can write the normalized flux predicted by Kozeny-Carman as follows, Note that in order to derive Eq. ( 35) the parameter C k has been taken constant.However, since the porous medium channels are changed by the non-uniform accumulation of biomass, the assumption of taking C k constant is probably inappropriate.Hence, our results may deviate from the results predicted by the Kozeny-Carman model.
In Fig. 8 the numerical results for the porosity φ versus the normalized flux are shown for the two different detachment rates studied in this work k 2 = 10 −7 (1/s) and k 2 = 0 (1/s), the two cases of quasi-steady-state biofilm growth models and the Kozeny-Carman relation [Eq.(35)].
The uniform growth model and the full model overlap from the initial porosity to 0.35 approximately where a sudden decrease in the normalized flux is described in the full-model for k 2 = 0 (1/s) and k 2 = 10 −7 (1/s).This is explained as follows: in the beginning in the full model, the biomass starts spreading through the network and since the thickness of the biomass in the tubes is still small, the influence of the biomass on the permeability is insignificant at this stage.However, since the nutrients are transported through the network and the biomass is spread continuously, uniform biofilm growth is stimulated in the network, causing a decrease in the permeability due to the accumulation of biomass.Afterwards, the nutrients are consumed by the bacteria and the biofilm starts growing and clogging the pores, therefore there is a reduction of the flux of the nutrients in the whole network.Hence, the nutrients are present preferentially near the inlet which causes a preferential growth of biofilm near the inlet and at the final stage causes the total decrease of the flux.The random growth model shows a linear decay of the normalized flux.For high porosity the slope of the decay of the normalized flux predicted by the random growth model is similar to the slope of the normalized flux predicted by the full model.The linear behavior of the random growth model deviates from the full model for lower porosity.The random biofilm growth predicts a plugging of the network when the porosity is about 0.2.The porosity is approximately half the initial porosity, which is in accordance with the percolation threshold for a rectangular network, [11].The fact that the full model stays in accordance with the uniform growth model seems to indicate that at the beginning the time evolution of the flux is predominantly determined by the localized growth kinetics of the biofilm, rather than the kinetics of spreading over the network.Finally, the Kozeny-Carman relation shows a behavior that is similar to the uniform biofilm growth, but the decrease of the normalized flux with the decrease in porosity is faster than the uniform growth.

Conclusions and outlook
In this work, we simulate biofilm growth, in particular its effects on the porous medium characteristics such as porosity and permeability.We use a two-dimensional pore network model to represent the porous medium.We develop a new model for biofilm growth, which predicts that the nutrients are not able to penetrate fully in the biofilm if the reaction term is dominant over the diffusion of nutrients within the biofilm.In addition, our model is able to simulate the spreading of the biofilm through the whole network which is a phenomenon that has been observed experimentally [15].The proposed model shows that at early stages biofilm growth is mostly uniform through the whole network, however eventually the biofilm will grow preferentially near the inlet of the network, plugging the pores at the inlet and causing a cease of the flux through the network.The modifications in porosity and permeability caused by biofilm growth might be beneficial for a microbial enhanced oil recovery technique, especially in the first stage before the plugging of the network.Since we see that uniform growth provides a relatively good correspondence with the full model for high porosity, we conclude that the clogging of the porous medium in high permeability layers is feasible without blocking the inlet.For this reason, we propose to stop injection of nutrients in order to avoid plugging the medium.This behavior is not described by the uniform growth model, the random growth nor the Kozeny-Carman relation.
Since we consider a 2D rectangular pore network model consisting of cylindrical tubes with the same radii, this model could be too simplified to describe a real reservoir field.Interesting further research is to find the representative elementary volume in order to upscale these results to the macroscale.In addition, future plans entail the study of the effects of biofilm growth in porosity and permeability in more complex topologies in 2D and 3D.

Fig. 1
Fig. 1 Pore network and biofilm thickness within a tube

Fig. 2
Fig. 2 Biofilm growth in the interior and in the extremes of the tubes

Fig. 3
Fig. 3 Network discretization and domain of computation

Fig. 5 2 3 )
Fig. 5 Comparison of the solution of the advection difusion equation for our model, CTRW model and an analytic solution

Table 1
Parameters for the simulation without biofilm growth Comparison of the normalized flux vs porosity for the full model, the random growth model, the uniform growth model and Kozeny-Carman relation