Access the full text.
Sign up today, get DeepDyve free for 14 days.
Mathematical Modeling, Laboratory Experiments, and Sensitivity Analysis of Bioplug Technology at Darcy Scale David Landa-Marb an, Gunhild Bdtker, Bartek Florczyk Vik, and Per Pettersson, Norwegian Research Centre; Iuliu Sorin Pop, University of Hasselt; Kundan Kumar, Karlstad University; and Florin Adrian Radu, University of Bergen February 4, 2020 Summary In this paper we study a Darcy-scale mathematical model for bio lm formation in porous media. The pores in the core are divided into three phases: water, oil, and bio lm. The water and oil ow are modeled by an extended version of Darcy's law and the substrate is transported by diusion and convection in the water phase. Initially there is bio lm on the pore walls. The bio lm consumes substrate for production of biomass and modi es the pore space which changes the rock permeability. The model includes detachment of biomass due to water ux and death of bacteria, and is implemented in MRST. We discuss the capability of the numerical simulator to capture results from laboratory experiments. We perform a novel sensitivity analysis based on sparse-grid interpolation and multi-wavelet expansion to identify the critical model parameters. Numerical experiments using diverse injection strategies are performed to study the impact of dierent porosity-permeability relations in a core saturated with water and oil. Introduction After primary and secondary production, up to 85% of the oil remains in the reservoir [23]. Microbial improved and enhanced oil recovery (MIEOR) is one of the secondary and tertiary methods to increase the oil production using microorganisms [35]. Bioplug technology is a MIEOR strategy that consists in plugging the most permeable zones in the reservoir, which provokes water to ow through new paths and recovering the oil in these new zones. Experiments in microsystems allow us to observe processes in more detail, which leads to improvement of the experimental methods in core-scale experiments. For example, in [21] the eects of ow velocity and substrate (also referred to as nutrients) concentration on bio lm in a microchannel was studied, nding values of substrate concentration and ow velocity for a strong plugging eect. Core samples from reservoirs can be used to study changes in permeability due to bio lm formation, e.g., in [30] two-phase ow experiments were performed to study the selective plugging strategy for MIEOR. In that study, the MIEOR eects increased the oil recovery around 25%. Mathematical models of bioplug technology are important as they help to predict the applicability of this MIEOR strategy and to optimize the bene ts. In [10] a mathematical model for single-phase ow was proposed which includes changes of rock porosity and permeability as a result of bio lm growth. Li et al. [18] built a mathematical model for two-phase ow including the eects of bio-surfactants and biomass on improving the oil recovery. The authors also compare the numerical results for dierent porosity-permeability relations. These porosity-permeability relations can also include the permeability of bio lm and be derived as a result of upscaling pore-scale models [32, 8, 15]. In this work, we present a two-phase core-scale model of bioplug technology and study the oil production for dierent porosity-permeability relations and injection strategies. Sensitivity studies of mathematical models are of great interest because they provide estimates of the in uence of physical parameters on a quantity of interest, e.g., bio lm formation. In [3] a regional steady-state sensitivity analysis was performed to identify parameters with the largest impact on a mathematical model for deammoni- cation in bio lm systems. Sensitivity analysis by means of Sobol decomposition provides rigorous estimates of parameter dependencies, but are prohibitively expensive to compute if the number of parameters is large. This is remedied for smooth problems by rst computing spectral (generalized polynomial chaos) expansions in the parameters, which then leads to ecient evaluation of the sensitivity indices via post-processing of spectral co- ecients [29]. The latter method was employed in [14], where a global sensitivity analysis was performed using Sobol indices to identify the critical parameters of a pore-scale model for permeable bio lm. In this paper, we consider nonsmooth models in the dependent parameters, hence spectral expansions with global smooth basis arXiv:2002.00090v1 [physics.app-ph] 18 Dec 2019 functions are not a robust choice. Instead, we propose a two-stage method where we rst use sparse grids to estimate a piecewise linear interpolant of the function of interest, which yields a surrogate that can be further sampled at negligible cost. Secondly, we compute a multi-wavelet representation of the interpolant, from which the sensitivities can be directly evaluated. The mathematical model consists of coupled nonlinear partial dierential equations. Two-point ux approx- imation (TPFA) and backward Euler (BE) are used for the space and time discretization respectively. To solve the resulting nonlinear algebraic system, Newton's method is used. The scheme is implemented in the MATLAB reservoir simulation tool (MRST), a free open-source software for reservoir modeling and simulation [19]. To summarize, the new contributions of this work are: • The comparison of experimental results to numerical simulations of a core-scale porous medium including bio lm formation. • Performing a novel global sensitivity analysis of the model parameters to identify the critical parameters. • Performing numerical simulations for dierent porosity-permeability relations to study their impact on the predicted oil recovery. • The study of diverse substrate injection strategies for the bioplug technology. The paper is structured as follows. Firstly, we introduce and describe the implementation of the core-scale mathe- matical model for two-phase ow including the eects of bio lm formation. Secondly, we present a comparison of numerical simulations to laboratory experiments of this core-scale model. Thirdly, we introduce a novel method for global sensitivity analysis and apply it to the mathematical core-scale model. Diverse numerical experiments for dierent porosity-permeability relations and injection strategies are also explained. Finally, we present the conclusions. Core-Scale Model We consider a core sample of radius r and length L initially lled with oil and water and a given bio lm distri- bution. Fig. 1 shows the schematic representation of this system. As water and substrate are injected, some biomass is detached due to erosion (shear forces caused by the water ow). The bio lm consumes substrate to produce metabolites (i.e., gases) and to grow which modi es the rock porosity and hence the rock permeability. Consequently, the ow pattern is modi ed in the sense that pores where oil was replaced by water, and which were forming a preferential path for water ow, become less permeable. Therefore water enters other pores, mobilizing the oil present there and leading to improved oil production. The mathematical model presented in this section aims to describe the following processes in a core after bacterial inoculation: two-phase ow, substrate transport, permeability changes due to porosity modi cation, and bio lm growth, detachment, and death. Oil r Water Permeability Water Biomass Substrate Porosity Substrate Metabolites Fig. 1|Core sample for laboratory experiments after inoculation of bacteria. We assume that the uids are immiscible and incompressible. We consider that the bio lm has water and biomass as the only two components, and that the bio lm porosity and bio lm permeability are constants. In order to determine the amount of uid outside the bio lm in the representative element volume (REV), we introduce the saturation of a uid S (for oil = o and for water = w) given by the ratio of volume of uid outside the bio lm over the volume of voids outside the bio lm (in REV). For this two-phase ow model, we have that S + S = 1. o w The mass conservation and extended Darcy's law for the oil phase are given by @ k r;o ( S ) +r ( ~v ) = Q ; ~v = k(rp ~g); : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (1) o f o o o o o o o @t and for the water phase @ k r;w [ ( S + )] +r ( ~v ) = Q ; ~v = k(rp ~g); : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (2) w f w b w w w w w w w @t 2 where is the porosity outside the bio lm, ~v the ow velocity, k the relative permeability, the uid density, f r; k the absolute rock permeability, ~g the gravity, the viscosity, the volume fraction of bio lm in the REV, the bio lm water content, and Q source/sink terms. The minimum value of is zero while the maximum w b value of is equal to the initial porosity of the rock , corresponding to the case when the pores are lled with b 0 bio lm. In this work, we assume that there are no uid sources/sinks, neglect the gravity eects, set to zero the residual oil saturation and irreducible water saturation, assume oil cannot penetrate into the bio lm, and neglect the capillary pressure (p = p ). The previous assumptions are commonly used to reduce the complexity and w o number of parameters of mathematical models for porous media. Consider the mass conservation for the water in Eq. 2. The rst term gives the changes on time for the total water mass in the REV, i.e., the water outside the bio lm ( S ) plus the water inside the bio lm ( ). The Brook-Corey relations are commonly used w f w w b w to model the relative permeability relations. These relations are given as a function of the water saturation and experimental parameters such as exponents for the saturation, values for the endpoint relative permeabilities, and residual phase saturations. Simpli ed expressions of these relations are given by k (S ) = S ; k (S ) = (1S ) ; : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (3) r;w w r;o w w where and are experimentally tted factors. We assume that water can ow inside the bio lm even when there is only oil outside the bio lm (S = 0). Then, the previous relation is not a good model because for S = 0 w w it results in k (0) = 0 which leads to zero water velocity. Recalling that the relative permeabilities are used r;w to extend Darcy's law to multiphase ow, we look for relationships which account for the water saturation and volume fraction of bio lm. These relationships need to ful ll the following criteria: when there is only oil outside the bio lm, then k (S = 0; ) = 1; when there is only bio lm in the REV, then k (S ; = ) = 1. Then, r;w w b r;w w b 0 we propose the following two relative permeability relations: b b b k (S ; ) = S 1 + ; k (S ; ) = (1S ) 1 : : : : : : : : : : : : : : : : : : : : : : : : : : : : (4) r;w w b r;o w b w 0 0 0 To describe the movement of substrate, we consider the following convection-diusion-reaction transport equations: ~ ~ [C ( S + )]+rj = R ; j = D ( S + )rC +C ~u : : : : : : : : : : : : : : : : : : : : : : : : : : : : (5) n f w b w n n n n f w b w n n w @t In the previous equations, C is the substrate concentration in the water, j the substrate ux in water, and D n n n the substrate dispersion coecient which includes mechanical dispersion plus diusion. The reaction term R is given by R = ; : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (6) n n b b K + C n n where is the maximum rate of substrate utilization, the biomass density, and K is the Monod half-velocity n b n coecient. The following equation describes the bio lm evolution: @ C ( ) = K K jjrpjj : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (7) b b n b b d b b str b b @t K + C n n Here, we assumed that all consumed substrate is used to produce biomass, a linear death of bacteria given by K , and erosion of bacteria given by the magnitude of the pressure gradient times a constant depending on the bio lm K . After biomass has been detached, in this model we assume that the detached biomass ows out of str the core and does not aect the rock properties; therefore, we do not include a transport equation for the detached biomass. Diverse porosity-permeability relations have been proposed for the last decades. We refer to [8] for a recent review of these relations. Three porosity-permeability relations commonly used in modeling are the following: 3 2 k k k p f vp f crit h f crit f crit = ; = ; = a + (1a) : : : : : : : : : : : : : : : : : : : (8) k k k 0 0 0 0 crit 0 0 crit 0 crit The rst one is called the power law [4, 9], where is a tting factor calibrated either from experimental data or taken from a process-speci c literature. The second relation is a variation of the power law known as the Verma-Pruess relation [34], where is a critical porosity when the permeability becomes zero, which value is crit between 70 and 90% of the initial porosity. The third relation is proposed by [31], where a is a weighting factor 3 between -1.7 and -1.9. These three relations do not include the bio lm permeability. Vandevivere [33] proposes the following relation of permeability and porosity for a plugging model " # ( " #) 2 2 2 k 1 B 1 B k v r f r b 0 = exp + 1 exp ; : : : : : : : : : : : : : : : : : : : : : (9) k 2 B 2 B k (k k ) 0 c 0 c 0 0 0 b f where k is the bio lm permeability, B is a relative porosity given by B = 1 = , and B is the critical point b r r f 0 c where bio lm begins to detach and form plugs. Thullner et al. [31] presents the following relation which includes the bio lm permeability: k k k th f crit b 0 = + : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (10) k k k + k 0 0 crit 0 0 b All these relationships are postulated at the Darcy-scale, based on experimental observations. A dierent ap- proach is followed in [32], [26], and [15], based on homogenization. Starting with models at the pore scale, one applies expansion methods to derive the mathematical models valid at the Darcy scale. In [15] eective porosity- permeability relations are derived for two dierent pore geometries: thin channels k and tubes k . In this work c t we compare the following four porosity-permeability relations: Vandevivere k (9), Thullner k (10), channel k , v th c and tube k . We refer to Appendix A for the mathematical expressions of k and k . t c t The porosity in the porous medium changes in time as a function of the bio lm volume fraction (t). When there is no bio lm, the porosity (t) is equal to the initial porosity. In the case when the porous medium is lled with bio lm ( (t) = 0), the porosity in the REV is equal to the bio lm porosity . The porosity-permeability f w relations (9-10) are given as a function of the rock porosity and do not include the bio lm porosity. The porosity- permeability relations k and k do include the bio lm porosity. Using the de nitions of (t) and (t), we have c t f b the following relation for the void space outside the bio lm (t) = (t): : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (11) f 0 b This mathematical model is implemented in a 3D domain with a uniform cell-centered grid. The TPFA is used for the spatial discretization, while BE for the time discretization. The Newton's method is used as linearization scheme. The implementation of this core-scale mathematical model is done in MRST, based on the modi cation of the polymer example in the module ad-eor (see [2]). Model Test Core-scale experiments under controlled conditions are performed in the laboratory for studying the eect on bio lm growth in porous media. These experiments aim to provide a better understanding of dierent features, i.e., the relation between bio lm composition and growth conditions, the plugging potential of dierent bacteria, and the adaptability of bio lm at diverse substrate ux. Fig. 2 shows a photograph of a typical core sample to study these eects. Fig. 2|Core sample for laboratory experiments. In this section, we describe an experiment performed by NORCE for Equinor. The aim of this experiment is to test dierent brine qualities on the same system and the same bio lm. The core has a length of 29.50 cm and 4 a radius of 1.90 cm. Initially, the core has an approximately homogeneous porosity of =0.23 and permeability k =1528 md. The core sample is introduced inside a core holder, fully saturated with brine. Heating cables wrapped around the exterior of the core holder are used to conduct the experiments at a constant temperature (30°C). A backpressure regulator (BPR) is used to control the outlet pressure. Bacteria are injected into the core and the core is left standstill overnight to allow the bacteria to attach themselves to the pore walls. Fig. 3 shows a diagram of the experimental set-up. Δp Δp Δp Δp Δp Δp r=1.90 cm φ =0.23 L=29.50 cm k =1528 md Ambient room temperature Medium P Water Eﬄuent BPR Fig. 3|Scheme of an experimental set up for laboratory experiments. The core is injected with two dierent brines with dierent substrate concentrations, where glucose is added as a 3 3 source of carbon. We refer to these brines as B =30:81 kg/m and B =42:63 kg/m . The two dierent brines are 1 2 10 3 injected at a constant ow rate of 8.3310 m /s (Darcy velocity of 6.4 cm/day). The injection strategy is the following: B is injected during 53 days, after B is injected during 82 days and nally B is injected again during 1 2 1 61 days. The core permeability is estimated at dierent times, using the measurement of the pressure dierence along the core divided by the pressure drop on the clean core. This ratio is known as resistance factor R and it is given by R (t) = p (t)=p . Assuming that there is not change in the uid rate, viscosity, and density, the f w 0 resistance factor gives an estimate of the current rock permeability k(t) = R (t)k . f 0 In the previous section we introduced the model equations for bio lm growth and two-phase ow in porous media. We simplify the mathematical model to compare to the experiment. This experiment can be modeled as a 1D single-phase ow system, where in addition we consider that the bio lm is impermeable and only diusion for the substrate to reduce the number of parameters. Table 1 presents this simpli ed version of the mathematical model which includes only six variables: the water velocity v , water pressure p , permeability k, porosity , w w substrate concentration C , and volume fraction of biomass . n b Name Equation Darcy v = (k= )@ p ; @ + @ v = 0 w w x w t f x w Permeability k = k ( = ) 0 f 0 Porosity = f 0 b Substrate @ (C ) + @ (v C D @ C ) = C =(K + C ) t n f x w n n f x n n b b n n n Bio lm @ = C =(K + C ) K K j@ p j t b n b n n n d b str x w b Table 1|Core-scale single-phase equations. Nine parameters are needed to solve the mathematical model in Table 1. To obtain a better estimate of the model parameters, it is necessary to perform various experiments under controlled input quantities. However, these experiments are expensive and time consuming. In [14] a pore-scale mathematical model is calibrated based on the experiments performed by [21], where measurements of the bio lm amount over time are taken for dierent ux velocities. For the core-scale system described in this section, only one experiment is performed. Then, we select parameter values from the literature in order to run numerical simulations and compare qualitatively the results to the experimental observations. Table 2 shows the selected values of these parameters. These parameter values have the same order of magnitude as the ones used in mathematical modeling [1, 6, 8, 14]. The stress coecient K is neglected before the simulations and its value is chosen to better t the experimental str measures. In addition to the model parameters and core dimensions, initial, and boundary conditions are needed 5 Parameter Notation Value Reference 6 1 Bacterial death rate K 3.210 s [12] 6 1 Maximum growth rate 510 s [18] Monod-half velocity K 0:915 kg=m [20] 10 2 Substrate diusion coecient D 510 m =s [7] Biomass density 20 kg=m [25] Power law constant 2:5 [8] Water dynamic viscosity 10 Pas Standard 3 3 Water density 10 kg=m Standard Table 2|Model parameters for the veri cation study. to complete the model. The initial porosity and permeability are =0.23 and k =1528 md respectively. We 0 0 assume that initially there is no substrate in the brine. The pressure is set to zero initially (in this system, we are only interested in pressure dierences that aects the detachment and transport of substrate). Recalling that the time after inoculation and before starting the substrate injection is roughly half day, we assume that the initial volume fraction of bio lm is distributed uniformly along the core and has a value of (xxx; 0) = 10 . Then, we performed numerical simulations using the same injection strategy of substrate. After simulations, the value of the stress coecient K that ts best the experimental data is 1.510 m/Pas. Fig. 4 shows the experimental str measurements of the average resistance factor for the core sample and the numerical results over time. Experiment Simulation B1 B2 B1 0 20 40 60 80 100 120 140 160 180 Time (days) Fig. 4|Experimental and simulated resistance factor R over time. From the numerical simulation and experimental data, we observe that the resistance factor changes over time, which means that the permeability is aected by the bio lm. We can distinguish between the three dierent 6 periods of substrate injection. For the injection of B , both simulations and experimental measures show that the resistance factor does not increase signi cantly. When increasing the substrate concentration to B , we observe that the resistance factor increases due to major bio lm activity. When the substrate injection is set back to B , we observe a steady state where the resistance factor has increased by a factor of three. From the simulation, it is possible to observe how the resistance factor decreases from changing the substrate injection B to B , which means 2 1 a sensible response of the bio lm to substrate input. One of the outcomes of this experiment is that the plugging can be controlled by modifying the substrate ux. In addition, the growth conditions of the microorganism are aected by the dierence physical environments. From the experimental data we observe uctuations on the measurements. This behavior is attributed to the dynamical attachment and detachment of biomass, a complex process not included in the mathematical model. However, this mathematical model is simple enough to have few parameters and complex enough to capture processes such as dynamical changes on the resistance factor and a steady state where the bio lm growth is balanced with the detachment and death of bacteria. Fig. 5 shows values of the simulations for pressure, substrate, and biomass along the core. We observe from subplot (a) that after a day of substrate injection the constant ow rate B requires an inlet pressure of approximately 200 Pa. From subplot (b) we observe that the inlet pressure has increased more than twice in order to keep the ow rate, as a result of permeability reduction due to the bio lm formation. From the substrate pro le (c), this decreases from the inlet to the outlet as bacteria consume them. Given the initial homogeneous volume fraction of bio lm, this has increased more than three order of magnitude in comparison with the initial volume fraction, but it has stopped growing due to the shear forces. Fig. 5|Pressures, substrate, and biomass pro les along the core. Global Sensitivity Analysis We keep the exposition on sensitivity analysis general to emphasize that the methodology is not restricted to the problems presented in this paper. Let q(~y) be a scalar multidimensional function with ~y = (y ; : : : ; y ) 2 R , 1 n de ned by a range of independent parameters, i.e. = (a ; b ) ::: (a ; b ). We associate a nonnegative weight 1 1 n n w (y ) = 1=(b a ) with each parameter y 2 [a ; b ] for j = 1; : : : ; n, and use the product measure w = w . j j j j j j j j j=1 More general weight functions are possible, including non-product measures corresponding to inter-dependence between the parameters [24]. In this paper we do not aim to quantify uncertainty by stochastic models, but sensitivities in outputs given bounds on the input parameters. To determine the relative eect of each of the n input parameters on the output quantity of interest q, we perform a global sensitivity analysis in terms of a Sobol decomposition [28]. The function q is decomposed as a series expansion in all subsets of variables; the sum of contributions from the individual variables in isolation, all combinations of pairs of variables and so on, leading to the expression: n n X X f;g fig fi;jg f1;:::;ng q(~y) = q + q (y ) + q (y ; y ) + : : : + q (~y); : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (12) i i j i=1 i=1;j>i 7 where the terms are de ned recursively by f;g q = q(~y)w(~y)d~y; : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (13) fig f;g q (y ) = q(~y)w (~y )d~y q ; 1 i n; : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (14) i i i i fi;jg fig fjg f;g q (y ; y ) = q(~y)w (~y )d~y q (y ) q (y ) q ; 1 i < j n; : : : : : : : : : : : : (15) i j i;j i;j i;j i j i;j and so on for higher-order terms. The subscript i means that the ith index is omitted, e.g., w = w . i j j6=i The Sobol index for the s-parameter combination fy ; y ; :::; y g is given by i i i 1 2 s fi ;:::;i g 2 1 s S = [q (y ; :::; y )] w (y ):::w (y )dy : : : dy : : : : : : : : : : : : : : : : : : : : : : (16) fi ;:::;i g i1 is i1 i1 is is i1 is 1 s Var(q) i ;:::;i The total variability in q due to variable i is obtained by summing over all subsets of parameters including parameter i, denoted I , which yields the total Sobol index for parameter i, X X S = S : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (17) fig fi ;:::;i g 1 s s=1 i2fi ;:::;i g 1 s Directly computing global sensitivities of a computationally complex function, e.g., a computer model, over a domain of only moderately high dimensionality is often infeasible due to the computational cost. As an additional challenge, we are interested in nonsmooth q, excluding the use of Sobol decomposition based on generalized poly- nomial chaos for smooth problems [29]. As a remedy, we approximate the function of interest with an interpolant on a Clenshaw-Curtis type sparse grid, using the software in [11]. Subsequently, a multi-wavelet decomposition of the interpolant is performed, and the global sensitivity indices are evaluated directly from the multi-wavelet coef- cients. The accuracy in the sensitivity indices as determined by the computer model is thus largely determined by two kinds of errors: interpolation error and basis truncation error introduced when replacing the interpolant by a series expansion in a nite set of multi-wavelets. Multi-Linear Collocation. The outline of multi-linear collocation closely follows the exposition in [11], which also provides the source code for the numerical implementation. We perform interpolation of possibly nonsmooth (but continuous) functions and therefore rely on localized basis functions. For discontinuous functions, we refer to the adaptive sparse grid methods introduced in [22]. Multidimensional collocation on sparse grids is built from tensor products of low-order single dimensional collocation rules. Temporarily assume that n = 1, denote y = ~y and let Y be a set of interpolation nodes in the parameter space, for which the function q is evaluated. The index ` refers to a re nement level: the higher `, the bigger the set of interpolation nodes. An interpolant of q is given by I (q) = q(y ) (y); : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (18) ` ` y y 2Y ` ` where (y) are nodal basis functions with the property 1 if ` = ` ; (y 0 ) = (19) y ` 0 otherwise. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : In this work we employ the piecewise linear `hat' functions for level ` > 1 i1 1i 1 2 jy y j if jy y j < 2 ; ` ` (y) = (20) 0 otherwise, : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : supplemented with the unit function (level ` = 1). The collocation nodes are nested, i.e., Y Y . Rather than ` `+1 expressing the interpolant as a function of the nodes in some suciently re ned set Y , it can be expressed as a sum of hierarchical surpluses, i.e., dierences between the interpolants at successive levels, (q) = I (q) I (q); I (q) 0: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (21) ` ` `1 0 Then, Eq. (18) can be expressed ` ` X X X (j) I (q) = (q) = q(y ) (Y ); : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (22) ` j j j=1 j=1 y 2Y Y j j j1 8 where the superscript (j) is used to emphasize that the basis functions are level speci c. Using this construction, one may use the local dierence between the interpolant and the true function as a measure to determine where further re nement in terms of more nodes are needed. In multiple dimensions (n > 1), a tensor-product formula would lead to very large sets of nodes even for moderate n. As a remedy, a sparse grid based on Smolyak's construction [27] is used, adding only a subset of the tensor product Clenshaw-Curtis nodes at each new level. The hierarchical multidimensional interpolant of q = q(y ; : : : ; q ) on total level p n, is given by 1 n I (q) = I (q) + I (q); : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (23) p;n p1;n p;n with I = 0 and n1;n X X X ~ ~ (i ) (i ) (i) (i) 1 n I (q) = ( )(q) = q(y ) I (q)(y ) : : : : : : : (24) p;n i i p1;n 1 n j j ~ ~ 1 n j j | {z } | {z } ~ ~ ~ jij=p jij=q j ~ ~ (i) (i) ~ ~ j j For more details on the multi-linear collocation method we refer to [11]. The hierarchical multi-linear inter- polant (23) is a surrogate from which we can obtain approximations of q. The Smolyak algorithm yields the multidimensional interpolant at low computational cost, but the global sensitivity indices are not directly avail- able. Next, we introduce a series expansion of I in multi-wavelets to obtain the Sobol decomposition of I . p;n p;n Sobol Indices via Multi-Wavelet Expansions. An alternative representation of q is via a spectral expansion in a set of multidimensional orthogonal basis functions f' g, q(Y ) = c ' (Y ); : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (25) j j j=1 where the in nite sum is truncated to some nite set in practical application. This is known as generalized polynomial chaos expansion if the basis functions are orthogonal polynomials [36], and it has been demonstrated that the Sobol indices can be directly identi ed from the expansion coecients, once these have been computed [29]. Nonsmooth functions can be represented by means of multi-resolution analysis where the basis functions are piecewise polynomial multi-wavelets [17], and the identi cation of the Sobol indices is straightforward. A piecewise linear multi-wavelet expansion of sucient resolution can exactly represent the piecewise linear interpolant, provided that the support nodes are aligned with the support of the multi-wavelets. The advantage of using the multi-wavelet expansions rather than trying to directly evaluate the Sobol component functions, is that we may rely on a single interpolation of the function of interest itself (and not its square, conditional on some subset of the variables, and so on). Let f' (~y)g be a nite-dimensional basis of multi-wavelets indexed by some set of nonnegative integers ~ ~ k k2I MW I N . Ideally, the number of basis functions should be limited but chosen in order to be a good approximation MW of the interpolant I (q) of the function of interest. The approximation properties of the multi-wavelet basis are p;n determined by the order of the piecewise polynomial multi-wavelets, and the resolution level that governs the localization in parameter space. Analogous to the hierarchical surplus de ned for the multi-linear interpolant, the hierarchical surplus of the multi-wavelet expansion is given by the contribution captured by the nest resolution level. The multi-wavelet basis is hierarchical, and can be enriched by adding multi-wavelets to regions labeled important due to large surpluses. We may now either project the interpolant I (q) onto the basis f' (~y)g , or interpolate the product p;n ~ ~ k k2I MW q' , and then perform the projection by computing the expected value. The latter is simpler, as it only involves integration of a piecewise linear function to compute each multi-wavelet coecient. However, we do not know a priori what multi-wavelets should be included in the truncated basis, and will therefore settle for a simple strategy that also admits adaptivity. Unlike q itself, the interpolant I (q) can be sampled repeatedly at moderate p;n computational cost. Given N random samples of the parameters ~y, and a tentative multi-wavelet basis of size P , we form the linear system for the P -vector of multi-wavelet coecients ~c ~c = ~q; : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (26) NP (k) N where the matrix 2 R is de ned by [] = ' (~y ), and ~q 2 R contains the evaluations of I (q) in the k;j j p;n samples y , j = 1; : : : ; N . The ordinary least squares solution to (26) gives the multi-wavelet coecients of the tentative bases. If the hierarchical surplus (de ned as the dierence between the interpolants of two successive grid levels) is above a user speci ed threshold, the approximation of I (q) is not suciently resolved, and we p;n may increase the tentative multi-wavelet basis functions in unresolved regions of parameter space. This is done by solving (26) with an extended basis. As a further measure of how well we represent I , we may compare the p;n variance predicted by the multi-wavelet expansion with the sample variance of I (q). p;n 9 Numerical Experiments For the numerical experiments, we consider a core-sample saturated with water and oil. We impose a more permeable zone (10 k ) in the middle of the core, as shown in Fig. 6a. This more permeable zone is set to study the eects of bio lm on oil recovery when it modi es the rock properties in a thief zone. We consider that initially there is only bio lm in the middle part of the thief zone with a volume fraction of 2.510 , as shown in Fig. 6b. Fig. 6|Initial rock permeability (a) and initial volume fraction of bio lm (b). The core has an initial water saturation of 0.2 and an oil saturation of 0.8. The substrate is injected on the left side of the core at a velocity of 1 m/day. The initial water pressure is set to 10 Pa. Table 3 lists the parameters for the numerical simulations. These boundary conditions, initial conditions, and parameters are selected to perform the diverse numerical studies in a computational time of order of minutes. Parameter Notation Values Reference 6 1 Bacterial death rate K 3.210 s [12] 4 1 Maximum growth rate 1.5910 s [20] Monod-half velocity K 0:915 kg/m [20] 10 2 Substrate diusion D 510 m /s [7] Biomass density 20 kg/m [25] Factor Brooks-Corey ; 2 [13] Critical point B 0.35 [33] Fitting factor Power law 2.5 [8] Stress coecient K 10 m/Pas [16] str Critical porosity 0 Assumed crit 12 2 Initial permeability k 2.4510 m Assumed 13 2 Bio lm permeability k 2.4510 m Assumed Initial porosity 0.21 Assumed Core length L 0.30 m Assumed Core radius r 0.19 m Assumed Injected substrate concentration C 5 kg/m Assumed Injected water velocity v 1.1610 m/s Assumed Bio lm water content 0.9 Assumed Water viscosity 10 Pas Standard Oil viscosity 3.9210 Pas Standard 3 3 Water density 10 kg/m Standard Oil density 800 kg/m Standard Table 3|Table of input variables and model parameters for the veri cation study. Impact of Porosity-Permeability Relations. We compare the oil recovery for four porosity-permeability relations: Vandevivere k (9), Thullner k (10), channel k (Appendix A), and tube k (Appendix A). These four v th c t relations include the bio lm permeability. Fig. 7 shows the pro les of these four porosity permeability relations for a highly permeable bio lm and a less permeable bio lm respectively. From Fig. 7 we can observe that these permeability-porosity relationships give dierent values of permeability for the same values of rock porosity. Thus, we expect to observe dierent oil recovery predictions for each of these 10 1 1 Power law Power law Vandevivere Vandevivere 0.9 0.9 Thullner Thullner Channel Channel 0.8 0.8 Tube Tube 0.7 0.7 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Ratio of Initial to Reduce Porosity Ratio of Initial to Reduce Porosity (a) (b) 1 3 Fig. 7|Porosity-permeability relationships of kkk =10 kkk (a) and kkk =10 kkk (b). bbb 0 0 0 bbb 0 0 0 relationships. Fig. 8 shows the comparison of percentage of oil recovery in comparison with the initial oil in the core for the four dierent porosity-permeability relations. The bio lm permeability is set to k = 0:1k . b 0 We observe that after one day of water injection, the Vandevivere porosity-permeability relation predicts larger oil recovery than the tube relation. In addition, the channel relation predicts slightly more oil recovery than the Thullner relation. These results are expected from Fig. 7b, where the Vandevivere relationship gives the fastest permeability reduction as a function of the porosity reduction in contrast to the tube relationship which gives the slowest permeability reduction. From Fig. 8, we conclude that the predicted oil recovery diers for the dier- ent porosity-permeability curves. The oil recovery could be over- or underestimated depending on the assumed porosity-permeability relation. Injection Strategies. The study of injection strategies is important for the optimization of the bioplug technol- ogy. Comparison of oil recovery for dierent injection techniques is possible through numerical experiments. In this work, we compare the oil recovery for dierent substrate injection strategies for the same initial conditions. For this study, we use the k porosity-permeability relation. Fig. 9 shows a reference pro le where only water is injected and the ve dierent strategies tested. The reference strategy consists in injecting only water at a constant ux rate. Strategy A consists in injecting water at a higher ux rate after half the time of the experiment. The second injection strategy is found commonly in the literature. Water and substrate are injected continuously at a xed ux rate and after increasing the water ow to diverge the water ow. Strategy C consists in inverting the direction of the high ux to study the eect on the oil recovery in comparison with Strategy B. For Strategy E, rst substrate is injected at a xed ux rate and double concentration. The system is then closed in order to let the bacteria consume the suspended substrate. Afterwards, the water ux is reactivated at the high rate. We observe that Strategy B and Strategy E use the same amount of substrate. Fig. 10 shows the three dierent oil predictions for the dierent injection strategies. From the previous plot we observe that Strategy E predicts the largest oil recovery. Changing the ow direction at half of the injection (Strategy C) does not result in an improvement on the oil recovery in comparison to Strategy B. Strategies C and D give similar predictions of oil recovery. Fig. 11 shows the simulation results for pressure, substrate, and biomass along the core for Strategy D. The cells with water saturation above 0.5 are shown in subplot (a). We observe that after three hours of injection the water has displaced most of the oil in the thief zone. In subplot (b) the cells where oil saturation is above 0.5 after 6 hours of injection are shown. We observe that most of the oil recovered is from the thief zone, but there is still signi cant amount of oil in the core. Subplot (c) shows the substrate concentration on the lower half part of the core after 12 hours. As described in Strategy D, the ux direction has changed, so now substrate is being injected on the right side of the core, as shown in the gure. From this subplot we observe that the substrate concentration decreases along the core, where the lower Ratio of Initial to Reduce Permeability Ratio of Initial to Reduce Permeability 100 90 Water (+ Substrate) High water flux Water flooding Vandevivere Thullner Channel Tube 0 2 4 6 8 10 12 14 16 18 20 22 24 Time [h] Fig. 8|Percentage of oil recoveries for the dierent porosity-permeability relationships. Reference Strategy A Strategy B W W W W+S W t t t Strategy C Strategy D Strategy E W+S W+2S W W+S W+S W W t t t Fig. 9|Reference and injection strategies. values are located in the thief zone. Given the initial homogeneous volume fraction of bio lm, this has increased more than two orders of magnitude in comparison with the initial volume fraction, showing a greater value on the left side. This result shows that changing the substrate ux direction will not lead to a symmetric bio lm formation. Oil recovery [%] 100 90 Water (+ Substrate) High water flux Reference Strategy A Strategy B Strategy C Strategy D Strategy E 0 2 4 6 8 10 12 14 16 18 20 22 24 Time [h] Fig. 10|Comparison of the oil recoveries for the dierent injection strategies over time. Fig. 11|Saturations, substrate, and biomass pro les along the core. Oil recovery [%] Sensitivity Analysis: Numerical Results. Sensitivity analysis is performed for two test cases based on the methodology presented in a previous section. Case I: Mean Permeability with Variability in Seven Parameters. We consider the single-phase core- scale mathematical model in Table 1 and the same boundary and initial conditions as in Model Test. The quantity of interest q is the mean permeability of the core after 100 days of water injection, divided by the initial mean permeability of the core, and we investigate its sensitivity with respect to the seven parameters (i.e., n = 7) with ranges shown in Table 4. The interpolant I (q) is computed for p = n + 0; : : : ; n + 6. The relative interpolation p;n error is estimated as the L norm of the hierarchical surplus at the nest level over the norm of I (q) itself. The 2 p;n interpolant at level p = n + 6 requires 30,241 evaluations of q, and results in an estimated relative error of less than 0.011. Parameter Notation Range Total Sobol Index 6 1 Maximum growth rate [4.5, 5.5]10 s 0.76 6 1 Bacterial death rate K [2.88, 3.52]10 s 0.33 Stress coecient K [1.35, 1.65 ]10 m/Pas 0.017 str Power law constant [2.25, 2.75] 0.0015 Monod-half velocity K [0.82, 1.01] kg/m 0.0017 Biomass density [18, 22] kg/m 0.0012 10 2 Nutrient diusion coecient D [4.5, 5.5]10 m /s 0.0012 Table 4|Total variability contribution of input parameters to percentage of oil extraction. The surrogate function I (q) is subsequently sampled N = 10 times and the multi-wavelet coecients are p;n computed by solving (26) for varying polynomial orders o and resolution levels L. The fraction of the total sampling variance explained by the retained multi-wavelet coecients are shown in Table 5. For this problem, the convergence is faster in the order o of the piecewise polynomials than in the resolution level L of the multi- wavelets, suggesting that q is relatively smooth. Only combinations of multi-wavelets of total polynomial order o and total resolution level L have been included in the bases used to generate the numerical results. The observed accuracy with respect to capturing the total sample variance suggests that this basis truncation is indeed a suitable strategy to reduce the computational cost of solving the linear system (26). 0 1 2 3 4 0 - 0.883 0.931 0.974 0.994 1 0.695 0.901 0.949 0.980 0.997 2 0.883 0.917 0.952 0.982 0.999 Table 5|Fraction of sample variance in test Case I represented by multi-wavelet expansion of polynomial order o and resolution level L. The total contribution of each material parameter, isolated, and in combination with the others, are shown in Table 4. Only maximum growth rate and bacterial death rate exhibit signi cant eect on the variability in q for the parameter ranges considered here. Case II: Two-Phase Flow with Three Variable Parameters. We now consider the case of two-phase ow (oil and water). The injected water velocity changes after 8 hours from v = 1 m/day to a higher velocity of v = 50 m/day. Along the length of the core, we discretize with 30 elements, while on the transversal section we discretize with 9 times 9 elements. Due to the computational complexity, only three parameters are considered in the sensitivity study: initial volume fraction, position, and length of the bio lm, as presented in Table 6. The remaining parameters, boundary, and initial conditions are the same as described at the beginning of this section. The quantity of interest is the percentage of oil extraction compared to the initial oil in the core (0-100%). We expect nonsmooth parameter dependence due to sharp changes in velocity. Parameter Notation Range Total Sobol Index 4 3 Volume fraction of the bio lm [2.510 , 2.510 ] 0:86 Length of the bio lm L [2, 20] cm 0:18 Position of the bio lm centre along the core X [11, 19] cm 0:05 Table 6|Total contribution of each initial parameter. A sparse grid interpolant on 7 levels yields an estimated error of 0.014 between the two nest levels of resolution, where the relative error has again been estimated as the L norm of the hierarchical surplus divided by the norm 14 6 of the solution itself. A set of N = 10 samples are drawn from the interpolant surrogate model to t a multi- wavelet model through ordinary least squares. As shown in Table 7, essentially all variance is captured by a multi-wavelet expansion in total level 2 and with piecewise quadratic basis functions. The total Sobol indices for this multi-wavelet representation are shown in Table 6. For the parameter ranges investigated, the variability in oil extraction is dominated by the initial volume fraction of the bio lm. The length and position of the bio lm should however not be entirely ignored. 0 1 2 3 0 - 0.887 0.964 0.991 1 0.662 0.901 0.970 0.991 2 0.872 0.905 0.970 0.991 3 0.950 0.907 0.971 0.993 Table 7|Fraction of sample variance in test Case II represented by multi-wavelet expansion of polynomial order o and resolution level L. Conclusions In this work we discuss a core-scale mathematical model for single- and two-phase ow including the transport of substrate and changes on the permeability due to formation of biomass. The single-phase laboratory experiment shows that the substrate input changes the plug-potential. The single-phase mathematical model captured the observed response of permeability to changes in the substrate ux. For the two-phase mathematical model, we investigated the eects on the simulated oil recovery for dierent empirical and upscaled porosity-permeability relations. These results show that the predicted oil recovery could be over- or underestimated depending on the assumed porosity-permeability relation in the mathematical model. Numerical simulations are performed for dif- ferent injection strategies to study the oil recovery. After simulations, injecting substrate and stopping the water ow to let the bacteria consume the substrate and after reactivating the ow at a higher rate results in the largest oil recovery prediction. The sensitivity analysis for the single-phase core-scale model shows that two parameters are responsible for almost all variability in the mean permeability: maximum growth rate and bacterial death rate. Both parameters need to be estimated in the laboratory with sucient accuracy to lead to a reliable estimate of the permeability changes. The sensitivity analysis for the two-phase core-scale model demonstrates less impact on the total variability in oil extraction from the initial position of bio lm as compared to the initial volume fraction and length of the bio lm. Thus, the amount of initial bio lm has a higher impact on the oil recovery in comparison to the initial position of bio lm in the thief zone. Nomenclature a = weighting factor, dimensionless B , B = critical and relative porosity, dimensionless c r 3 3 B , B = injected brine concentrations, m/L , kg/m 1 2 3 3 C , C = injected substrate concentration and substrate concentration, m/L , kg/m i n d = core diameter, L, m D = substrate diusion coecient, L /s E, F , G, V , W , X = integration coecients, dimensionless 2 2 g = gravity, L/t , m/s h = bio lm thickness, dimensionless I = all subsets of parameters including parameter i for global sensitivity analysis, dimensionless I (q) = hierarchical multidimensional interpolant of q on total level p n, dimension dep. on q p;n I = multi-wavelet index set of nonnegative integers MW 2 2 j = substrate ux, m/t L , kg/sm J = Bessel function of order of rst kind, dimensionless 2 2 k, k , k , = rock permeability, initial rock permeability, and bio lm permeability, L , mdarcy [m ] 0 b 2 2 k , k , k = power law, weighted, and Verma-Pruess permeability relationships, L , mdarcy [m ] p h vp 2 2 k , k , k , k = channel, tube, Thullner et al., and Verma-Pruess permeability relationships, L , mdarcy [m ] c t th v k ; k = oil and water relative permeabilities, dimensionless r;o r;w 15 ~ k = bio lm permeability, dimensionless 1 1 K = bacterial death rate, t , s K = stress coecient, L t/m, m/sPa str 3 3 K = Monod half-velocity coecient, m/L , kg/m ` = re nement level, dimensionless L = number of resolution levels, dimensionless L, L = core and bio lm length, L, cm n = number of stochastic dimensions, dimensionless N = set of n-tuples of nonnegative integers, dimensionless o = polynomial order of the multi-wavelet expansion, dimensionless p , p = oil and water pressure, m/Lt , Pa o w P = size of multi-wavelet basis q = quantity of interest for global sensitivity analysis, varying dimension r = core radius, L, cm R = resistance factor, dimensionless 3 3 R = substrate reaction term, m/tL , kg/sm S ; S = saturation of oil and water, dimensionless o w S = total Sobol index for parameter i, dimensionless fig t = time, t, days [hours] T = temperature, T, °C v , v , v = injected water velocity, oil, and water velocity, L/t, m/s i o w w = variable depending on the bio lm thickness, dimensionless X = position of the bio lm centre along the core, L, cm y = general parameter for global sensitivity analysis, dimensionless Y = interpolation nodes in the parameter space Y = Bessel function of order of second kind, dimensionless Y = set of interpolation nodes in the parameter space, dimensionless = tting factor Brooks-Corey relationship, dimensionless = factor Brooks-Corey relationship, dimensionless = tting factor Power law, dimensionless = bio lm water content, dimensionless 1 1 = maximum speci c biomass production rate, t , s , = water and oil viscosity, m/Lt, Pas o w = variable dependent on bio lm permeability and porosity, dimensionless 3 3 ; ; = density of biomass, oil, and water, m/L , kg/m b o w = rock porosity, dimensionless , = volume fraction of bio lm and void space outside the bio lm, dimensionless b f ' = jth orthogonal basis function, dimensionless = piecewise linear interpolation functions, dimensionless = range of independent parameters for global sensitivity analysis, dimensionless Acknowledgments The work of DLM, GB, BFV, KK, PP, and FAR was partially supported by the Research Council of Norway through the projects IMMENS no. 255426, MICAP no. 268390, and CHI no. 255510. ISP was supported by the Research Foundation-Flanders (FWO), Belgium through the Odysseus program (project G0G1316N) and the Akademia grant of Equinor ASA. The authors also appreciate the support from Equinor ASA related to the experimental work reported herein. References [1] Alpkvist, E. & Klapper, I. 2007 A multidimensional multispecies continuum model for heterogeneous bio lm development. Bull. Math. Biol. 69 (2), 765{789. [2] Bao, K., Lie, K.-A., Myner, O. & Liu, M. 2017 Fully implicit simulation of polymer ooding with MRST. Comput. Geosci. 21 (5-6), 1219{1244. [3] Brockmann, D., Rosenwinkel, K.-H. & Morgenroth, E. 2006 Modelling deammoni cation in bio lm systems: Sensitivity and identi ability analysis as a basis for the design of experiments for parameter estimation. Comput. Aided Chem. Eng. 21, 221{226. [4] Carman, P.C. 1937 Fluid ow through granular beds. Trans. Inst. Chem. Eng. 15, 150166. 16 [5] Corey, A.T. 1954 The interrelation between gas and oil relative permeabilities. Prod. Mon. 19, 38{42. [6] Duddu, R., Chopp, D. L. & Moran, B. 2009 A two-dimensional continuum model of bio lm growth incorporating uid ow and shear stress based detachment. Biotechnol. Bioeng. 103 (1), 92{104. [7] Hardy, B., Sarko, A. 1993 Molecular dynamics simulation of cellobiose in water. J. Comput. Chem. 14 (7), 848{857. [8] Hommel, J., Coltman, E. & Class, H. 2018 Porosity{permeability relations for evolving pore space: A review with a focus on (bio-)geochemically altered porous media. Transp. Porous Med. 124 (2), 589{629. [9] Ives, K. & Pienvichitr, V. 1965 Kinetics of the ltration of dilute suspensions. Chem. Eng. Sci. 20 (11), 965{973. [10] Kim, S.B. 2006 Numerical analysis of bacteria transport in saturated porous media. Hydrol. Process. 20 (5), 1177{1186. [11] Klimke, A., & Wohlmuth, B. I. 2005 Algorithm 847: Spinterp: piecewise multilinear hierarchical sparse grid interpolation in Matlab. ACM Trans. Math. Softw. 31 561{579. [12] Kundu, S., Ghose, T.K. & Mukhopadhyay, S.N. 1983 Bioconversion of cellulose into ethanol by Clostrid- ium thermocellum { product inhibition. Biotechnol. Bioeng. 25 (4), 1109{1126. [13] Lacerda, E.C.D.S., Priimenko, V.I. & Pires, A.P. 2012 Microbial EOR: A Quantitative Prediction of Recovery Factor. Society of Petroleum Engineers. [14] Landa-Marban, D., Liu, N., Pop, I. S., Kumar, K., Pettersson, P., Bdtker, G., Skauge, T. & Radu, F. A. 2019 A pore-scale model for permeable bio lm: Numerical simulations and laboratory experiments. Transp. Porous Med. 127 (3), 643{660. [15] Landa-Marban, D., Bdtker, G., Kumar, K., Pop, I. S., & Radu, F. A. An upscaled model for permeable bio lm in a thin channel and tube. Submitted [16] Landa-Marban, D., Pop, I.S., Kumar, K. & Radu, F.A. 2019 Numerical Simulation of Bio lm For- mation in a Microchannel. In Numerical Mathematics and Advanced Applications ENUMATH 2017 (ed. Radu, F.A., Kumar, K., Berre, I., Nordbotten, J.M. & Pop, I.S.r), 126, 799{807, Springer International Publishing, Cham. [17] Le Matre, O., Najm, H., Ghanem, R., & Knio, O. 2004 Multi-resolution analysis of Wiener-type uncertainty propagation schemes. J. Comput. Phys. 197, 2 502{531. [18] Li, J., Liu, J., Trefry, M. G., Park, J., Liu, K., Haq, B., Johnston, C. D., & Volk, H. 2011 Interactions of Microbial-Enhanced Oil Recovery Processes. Transp. Porous Med. 87 (1), 77{104. [19] Lie, K.-A. 2019 An Introduction to Reservoir Simulation Using MATLAB/GNU Octave: User Guide to the MATLAB Reservoir Simulation Toolbox (MRST). Cambridge University Press. [20] Linville, J., Rodriguez, M., Mielenz, J. & Cox, C. 2013 Kinetic modeling of batch fermentation for Populus hydrolysate tolerant mutant and wild type strains of Clostridium Thermocellum. Bioresour. Technol. 147, 605{613. [21] Liu, N., Skauge, T., Landa-Marban, D., Hovland, B., Thorbjrnsen, B., Radu, F. A., Vik, B. F., Baumann, T. & Bdtker, G. 2019 Micro uidic study of eects of owrate and nutrient concentration on bio lm accumulation and adhesive strength in a microchannel. J. Ind. Microbiol. Biotechnol. 46 (6), 855{868. [22] Ma, X. & Zabaras, N. 2009 An adaptive hierarchical sparse grid collocation algorithm for the solution of stochastic dierential equations. J. Comput. Phys. 228 3084{3113. [23] Patel, I., Borgohain, S., Kumar, M., Rangarajan, V., Somasundaran, P. & Sen, R. 2015 Recent developments in microbial enhanced oil recovery. Renew. Sustain. Energy Rev. 52 1539{1558. [24] Rahman, S. 2014 A generalized ANOVA dimensional decomposition for dependent probability measures. SIAM/ASA J. Uncertain. Quantif. 2 670{697. [25] Ro, K. S. & Neethling, J. B. 1991 Bio lm density for biological uidized beds. Res. J. Water Pollut. Control Fed. 63, 815{818 17 [26] Schulz, R. & Knabner, P. 2016 Derivation and analysis of an eective model for bio lm growth in evolving porous media. Math. Methods Appl. Sci. 40 (8), 2930{2948. [27] Smolyak, S. 1963 Quadrature and interpolation formulas for tensor products of certain classes of functions. Soviet Mathematics, Doklady 4 240{243. [28] Sobol, I. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math.Comput. Simulat. 55 (1) 2001 271{280. [29] Sudret, B. 2008 Global sensitivity analysis using polynomial chaos expansions. Reliab. Eng. Syst. Safe. 93 (7) 964{979. [30] Suthar, H., Hingurao, K., Desai, A. & Nerurkar, A. 2009 Selective Plugging Strategy-Based Microbial-Enhanced Oil Recovery Using Bacillus licheniformis TT33. J. Microbiol. Biotechnol. 19 (10), 1230{ [31] Thullner, M., Zeyer, J. & Kinzelbach, W. 2002 In uence of microbial growth on hydraulic properties of pore networks. Transp. Porous Med. 49 (1), 99122. [32] van Noorden, T. L., Pop, I. S., Ebigbo, A. & Helmig, R. 2010 An upscaled model for bio lm growth in a thin strip. Water Resour. Res. 46 (6), W06505. [33] Vandevivere, P. 1995 Bacterial clogging of porous media: a new modelling approach. Biofouling 8 (4), [34] Verma, A. & Pruess, K. 1988 Thermohydrological conditions and silica redistribution near high-level nuclear wastes emplaced in saturated geological formations. J. Geophys. Res.: Solid Earth 93 (B ), 1159{1173. [35] Wood, D.A. 2019 Microbial improved and enhanced oil recovery (MIEOR): Review of a set of technologies diversifying their applications. Advances in Geo-EnergyResearch 3 (2), 122{140. [36] Xiu, D. & Karniadakis, G. E. 2002 The Wiener{Askey polynomial chaos for stochastic dierential equa- tions SIAM J. Sci. Comput. 24, 2 619{644. Appendix A | Eective Porosity-Permeability Relations A detailed description of the following two porosity-permeability relationships can be found in [15], where both relationships are derived by homogenization of a pore-scale model. To this aim we let h be the dimensionless thickness of the bio lm layer, the volume fraction of bio lm, the initial porosity, the bio lm porosity, and b 0 w k the bio lm permeability. The thickness of the bio lm h is given as a function of the volume fraction of bio lm b b ~ ~ and the initial porosity for the thin channels as h = = , whereas for the thin tubes h = 1 1 = . b 0 b b 0 b b 0 ~ ~ We use the notation w = 1 h and k = k =k . b b b 0 The eective porosity-permeability relation for a porous medium modeled as a stack of thin channels is given by h i h i ~ ~ 3 W exp () exp(h ) 1 X exp() exp(h ) 1 b b k w ~ ~ = wV +k h ; : : : : : : : : : : : : : : : : : : (A-1) b b k 6 where h i h i ~ ~ ~ ~ ~ ~ ~ w + 2k exp h + exp h + 2w k exp h exp h 4k b b b b w b b b h i V = ; : : : (A-2) ~ ~ 2 exp h + exp h b b ~ ~ k exp (w) w k exp () b b w W = ; : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (A-3) ~ ~ exp h + exp h b b ~ ~ k exp (w) + w k exp () b b w X = ; : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (A-4) ~ ~ exp h + exp h b b where = =k . w b 18 The eective porosity-permeability relation for a porous medium modeled as a stack of thin tubes is given by k w 2 (Y () F J () G wY (w) F + wJ (w) G) t 1 1 1 1 2 2 = w E + + k 1 w : : : : : : : : : : : : (A-5) k 8 where 2w [J () Y (w) J (w) Y ()] + k [J (w) Y (w) + Y (w) J (w)] w 0 0 0 0 b 0 1 0 1 E = 4 [J () Y (w) + Y () J (w)] 0 1 0 1 4k + w [J () Y (w) + Y () J (w)] b 0 1 0 1 ; : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (A-6) 4 [J () Y (w) + Y () J (w)] 0 1 0 1 2k Y (w) + w Y () b 1 w 0 F = ; : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (A-7) 2 [J () Y (w) + Y () J (w)] 0 1 0 1 2k J (w) + w J () b 1 w 0 G = ; : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : (A-8) 2 [J () Y (w) + Y () J (w)] 0 1 0 1 where = i =k and i is the imaginary number. Here, J (z) and Y (z) are the Bessel function of order of w b rst and second kind respectively. David Landa-Marb an is a post-doctor at the Norwegian Research Centre at the Energy Department. He is interested in multiscale modeling for bio lms. Landa-Marb an holds a PhD degree in applied mathematics from the University of Bergen. Gunhild Bdtker is a senior researcher at the Norwegian Research Centre, NORCE. Bdtker is research director for the group integrated microbiology, chemistry, and physics at the NORCE Energy Department. She is inter- ested in reservoir microbiology, MEOR, reservoir souring, and bio lm injectivity. Bdtker holds a PhD degree in microbiology from the University of Bergen. Bartek Florczyk Vik is a senior researcher at the Norwegian Research Centre at the Energy Department. He is interested in laboratory experiments for petroleum microbiology. Vik holds a PhD degree in physics from the University of Bergen. Per Pettersson is a senior researcher at the Norwegian Research Centre at the Energy Department. He is interested in hyperbolic problems, uncertainty quanti cation, numerical methods for subsurface CO2 storage, boundary conditions, and time-stability. Pettersson holds a double PhD degree in scienti c computing from the Uppsala University and in computational and mathematical engineering from Stanford University. Iuliu Sorin Pop is a professor at University of Hasselt. He is interested in ow and reactive transport in porous media, bio lm growth, and geothermal energy. Pop holds a PhD degree in mathematics from the Babes-Bolyai University. Kundan Kumar is a senior lecturer of mathematics at Karlstad University, Sweden and holds an associate pro- fessor II position at the University of Bergen, Norway. He is interested in upscaled models and computational tools for coupled multiphysics processes, domain decomposition techniques, and iterative and multi-rate algorithms for coupled ow and geomechanics. Kundan holds a PhD degree (cum laude) in applied mathematics from the Eind- hoven University of Technology, The Netherlands. Florin Adrian Radu is a professor of applied mathematics at University of Bergen. He is interested in mathe- matical modeling, analysis, and numerical simulation of multiphase ow and multicomponent reactive transport in porous media. Radu holds a PhD degree in mathematics and a habilitation degree from the University of Erlangen-Nuremberg, Germany.
Physics – arXiv (Cornell University)
Published: Dec 18, 2019
You can share this free article with as many people as you like with the url below! We hope you enjoy this feature!
Read and print from thousands of top scholarly journals.
Already have an account? Log in
Bookmark this article. You can see your Bookmarks on your DeepDyve Library.
To save an article, log in first, or sign up for a DeepDyve account if you don’t already have one.
Copy and paste the desired citation format or use the link below to download a file formatted for EndNote
Access the full text.
Sign up today, get DeepDyve free for 14 days.
All DeepDyve websites use cookies to improve your online experience. They were placed on your computer when you launched this website. You can change your cookie settings through your browser.