Conditions for upscalability of bioclogging in pore network models

In this work, we model the biofilm growth at the microscale using a rectangular pore network model in 2D and a cubic network in 3D. For the 2D network, we study the effects of bioclogging on porosity and permeability when we change parameters like the number of nodes in the network, the network size, and the concentration of nutrients at the inlet. We use a 3D cubic network to study the influence of the number of nodes in the z direction on the biofilm growth and on upscalability. We show that the biofilm can grow uniformly or heterogeneously through the network. Using these results, we determine the conditions for upscalability of bioclogging for rectangular and cubic networks. If there is uniform biofilm growth, there is a unique relation between permeability and porosity, K ∼ ϕ2, this relation does not depend on the volume of the network, therefore the system is upscalable. However, if there is preferential biofilm growth, the porosity-permeability relation is not uniquely defined, hence upscalability is not possible. The Damköhler number is used to determine when upscalability is possible. If the Damköhler number is less than 101, the biofilm grows uniformly and therefore the system is upscalable. However, if the Damköhler number is greater than 103, the biofilm growth exhibits a deviation from uniform biofilm growth and heterogeneous growth is observed, therefore upscalability is not possible. There is a transition from uniform growth to preferential growth if the Damköhler number is between 101 and 103.


Introduction
In primary oil production, a wellbore is drilled from the surface to the ground and oil is extracted from the reservoir by natural mechanisms such as the internal pressure of the reservoir. When the initial production declines secondary oil recovery techniques such as waterflooding or gas injection are implemented. However, two thirds of the oil are still trapped in the ground even after primary and secondary recovery [5]. Microbial enhanced oil recovery (MEOR) is a tertiary oil recovery technique which aims at increasing the mobility of the remaining oil using the growth of bacteria and the resulting by-products. Bacterial growth enhances oil recovery by increasing the efficiency of the waterflooding process, by clogging highly permeable layers such that the flow through the oil-containing regions, with low permeability, is enhanced. Furthermore, bacterial growth reduces interfacial tension and changes the rock wettability [1,12] which enhances oil mobility.
The development of computational models is of vital importance to design a proper field strategy for MEOR in oil reservoirs. These models describe bacterial growth and predict the changes in the characteristics of the porous media like the permeability and porosity [22]. Among bacterial growth models in porous media, there are the Darcy continuum models [24,29], bacterially based models [16], Lattice Boltzmann-based simulations [9,17], and pore network models (PNM) [6,8,18,21,23,26,28]. The secretion of extracellular polymeric substances (EPS) by bacterial population causes the formation and growth of biofilm on the walls of the porous media. In biofilm growth models, it is usually assumed that the porous medium consists of three components: the grains, the biofilm, and the fluid which contains the nutrients needed for the biofilm growth. The equations that describe biofilm growth are written only for the fluid and biofilm since the grains are assumed to be impermeable to the liquid and the nutrients [17].
Cunningham et al. [7] studied experimentally the effect of biofilm growth on the porosity, permeability, and friction factor of the porous medium. They reported a decrease in the porosity between 50% and 96% and a decrease in permeability between 92% and 98% due to the accumulation of biofilm. In the continuum scale, Taylor et al. [25] obtained an analytic expression that describes the relation between porosity and permeability. However, they assumed that the biofilm grows uniformly through the domain of computation which not always occurred according to laboratory experiments [15]. Therefore, microscale biofilm growth models such as pore network models (PNM) and pore-scale models are used to describe a non uniform biofilm growth [8,32]. It is needed to state the conditions such that the biofilm grows uniformly.
In PNMs, the porous medium is usually represented as a two or three-dimensional lattice of cylindrical interconnected tubes in which water or any fluid can flow [4]. The biofilm development is caused by the injection of nutrients into the network which are transported within an aqueous phase. The injection and consumption of nutrients are described by a convection-diffusion-reaction equation in which the reaction term models the consumption of nutrients caused by bacteria which results in biofilm growth. The biofilm grows and it adheres to the walls of the cylinders. Thereby, the biofilm changes the radii of the pores, which consequently leads to porosity and permeability reduction [6,23,28].
Even though the bacterial population and the EPS are two different phases, they are usually lumped together and are represented as a continuous uniform layer of biomass attached to the walls of the pores [8,23,28]. This uniform layer of biomass is referred to as biofilm.
Biofilm growth models at the microscale are needed to account for heterogeneities in the biofilm growth, which are typically ignored in large-scale models. In PNMs, the information obtained at microscale is averaged over the network to get a macroscopic description at the continuum scale [13]. The influence of network characteristics such as the coordination number on macroscopic transport phenomenon has been shown in [30]. They showed that the dispersivity decreases if the coordination number increases.
In general, the purpose of upscaling is to get an effective description on a macro level when there exists a good description on a small-scale level [10]. Hese et al. [10] use upscaling methods to obtain an effective one dimensional representation based on a system of two-dimensional partial differential equations. Their study is focused on the scaling behaviour of Monod-type reaction kinetics. They showed that the upscaled description of Monod kinetics leads to a concentration dependent transition between a reaction and a diffusion limited regime. Wu et al. [31] computed the upscaled grid-block permeability from fine-scale solutions of the flow equation. They studied the upscaling of single phase flows through media with a periodic small amount of heterogeneity. They claim that their results are also useful for the understanding of the upscaling of random media. The equivalent permeability is a constant permeability that represents an heterogeneous medium. However, it is impossible to obtain a one-to-one mapping as a complete mapping between the real heterogeneous medium and the homogeneous upscaled medium. Therefore the equivalence, that is the one-to-one mapping, is defined in a limited sense [19]. Battiato and Tartakovsky [3] studied the transport of a solute in a porous medium which is subjected to a nonlinear heterogeneous reaction. This solute precipitates on the solid matrix to form a crystalline solid. They investigated the sufficient conditions under which the macroscopic advection-dispersion-reaction equations provide an accurate description of the pore-scale processes. Despite their relevant findings, they did not consider any change on the morphology of the porous medium. In the present study, we investigate under which conditions we can upscale a small-scale heterogeneous medium to a largescale homogeneous medium, in which we can apply models for uniform growth.
In this work, we study the effects of biofilm growth on the porosity and permeability of the network. We use the model for biofilm growth described in [14], which models incomplete transmigration of nutrients through the biofilm as a result of high bacterial consumption rate and a low diffusion rate of the nutrients. In particular, we study the process of bioclogging which features inhibition of the flux through the network due to biofilm growth. We compute the amount of biomass per volume needed to block the network for different number of nodes in the network, different network sizes and different inlet concentrations of nutrients in the network. Furthermore, we describe the conditions for uniformity and upscalability of the pore network biofilm growth model.
As long as the medium, in this case the pore network model, is evolving in a spatially homogeneous manner, upscaling can be performed on the basis of computing an effective porosity and permeability. In the current paper, however, we are dealing with the injection of nutrients on the inlet boundary. The nutrients are being consumed by the bacteria in the porous medium and thereby, effectively, converted into biomass that clogs the network tubes. If the concentration of the nutrients at the inlet boundary is not sufficiently high, then, the nutrients will all be consumed and converted before they are able to reach the regions in the domain that are further away form the inlet. If this happens in the current network, then the network cannot be used for upscaling purposes in which one determines the effective permeability and porosity relation. The current paper will address this issue in terms of derivation of a relation between permeability and porosity. For upscalable conditions we will derive a tractable, functional relation between the effective porosity and permeability of the network, which can be used as an alternative relation to standard Kozeny-Carman relation. We will analyze the applicability of such a relation by varying the characteristics of the pore network. The current network provides a computational approach to classical upscaling that is carried in more mathematical rigor.
The paper is organized as follows. In Section 2, we describe the equations for the transport of nutrients and the model used for the biofilm growth. In Section 3, we describe the effects of biofilm growth on porosity and permeability when we vary the number of nodes in the 2D network, the 2D network size and the inlet concentration of nutrients. In addition we present the results of the biofilm growth in a 3D cubic network when we vary the number of nodes in the z direction. Finally, in Section 4, we present the discussion, draw the conclusions, and present the outlook.

Mathematical model
In this section, we present the equations that describe the transport of nutrients and the biofilm growth in the porous medium. We consider two type of networks. Firstly, the porous medium is represented as a 2D rectangular network composed of interconnected cylindrical tubes whose radii and length are the same (see Fig. 1). Secondly, a cubic network is used to study the influence of the number of nodes in z direction on the porosity-permeability relation (see Fig. 2). We assume that the bacteria and biofilm are lumped together and hence we refer to them as a single phase: biofilm. The growth of biofilm is initiated by the nutrients which are injected into the network and transported within a fluid phase. The thickness of the biofilm in the tube t ij is represented by r b ij , the radius available for water by r w ij , and the total radius of the tube by R (see Fig. 3a). The volumetric flow of the aqueous phase q ij in the tube t ij is described by a modified form of the Poiseuille equation [27], 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 [27] is a good approximation for an impermeable biofilm. Despite the β term is only a negligible addition in Eq. 1, we incorporate it to keep the description general.
In each of the nodes, mass conservation is required. Therefore, for the node n i , we have where S i = {j | n j is adjacent to the node n i } and where q ij is the flux through the tube that connects node n j to node n i .
The transport of nutrients is described by an advection diffusion reaction equation. The concentration of nutrients is denoted by C, where D is the diffusion coefficient of nutrients through water and u is the average advection velocity which is related to the local flux q by u = q/A, where A is the area of the cross-section of the tube. Additionally, b + is the concentration of the biofilm produced as a result of the consumption of nutrients (no detachment of biofilm is taken into account in this term). The concentration of biofilm b is related to the volume of biofilm by, where ρ denotes the mass density of biofilm, V T the total volume of the tube and V bf the volume of biofilm. The reaction rate is described in the following paragraphs. In this work, we use the model for the biofilm growth reported in [14]. Since the reaction rate of nutrients is higher than the diffusion rate within the biofilm, it is assumed that nutrients interact with the biofilm only in a thin layer adjacent to the water biofilm interface, p . This layer defines implicitly a volume which is called the penetration volume of the nutrients V p ≈ 2πR p l, and it is assumed to be constant during the whole process of biofilm growth. In general, in each of the tubes, there are two different water biofilm interfaces. Therefore, we consider two modes of biofilm growth: internal biofilm growth and biofilm growth at the extremes of the tube. The interior biofilm growth takes place within the tube and is described as follows. If the volume of biofilm is smaller than the penetration volume V p , the nutrients are present in the whole biofilm volume and hence the biofilm growth rate is proportional to the volume of biofilm (see Fig. 3b). However, if the biofilm volume is much larger than the penetration volume, the nutrients are consumed only within this volume and the biofilm growth rate is proportional to the area between water and biofilm interface (see Fig. 3c).
The biofilm interior growth rate in the tube t ij , V i bf ij can be written as, In this equation, f (V bf ij ) ≥ 0 is a sigmoid-like function for V bf ij that depends on the penetration volume V p . Further, A i wbf = 2πr w l is the internal interfacial area between water and biofilm, C ij is the concentration of nutrients within the tube, E s is a saturation constant, k 1 is a growth rate constant and A i T = 2πRl is the external area of the tube. The ratio between the interior interfacial water biofilm area A i wbf and the external area of the tube A i T is a measure of the biofilm growth within the tube. If this ratio is zero, then there is no biofilm in the tube or the tube is full with biofilm. This means that if the tube is entirely filled with biofilm, then, interior growth stops since there is no more space in the tube. The sigmoid-like function is defined as, Note that the functionf = V p f (V bf ij ) is an increasing function of V p . This means that if V p increases the biofilm growth rate increases.
We write the area A i wbf in terms of the total volume of the pore V T and the volume of biofilm V bf , which gives, If there is no initial biofilm in the tube, the interfacial area between water and biofilm area is zero, therefore there is no biofilm growth in the interior of the tube.
Note that Eq. 7 represents a continuous relation between the biomass growth rate and the volume of biofilm V i bf ij . Secondly, we describe the biofilm that grows in the extremes of the tube. Since the penetration layer in the extremes is very small compared to the whole volume of biofilm, the biofilm growth rate is approximately proportional to the interfacial area between water and biofilm in the extremes A e wbf (see Fig. 4). We assume only interactions between nearest neighboring tubes. The interfacial area between water and biofilm A e wbf between the tube t ij and the tube t jk can be written in terms of the difference between volumes of biofilm of these neighboring tubes. If the volume of biofilm V bf jk in the tube t jk (connected to the node n j ) is larger than the volume of The side view of two neighboring tubes is shown. b The cross-section of the boundary between the neighboring tubes t jk and t ij is shown. The biofilm growth rate from t jk to t ij is proportional to A e wbf biofilm V bf ij in the tube t ij , then the biofilm grows in the extreme of the tube t jk and it is given to the neighboring tube t ij . The biofilm growth in the extreme of the neighboring tube t jk is given by Here, A T is the cross-sectional area of the tube. The ratio between the external interfacial water biofilm area A e wbf and the cross-sectional area of the tube A 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 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 ij . 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 biomass is lost in the tube. 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 wbf between the tube t ij and the tube t jk can be written in terms of the volume of the biofilm of the tubes. Hence, the equation for the biofilm growth at the extreme of the tube t jk reads as We take into account all the neighboring tubes whose volumes of biofilm are larger than the volume of biofilm in the tube t ij . To this extent, we introduce the following index set notation for the tube t ij which connects nodes n i and n j . Consider node n j , then we define the set of neighboring nodes of it, except n i by Λ ji . Therefore, the equation for the biofilm growth in the tube t ij due to biofilm growth in the extremes of the neighboring tubes is written as where We assume that the detachment of biomass is proportional to the area of the interface between water and biofilm, hence, detachment rate can be written in terms of the volume of biofilm as where k 2 is the detachment rate coefficient. Finally, when we take into account the interior growth, the growth in the neighboring tubes and the detachment of biofilm, the equation for the biofilm growth in the tube t ij can be written as Further, H (V bf ij ) 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 describe 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 nutrients consumption in the tube t ij is the result of the interior biofilm growth and biofilm growth in the extremes of the tube.
The equations that describe the transport of nutrients (3) and the biofilm growth (12) hold for the rectangular and the cubic network. The boundary conditions are different due to the different domain of computations. We present the boundary conditions for 2D case and for the 3D case in the following paragraphs.
When the mass conservation (2) is applied to every node in the network, a linear system for the pressure arises. The boundary conditions for the pressure in the 2D case are: Here, L x is the size of the network in the x direction and L y the size in y direction.
The initial condition for the concentration of nutrients is, The boundary conditions for the concentration of nutrients are: The initial condition for the biofilm growth is: Our algorithm randomly chooses 4% of the tubes. We have chosen b 0 = 10 −4 [kg/m 3 ] for the initial tubes that were seeded with biofilm.
links the concentration of nutrients in the nodes and the concentration in the tubes.
For the 3D case, the boundary conditions for the pressure are: The initial and boundary conditions for the concentration are:

Numerical method
The equation for the transport of nutrients is solved for the concentration C i at each node of the network. The advection part is solved using first order upwind scheme and a time-implicit method for the time integration. Hence, the discretization of the advection part reads as, where Ω i = {j | q ij is directed towards the node n i }. The diffusion part is discretized using a time-implicit method for the concentration. However, the area used is from the previous time step, 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, Therefore, 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. 5).

Results
In this section, we describe the overall mechanism of biomass growth and its implications on network characteristics such as the porosity and permeability. In addition, in order to determine the conditions for upscalability to real reservoir dimensions, we study the influence of the numerical parameters such as the number of nodes and physical parameters such as the size of the network and the inlet concentration of nutrients on the dynamics of transport of nutrients and biofilm growth. In addition, we use a 3D cubic network to study the influence of the nodes in the z direction. At the end, we mention the requirements in terms of the Damköhler number for upscaling this microscopic model to a continuum-scale model.
Firstly, we present the results of the transport of nutrients and the biofilm growth process in the pore network. Initially, there is a biomass concentration b 0 = 1 × 10 −4 [kg/m 3 ] in 4% of tubes (Fig. 6a). When the biomass gets into contact with the nutrients, biofilm starts growing and Pressure per node Flux in the tubes

Concentration of nutrients
Concentration of biomass spreading to the neighboring tubes (Fig. 6b, c). Due to a high injection rate, the nutrients are distributed over the whole network shortly after the beginning of the process. Therefore, biofilm grows uniformly through the network and hence the nutrients are consumed (Fig. 6d). After several minutes, depletion of nutrients near the outlet of the network is observed. This is because the consumption of nutrients by the bacteria near the inlet is very high, most nutrients are unable to reach the outlet. Hence, preferential biofilm growth is observed near the inlet (Fig. 6e, f) and the biofilm developed in this area causes the plugging of the network. This implies that a heterogeneous end-state is reached if the inlet concentration of nutrients is not large enough. The heterogeneous or preferential growth depends on parameters like the size of the network or the number of nodes in the network; therefore, the relation among the fraction of biomass and permeability varies with these parameters and upscaling is not possible in this case. However, if there is sufficient amount of nutrients, there is no depletion and the biofilm grows uniformly during the whole process. Therefore, the relation between fraction of biomass does not depend on the size of the domain of computation or the number of nodes and the problem is upscalable. This scenario is not shown in the Fig. 6.

Number of nodes
In this first set of simulations, we study the effect of the number of tubes in the network. The length in x and y direction and the initial porosity are constant for this simulations. We perform four simulations in which we use four different networks with different number of nodes (25 × 15, 50 × 30, 100 × 60, and 200 × 120). The length of the tubes decreases as the number of nodes increases. In order to keep the initial porosity constant throughout all four simulations, the radii of the tubes are adjusted in each of the simulations. Note that for each simulation, the radius of the tubes is constant over the network. The inlet concentration for these simulations is C in = 1[kg/m 3 ]. The value of the radii and the length of the tubes for this set of simulations are shown in Table 1. The complete list of parameters is shown in Table 2. 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). The evolution of the normalized flux through the network for the detachment rate k 2 = 0 [1/s] is shown in Fig. 7. We observe that there are some deviations among the curves; however, the network is plugged in all the cases around 300 min an even an S-shape is observed for very large number of tubes. The fraction of biofilm volume in the network is given by where V bf ij is the volume of biofilm in the tube t ij and V T is the volume of the tube t ij . In Fig. 8, we present the results obtained for the evolution of the fraction of biofilm volume in the network for k 2 = 0[1/s]. At the beginning of the process, the biofilm starts growing and spreading through the network, then, the biomass grows uniformly and in the final stage a preferential growth near the inlet causes the plugging of the network. We observe that the fraction of biomass necessary to block the networkS b decreases as the number of nodes increases. We can explain this by drawing an analogy with the minimum fraction of biomass needed to block one column of tubes in y direction. In order to keep the same initial porosity, if the number of nodes in the x and y direction doubles, then the radii of the tubes and the length reduces approximately to the half. Therefore, the volume of biomass in one column of the tubes decreases by 1/4 when the number of nodes in the y direction doubles. The sum of the volume of all the pores in the network decreases approximately by 1/2 when the number of nodes increases. Therefore the fraction of biomass needed to block the network decreases by a factor of 1/2. Note that in order to keep the same porosity, the total volume of the network, V nt = L x × L y × 2R decreases by 1/2 when the number of nodes increases. Hence, the fraction of biomass needed to block the network decreases when the number of nodes increases. In our simulations, the preferential growth is taken over more than one column; however, a similar argument to explain the decrease ofS b is valid. All configurations have the same initial porosity  Fig. 7 The normalized flux for different number of nodes for a rectangular network In Fig. 9, we show the normalized flux versus the fraction of biomass. We observe that if the fraction of biomass remains small, then the curves coincide for the four cases; however, as the fraction of biomass increases, the curves deviate from the uniform growth and exhibit a steeper descend as the number of nodes increases.
In order to determine when the preferential growth deviates from uniform growth, and therefore when upscalability is possible, we study the effects of biofilm growth on the Damköhler number, which is defined as

Da =
Reaction rate Advective transport rate . (28) In this case, we compute the Damköhler number related to the advective rate because the transport of nutrients is Nx=50 Nx=100 Nx=200 Uniform growth Fig. 9 The fraction of biomass vs the normalized flux for different number of nodes for a rectangular network mainly determined by this process, since the Peclet number is larger than 10 1 from the beggining of biofilm growth to 250 minutes approximately. The Damköhler number for the entire network, is obtained by dividing the average of the reaction rate by the average of the advective rate, Here, the sum is taken over all the tubes in the network. In Fig. 10, the Damköhler number for a various number of nodes in the network is shown. We plot two horizontal lines that enclosed the Damköhler number at which a transition from uniform growth to preferential growth  occurred. The transition from uniform growth to preferential growth occurs at different times for different number of nodes. However, the Damköhler number of the transition ranges between the two horizontal lines for all the cases. For a Damköhler number less than 10 1 the biofilm grows uniformly for all the cases. Further, for Damköhler number greater than 10 3 there is no uniform growth and upscaling is not possible.

Size of the network
In the second set of simulations, we study the impact of the size of the computational domain. We perform four simulations in which we enlarge the domain of computation.
The ratio between the length in the x direction and the length in the y direction is constant. In order to keep the same initial porosity, the radius increases slightly while we increase the size of the network. The inlet concentration for these simulations is C in = 1[kg/m 3 ]. In Table 3, the size of the network and the radius for each network are shown. For each simulation, the tubes in the network have the same radius. In Fig. 11, the normalized flux as a function of time for each simulation is shown. We observed that the network is plugged after 300 min for all the cases. In Fig. 12, the fraction of biomass as a function We observe that when the size of the computational domain is larger, the fraction of biomass needed to block the network decreases (see Fig. 13). Note that the total amount of biomass increases when the size of the network increases. The minimal amount of biomass required to block the network is the volume of all the tubes in one column. When we increase the size of the computational domain, the amount of tubes in one column is doubled while the total amount of tubes is four times higher, hence the relative contribution to the volume from one column decreases when we increase the size of the network. The minimal amount of biomass to block the network tends to zero as the size of the network increases and it is lower than the percolation threshold for a rectangular network. On the other hand, there is a maximum of biofilm growth when there is uniform growth because in that case all the tubes have to be filled with biofilm in order to plug the network completely. A transition from uniform growth to preferential growth is observed as we increase the size of the network. In Fig. 14, the Damköhler number is shown for different sizes of the network. We observe that there is uniform growth and therefore upscalability when the Damköhler number is less than approximately 10 1 . Above 10 3 , there is a preferential growth and upscalability is not possible.

Inlet concentration of nutrients
In the third set of simulations, we study the effects of the inlet concentration on the biofilm growth. In this set of simulations, the size of the network is 100 × 60 nodes and we use five different inlet concentrations, Fig. 15, the normalized flux is shown. We observe that the network is plugged at around 300 minutes for all the cases. In Fig. 16, the evolution of the fraction of biomass is shown. It is shown that the biomass saturation value increases as the inlet concentration increases. When the concentration C in = 1, the biomass saturation value is aroundS b = 0.70; however, if the inlet concentration is C in = 50, then the fraction of biomass necessary to block the network is S b = 1. Therefore, when we increase the concentration, the nutrients can reach the region near the outlet of the network and no preferential growth is observed; hence, the model predicts a uniform biofilm growth for concentrations larger than 25 [kg/m 3 ] (see Fig. 17). In Fig. 18, the Damköhler number is shown for different inlet concentrations of the network. We observe that there is uniform growth and therefore upscalability when the Damköhler number is lower than 10 1 . Further, we observe that when the inlet concentration C in = 25 and C in = 50, there is always uniform growth and hence upscalability.

Results for a 3D cubic network
In this section, we present the results obtained for the biofilm growth model in a 3D cubic network. The number of nodes in the x direction is N x = 50 and the number of nodes in y direction is N y = 30. The number of nodes in the z direction, N z , was varied to study the influence of the number of nodes on the porosity-permeability relation. The inlet concentration of nutrients was set constant and the radius is the same for all the tubes in the network. We performed five simulations with a different number of nodes in the z direction, N z = 2, N z = 4, N z = 8, In Fig. 19, the normalized flux as a function of time is shown. We observe that the flow through the network ceases after approximately 300 min for all the cases.
In Fig. 20, the evolution of the fraction of biomass in the network S b is shown. When the network is blocked,S b is around 0.7 for the 2D case. For the 3D cases, this fraction is larger: if N z = 2, thenS b is around 0.8 and as the number of nodes in the z direction increases, this fraction converges to 0.8510 approximately. The order of convergence computed is α = −0.93 approximately. Figure 21 shows the relation between the normalized flux and S b for the 2D case, the 3D cases and the case of uniform  biofilm growth in 2D. Note that the uniform biofilm growth in 2D and 3D yield similar results because there is no flux in the z direction. The biofilm growth in the 2D case deviates from uniform growth at around S b = 0.5. For the 3D cases, the deviation from uniform growth occurs at around S b = 0.7. We also observe this behavior in Fig. 22 where the average Damköhler number is plotted as a function of time. For the 2D case, there is a transition between uniform growth and preferential growth between 70 and 100 min. For the 3D cases, this transition is between 120 min and 160 min. However, the transition from uniform growth to preferential growth takes place in the same Damköhler regime.

Uniform biofilm growth
In the uniform biofilm growth, the biomass grows at the same rate in all the tubes of the network. Using a similar process as the equivalent resistance in electric circuits, we can easily obtain a relation between the total flux in the network and the fraction of biomass. Since the radii are reduced in all the tubes at the same rate and assuming that initially the radii are equal in all the tubes of the network, it follows that during the whole uniform biofilm growth process the radius remains the same in all the tubes of the network. This implies that there is only flux in the horizontal tubes of the network. Hence, the total flow in the network can be written as, in which N y is the number of tubes in the y direction and N x is the number of tubes in x direction. This expression is equivalent to consider the whole network as a single tube. If we disregard the second term in Eq. 30 and we express r 4 w ij as the fraction of volume of biofilm, the normalized flux for uniform growth is given by, In which S b0 is the initial fraction of volume of biofilm in the network. In terms of porosity, the normalized flux is given by the following equation, where φ 0 is the initial porosity. We use Darcy's Law to relate the permeability K with the flux Q, If the pressure drop ΔP , the cross-sectional area of the network A, the length in the flux direction L, and the viscosity μ are constant, we have that K K 0 = Q Q 0 = Q n . In which K 0 is the initial porosity and Q 0 is the initial flux. Therefore we can express (32) as, which relates the porosity and the permeability for the uniform biofilm growth.

Discussion and conclusions
In this work, we study the conditions for upscalability of bioclogging using a pore network model. We use a biofilm growth model that takes into account the spreading of the biofilm through the network and assumes that the consumption of nutrients is taken only in a small layer of biofilm adjacent to the water biofilm interface since the consumption of nutrients is faster than the diffusion of them through the biofilm in each tube. The biofilm growth requires the injection of nutrients through the network which are transported within a fluid phase. In general, it is shown that initially the biofilm grows uniformly across the network but afterwards there is a preferential growth near the inlet of the network, due to depletion of nutrients in the back of the network. This causes the plugging of the network and blocks the nutrient inflow. The amount of nutrients needed to clog the network depends on factors like the size of the network, the number of nodes, and the inlet concentration of nutrients. Therefore, the analytic relation between fraction of biomass and normalized flux may not be unique and upscalability is not always possible. However, if the inlet concentration of nutrients is about 25 [kg/m 3 ] there is no preferential growth and the biofilm grows uniformly through the network. This analytic relation does not depend on the size of the network, the number of nodes, or the inlet concentration; therefore, upscalability is possible. In the case of uniform growth, there is a unique relation between the fraction of biomass or porosity and the permeability of the network. We use the Damköhler number to determine when the biofilm grows uniformly through the network and therefore when upscalability is possible. We found that if Da < 10 1 the biofilm grows uniformly through the network. However, if the Da > 10 3 , there is preferential growth and therefore no upscalability is possible. If 10 1 < Da < 10 3 , there is a transition between uniform and preferential growth.
We performed four sets of simulations to determine the conditions for upscalability of bioclogging. In the first set of simulations, we vary the number of nodes in the network, we observe that at early stages the biofilm growth is mostly uniform for N x = 25 and N x = 50. However at the end of the process, a heterogeneous biofilm growth in the network is observed. This phenomenon does not allow upscalability for the whole process. For N x = 100 and N x = 200, the preferential growth is more significant than in the previous case. For this set of simulations, it is shown that if the Damköhler number is lower than 10 1 approximately, uniform growth is still observed in all the cases. However, if the Damköhler number is larger than 10 3 , the biofilm growth model starts to deviate from uniform growth and upscaling is not possible. Note that increasing the number of nodes leads to a decrease in velocities and hence to a increase in the Damköhler (see Eq. 29).
In the second set of simulations, we vary the size of the network. We observe that the fraction of biomass needed to block the network decreases when the size of the network increases. In this case, similar to the first one, there is not a unique relation between the fraction of biomass and the permeability of the network for the entire process. The Damköhler number predicts uniform growth below 10 1 approximately and preferential growth above 10 3 . We see from Eq. 29 that if we increase the size of the network, the Damköhler number also increases.
For the third case of simulations, we vary the inlet concentration of nutrients. We observe that for inlet concentrations C in = 25 [kg/m 3 ] and C in = 50 [kg/m 3 ], the fraction of biomass needed to block the network is approximately equal to one. For larger inlet concentration of nutrients, there is no depletion and hence upscaling is possible since there is uniform biofilm growth. In this case, if we increase the inlet concentration the Damköhler number decreases.
In the last set of simulations, we obtained that for the z = 1 (2D case), the depletion of nutrients occurs faster. This implies that the biofilm growth deviates from uniform growth earlier in the 2D case than in the 3D case. The depletion occurs faster in 2D since in 3D there are more ways in which the water can flow, hence, the nutrients are able to travel further in a 3D network and depletion takes longer. The amount of biomass needed to block the network is larger for the 3D cases than for the 2D case. However, as the number of nodes in the z direction increases, the amount of biomass needed to block the network converges to a limit value.
Even though the depletion of nutrients occurs faster for the 2D case, the transition between uniform growth and preferential (or heterogeneous) growth occurs for the same values of the Damköhler number, which is between Da = 10 1 and Da = 10 3 . Therefore, we conclude that the criteria for uniform growth and therefore upscalability in 3D are similar to the criteria used for the 2D network.
In 2D, we performed three sets of simulations in which we vary the number of nodes in the network, the size of the network and the inlet concentration of nutrients. For the first two cases, upscalability is not possible since there is no a unique relation between the amount of biofilm and the permeability of the network. However for the third case, when we vary the inlet concentration of nutrients, we observe that for concentrations larger than 25 [kg/m 3 ] the model describes uniform biofilm growth in the network, which allows upscalability of these results, since in uniform biofilm growth the relation between the fraction of biomass and the permeability does not depend on the volume of the network. In addition, we show that if the Damköhler number is less than approximately 10 1 , the biofilm evolves similarly to uniform growth and that if it is above 10 3 , preferential growth is observed; therefore, we can use the Damköhler number to determine whether upscaling is possible. For the first two cases, the Damköhler number is not always below this limit and therefore upscalability is not possible in these cases. However, for the third set of simulations, we observe that the Damköhler number is below this limit for the entire process for an inlet concentration of C in = 25 [kg/m 3 ] and C in = 50 [kg/m 3 ]. For the upscalable case, we obtain a relation between the permeability and the porosity, K ∼ φ 2 . This formula can be seen as an alternative to the classical Kozeny-Carman equation in the case of gradually clogging of the network.
Interesting research may be the extension of this model to 3D networks with different topologies to determine the effects of biofilm growth on the relation between porosity and permeability. This relation could differ from K ∼ φ 2 . In addition, it might be interesting to verify the relation between porosity and permeability, K ∼ φ 2 , in laboratory scale and obtain an appropriate Damköhler number regime for uniform growth. Forthcoming research might be the extension of this model to two phase flow for studying the possibility of flow diversion for MEOR. Finally, this model can be used in other problems like pore-elasticity or in physiological situations such as the modelling of clogging of arteries in the brains or on the heart due to arteriosclerosis.