Reduced Model for Properties of Multiscale Porous Media with Changing Geometry

: In this paper, we consider an important problem for modeling complex coupled phenom-ena in porous media at multiple scales. In particular, we consider ﬂow and transport in the void space between the pores when the pore space is altered by new solid obstructions formed by microbial growth or reactive transport, and we are mostly interested in pore-coating and pore-ﬁlling type obstructions, observed in applications to bioﬁlm in porous media and hydrate crystal formation, respectively. We consider the impact of these obstructions on the macroscopic properties of the porous medium, such as porosity, permeability and tortuosity, for which we build an experimental probability distribution with reduced models, which involves three steps: (1) generation of independent realizations of obstructions, followed by, (2) ﬂow and transport simulations at pore-scale, and (3) upscaling. For the ﬁrst step, we consider three approaches: (1A) direct numerical simulations (DNS) of the PDE model of the actual physical process called BN which forms the obstructions, and two non-DNS methods, which we call (1B) CLPS and (1C) LP. LP is a lattice Ising-type model, and CLPS is a constrained version of an Allen–Cahn model for phase separation with a localization term. Both LP and CLPS are model approximations of BN, and they seek local minima of some nonconvex energy functional, which provide plausible realizations of the obstructed geometry and are tuned heuristically to deliver either pore-coating or pore-ﬁlling obstructions. Our methods work with rock-void geometries obtained by imaging, but bypass the need for imaging in real-time, are fairly inexpensive, and can be tailored to other applications. The reduced models LP and CLPS are less computationally expensive than DNS, and can be tuned to the desired ﬁdelity of the probability distributions of upscaled quantities. and many-pore geometries. The mean and standard deviation of K are presented for several selected vales of the relative volume V o = | D o | | D v | , as shown. For the many-pore geometry and V o = 0.1, the method LP generates 20 out of I = 100 impermeable geometries. For V o = 0.15, most are impermeable. These facts result in a large standard deviation, and these results are not included.


Introduction
In this paper, we develop reduced computational models motivated by a particular class of multiscale applications in porous media. Porous media are ubiquitous in nature and technology, and flow and transport in porous media has important applications in groundwater management, oil and gas reservoirs, modeling human tissue including lung and bone, and in other materials including wood and paper products.
Porous media are made of solid grains surrounded by the fluid domain D f in pore space; see Figure 1 for illustration. The modeling of processes in porous media proceeded only at a Darcy scale of [cm] and above up until the year 2000. At Darcy scale, only the Darcy properties, such as porosity φ (average volume fraction of the voids) and permeability K (average proportionality coefficient in Darcy's law for momentum equation), were studied and experimentally obtained, and the details of D f were not recognized in practice, even if they were recognized as important [1]. In particular, there are well known correlations for K = K(φ) derived from idealizing the pore-scale geometry of D f as a bundle of tubes or space between, e.g., periodically distributed spherical grains. In a more general context, the wide variety of porous media and the uncertainty of the actual values of φ and K in heterogeneous reservoirs led to research on the stochastic modeling of φ(x, ω) and K(x, ω), such as in [2][3][4][5], where ω is the random variable representing the uncertainty. Since about the year 2000, advances in imaging led to an explosion of micro-images of D f and of its surrounding rock D r -see, e.g., [10][11][12] and the Digital Rocks portal [13]. See also efforts on pore network modeling [10,14,15], and [16] and references therein. Each D f imaged by µ-CT is segmented and digitized as a collection of voxels. The voxel geometry representing D f or its Representative Elementary Volume (REV) subset can be used as a computational domain for direct numerical simulations (DNS) of flow and transport, followed by upscaling of Stokes to Darcy's law, to deliver φ and K for this particular D f . This process is now fairly standard, and its many variants have been widely reported in the literature, with many interesting and important algorithmic and application questions resolved and many more raised. It is now recognized that modeling at the pore-scale and upscaling of the properties of porous media are important for the qualitative understanding and quantitative modeling of the processes at large scales. However, questions remain on how to make use of K found for a particular D f , since the actual geometry of each D f is completely random. One might, therefore, infer that, for a particular collection of samples or corresponding REVs, φ = φ(ω) and K = K(ω). For perspectives on the randomness of porous media, see [1,17]. See also our work on analysis and algorithms for anisotropic non-Darcy flow to upscaling from pore-scale to Darcy scale in [7,8,[18][19][20].
However, much less is known about φ and K in the important case when the pore geometry changes in time due to chemical reactions or similar processes; see Figure 2 (bottom). Here, the recent work includes [9,21,22], which emphasize the multifaceted challenges associated with pore-scale modeling and upscaling when the geometry changes. In principle, one could attempt to image the evolving porous media in time t, observe D f (t), find the corresponding φ(t), and K(t), and record information about some quantity ξ associated with the process creating the solid obstructions. However, obtaining these data are expensive in terms of both material and person-time efforts, as well as being computationally intensive. It might also be technologically infeasible: in the case of biomass, x-ray imaging requires stopping the reproduction, and if done continuously, the process alters their ability to grow. Therefore, a typical workflow involves taking images prior to and after the process is complete. In particular, Figure 2 shows pore-scale after more than 50% of D f has been plugged up. In addition, the values of φ = φ(ξ; t; ω), and K = K(ξ; t; ω) are further burdened with more uncertainty. These challenges motivate the alternative methods we explore in this paper. Aim of this paper. Our goal is to propose a practical method to approximate the connection between φ and K, specific to a particular process which obstructs the pore space, while allowing and taking advantage of the inherent randomness of the problem. We focus on the case when the pore space is being obstructed by a new solid phase, leading to the change in flow paths, clogging of the pores and the decrease in φ and K. In particular, we consider microbial growth leading to the plugging of porous media as well as the process of crystal formation and phase transitions. These two processes create obstructions of the pore-coating and pore-filling type. Other important scenarios not explored in this paper include the dissolution of the rock grains, which leads to the increase in φ(t) and K(t), or to the rearrangement and local changes in φ and K due to mechanical deformation.
Our approach to connect φ to K includes three steps: (1) generation of independent realizations of obstructions in D f , followed by (2) flow and transport simulations at porescale, and (3) upscaling.
Step (1), the generation of new solid obstacles appearing in some original unobstructed D f , can be done in several ways, which exploit computations and randomness in many different ways. Steps (2)(3) to calculate φ and K for the given D f or its REV subset are fairly standard, as we explain in Section 2. In all cases, the variable ξ responsible for the solid appearance is correlated with φ.
First, we consider (1A) Direct Numerical Simulations (DNS) of the actual physical process underlying the changing D f (t) for the microbial growth application and pore-coating obstructions. Each instance i = 1, . . . I of the (Monte Carlo ensemble of) DNS requires an initial condition, which is assumed to be random, and leads to an instance of D  f (t) gives φ (i) (t) and K (i) (t), which we collect in an experimental probability distribution for K(φ; ω) built from φ (i) (t), K (i) (t) i . Unfortunately, DNS are extremely complex and depend on a multitude of parameters which may or may not be possible to calibrate. Rephrasing, DNS give very specific answers to what may be a very non-specific question.
Because of this complexity, we propose two new computational evolution algorithms (1B, 1C) which are surrogates of DNS of lesser complexity and are designed to generate plausible independent realizations D (i) f for a given fixed φ. These eventually give an experimental probability distribution for K(φ; ω) from the collection K (i) (φ) i . The two new models produce results similar to the known experimental and imaging results as well as to DNS. Both models share a design: they generate D f to a particular design of interest. We find these local minima in one of two ways, which we call (1B) Constrained Localized Phase Separattion (CLPS) and (1C) LP.
(1B) is the Constrained Localized Phase Separation (CLPS), which seeks stationary solutions to a partial differential equation, a continuum gradient flow model, and is our proposed modification of the well-known phase separation Allen-Cahn model. In turn, (1C) is the (lattice at pore-scale (LP) model, a discrete model which finds the states of a Metropolis algorithm adapted to pore-scale. Both CLPS and LP are initiated with a random initial state and calibrated to the desired qualitative behavior, consistent with the images from the literature. When run multiple times, we find as many independent realizations I as desired.
Broader context. First, this work is an alternative to traditional multiscale models which use a fixed pore-scale geometry [7,8,16,20,23] to calculate the corresponding Darcy properties. Instead, we consider ensembles of geometries which change to reactive transport. Our entire process is equivalent to creating off-line reduced models of certain quantities of interest, and is a computationally attractive alternative to DNS, which itself is an in silico alternative to imaging and experiments.
This work therefore fits in the area of reduced models [24] as well as in the emerging area of multifidelity approaches [25]. These create a hierarchy of models which approximate the same output but with a different fidelity or modeling accuracy. Our paper addresses the specific challenges of multifidelity and reduced models in the area of pore-scale modeling and upscaling. Specifically, the outputs of our models are pore geometries produced by the DNS, CLPS, and LP models (high-, intermediate-, and low-fidelity models, respectively). While we do not define any metric to compare these outputs, i.e., the geometries, directly, we do so indirectly after upscaling when we compare the Darcy-scale quantities.
We find that, overall, the reduced models CLPS and LP produce geometries which are visually similar to those generated by DNS. Furthermore, the Darcy scale quantities calculated by upscaling are naturally also very similar between DNS and the reduced models. The LP method is very fast but is also somewhat sensitive to the grid resolution and may produce rugged structures at a low obstruction volume. The CLPS method is a bit harder to calibrate and more computationally involved but produces results closer to the DNS. All methods produce upscaled quantities which agree well with those known from the literature for pore coating (mild decrease in K) and pore filling scenarios (dramatic decrease in K).
For upscaling, we draw upon our prior work in collaboration with physical scientists [6,18]. In particular, we set up flow experiments with flow rates in the linear laminar regime, as established in [1,7], and we compare the resulting permeability reduction to those cited in [26,27].
We see possibilities to expand on this work in the future and, in particular, the need to define precise metrics of comparison between the geometries with some form of shape analysis, as well as to connect these results to the data from imaging. Finally, the calibration of the LP and CLPS models was done heuristically, but we envision that it could be guided by some of the emerging data science approaches.
Overview. The outline of this paper is as follows. In Section 2, we define the notation and motivate our interest in the obstructions of two particular types, called "pore filling" and "pore coating". These are important in two separate fields, respectively, hydrate crystal and biofilm modeling. We also explain how to post-process our geometry with steps (2) and (3). In Section 3, we present our (1A) DNS model for pore-coating associated with microbial growth. In Section 4, we present the (1B) CLPS model, and in Section 5, we present (1C) the LP model. In Section 6, we present our results and discuss the complexity of the three methods and sensitivity to the parameters. We close with a Summary and an Appendix with auxiliary tools and details.

Motivation and Notation on Modeling Obstructions
In this section, we explain the specifics of porous media and define the notation. A physical process taking place in porous media can be modeled at Darcy scale and at pore scale. The Darcy-scale models rely primarily on empirical data obtained for cores extracted from soil and rock, and the pore-scale models rely primarily on imaging and micro-fluidics experiments. In this paper, we focus on fluid flow of water, with viscosity µ, and reactive transport with phase transitions.
Numerical simulation of flow and transport aids in the understanding of and helps to optimize human-engineered systems involving porous media, including oil and gas recovery, remediation of contaminated aquifers, and biomedical applications in human body. Traditional practical models of porous media such as aquifers or hydrocarbon reservoirs [1,28] are posed at the Darcy scale.
In principle, one can also simulate various processes at the pore-scale, which delivers more precision. However, the complexity of the transient flow and transport computations is enormous; see [6] for various simulations of flow alone in realistic 3D pore geometries of synthetic porous media and sandstone. See also [29] for the repository of 3d porescale images.
In this paper, we propose reduced models for the multiscale connection between the pore scale and Darcy scale. We consider a porous domain ) at which one recognizes the pore walls, e.g., the interface Γ = ∂D f ∩ ∂D r between the pore (void) space and the solid matrix. We have that where The Darcy scale is the scale of x = O(1[cm]) or above, at which the study of flow and other properties does not recognize D r and D f but rather relies on average properties found experimentally in a laboratory using core samples. In particular, porosity φ (volume fraction) is of interest, and permeability K, the coefficient of proportionality in Darcy flow between volumetric flow and pressure gradient. The coefficients φ and K lump the geometrical information about D r and D f . In practice, one uses the probability distributions K(φ; ω) around some known mean; here, ω refers to the uncertainty.
Imaging. In this paper, we work with geometries of D f and D r found from imaging. See Figures 3 and 4 for illustration of a d = 2 slice D from x-ray ct or (micro) µ-CT tomography. In such an image, D is covered by a union of cells d jl , also called voxels. For simplicity, we will assume henceforth that every voxel is a cube of volume h d , with the side length h fixed by the µ-CT equipment used. For simplicity, we also assume that every D is rectangular. In particular, as in Figure    Obstructions and their volume. In many important applications, the pore-scale geometry changes due to various bio-chemical processes, e.g., due to microbial growth, or crystal precipitation. The flow domain D f itself is thus made of the void space D v , which is filled with ambient fluid, and of the obstructions in D o , which are impermeable or partially permeable to the flow and transport. We have that where Γ ov = ∂D o ∩ ∂D v . In this paper, we consider solid obstructions which form due to the particular bio-geochemical reactions such as crystal precipitation, or biofilm growth. We motivate and describe these below, and illustrate these first in Figure 3.
To quantify the presence of the solid obstructions, we define The reference porosity φ 0 is the initial value of porosity before obstructions form, the current porosity φ is the volume ratio of space available to fluid, excluding obstructions to the entire porous domain, and V o is the obstruction volume ratio relative to the original volume of D f .

Motivation and Models for Pore Coating and Pore Filling Obstructions
Our main interest is to consider pore-coating and pore-filling obstructions. These two completely opposite characteristics are motivated by our long-term modeling projects with an impact on environmental and energy resource studies involving computational models of microbial growth [6,30] and methane hydrate modeling [31][32][33] in porous media.

Pore Coating Obstructions: Biofilm
Scientists and engineers use microbes to alter the flow paths in porous media in various applications. The understanding and quantitative modeling of microbial growth in porous media is important, e.g., in selective plugging for the needs of enhanced oil and gas recovery, as well as in carbon sequestration [34][35][36][37]. The microbial growth is promoted by providing nutrients. The microorganisms form a special extracellular polymer substance (EPS) of high density, which is impermeable or almost impermeable to the external flow, and which typically clings to the pore walls. EPS protects the microbes but clogs the pore void space. In particular, a reduction in permeability is significant; see, e.g., [6,[38][39][40].
The experiments and DNS show that the biofilm growth and the Darcy-scale properties depend on various conditions at which the growth occurs, e.g., on flow rate. However, imaging biofilm is intrusive and harmful; thus, it is close to impossible to predict K(φ; ω) for the growth in a variety of different conditions. We present a DNS model of biofilm-nutrient (BN) dynamics at pore scale in Section 3.

Pore Filling Obstructions: Hydrate Crystals
Methane hydrate is a naturally occurring and highly concentrated form of methane: a crystal lattice comprised of frozen water cages trapping gas molecules. It holds significant quantities of carbon in the global system but forms in areas of low temperature and high pressure and is abundant in permafrost and marine sediments. It is important in energy and environmental studies as an energy source and potential climate accelerant, and some modeling and experiments have been forthcoming [27,31,33,[41][42][43]. In particular, it is known that, at low V o , the hydrate grows in the pore centers and causes significant changes to the permeability and sediment strength; see, e.g., the experiments in [44].
However, observations and experiments on the impact of hydrate on permeability are challenging due to the requirements of low temperature and high pressure needed for their stability [41,[45][46][47]. Thus, it is difficult or very expensive to obtain K(φ; ω) experimentally or to image D v for hydrates.
In turn, DNS for phase field models for hydrate are very complex and dependent on a multitude of parameters, as shown in [48,49]. We have forthcoming work on the pore-scale models of pore-filling crystal formation, but providing DNS here is out of scope.

Obtaining Darcy Permeability K with Computations
K depends on the connectivity within D f in a nontrivial way. Some experimental algebraic correlations for K(φ) assume a particular simple form of D f . For example, the Carman-Kozeny relationship assumes that D f is a bundle of capillaries and postulates K(φ) ∼ φ α , where α ≥ 2; see [1,28].
For more general D f , the workflow to obtain φ, K from fluid flow simulations in D f is called upscaling or numerical homogenization. For laminar flow, the permeability K is found by (2) finding the solutions (u, p) to a viscous flow model at the pore-scale, and (3) upscaling (u, p).
One additional consideration in pore-scale simulations is whether D used in Algorithm [Pore2K] (Algorithm 1) is a Representative Elementary Volume (REV) large enough for it to make sense to consider its Darcy scale properties as representative. This important issue was considered in, e.g., [20]. For example, consider the illustration in Figure 5, in which we show a realistic image from µ-CT as well as an idealized single-pore geometry. With studies in [50], we find that the one-pore domain is sufficient for the study of quantities of interest related to the flow, but may not be adequate for transport-related quantities.
(a) single-pore 50 × 50 (b) many-pore 152 × 114  Given D v we calculate φ from (3). Next, we solve the velocity u(x) = (u 1 (x), u 2 (x)) and pressure p(x) satisfying the Stokes system in D v with the no-slip boundary condition on the pore walls, inflow condition on the portion Γ in of ∂D, and a natural outflow condition on the outflow part Γ out where n is the unit outward normal to Γ out . In simulations, we typically choose Dirichlet inflow velocity u in of parabolic shape with averageū D .
To approximate the solutions to (4), we use an MAC scheme extended to Brinkman flow, which we implemented in MATLAB as described in [30], or the computational environment HybGe-Flow3D [16,51].
The approximations to (u, p) are averaged over D to U and P, respectively, and fit to the Darcy model In particular, we find U by averaging over D f From these, the permeability tensor K is calculated. See, e.g., [7,19,20] for a detailed description of this workflow, including setting up multiple flow scenarios to provide sufficient information to calculate the full tensor K. In this paper, we report only on one scalar value of K, the component corresponding to the flow from left to right. Example 1. We simulate flow with (4) in a "many-pore" geometry D = (0, 1) 2 [mm 2 ]. We consider the case with and without obstructions, and consider D o of both pore-filling and porecoating type. We set the flow to go from left to right, with the average u in = 3.6 × 10 −3 [mm/h] of u in through Γ in on the left, and Γ out on the right; see Figure 6. The flow rate is in the linear laminar regime and is similar to that used in [6]. The parameters used are in Table 1.

Flow: from Pore-Scale to Darcy Scale with Obstruction Formation
When obstructions appear in D o ⊂ D f , they clog the pore space and significantly alter the flow and transport paths in D v , and the quantities φ and K decrease from φ 0 and K 0 , respectively, see illustration in Example 1 and Figure 6.
To quantify this change, one can run physical flow experiments or imaging in parallel to the obstruction formation to obtain φ(t) and K(t). However, to our knowledge, such a set-up is extremely challenging, since it requires stopping the reactions somehow at fixed time intervals to run the flow experiment. In consequence, it involves large uncertainty. Therefore the sampling is frequently limited to only few time steps: the beginning and the end of experiment; see, e.g., [6,41] for the study of (φ, K) for hydrate crystal formation and of biofilm clogging, respectively. However, such imaging efforts in real time are very costly and time consuming and pose additional technical challenges, including how to stop the process for the purposes of imaging. For the applications of interest, including both the biofilm and hydrate applications, as detailed in Section 2.2, it is close to impossible to image D o (t) at a fine resolution of t.
Instead, an unobtrusive way to understand how the pore-scale geometry D v is changing with time is to use Direct Numerical Simulations (DNS) of the underlying bio-chemical process at the pore-scale; see, e.g., the models in [52] or [6,30]. DNS are very appealing because they seem to closely represent the physical reality. However, they involve a multitude parameters with large uncertainty, including the initial conditions. Unfortunately, their complexity is also enormous: they require 3D simulations with O (10 M) cells and many time-steps, and thus cannot be used easily in multiscale scenarios involving simulations of complex coupled systems.
The general strategy to overcome this complexity is to consider reduced or surrogate models [24]; see also our pore-scale-related work in [16]. Reduced models exploit data collected or produced off-line in a library of simulation results to provide answers to queries on some desired quantities of interest. For the case considered in this paper, such a library can provide answers to a query on, say, sampling K from φ for a given φ 0 , K 0 and a given type of process, e.g., microbial growth or hydrate growth. More generally, answers can be provided for pore-coating or pore filling scenarios, with an associated measure of uncertainty. In the current era of data science, these models or queries can be replaced by neural networks.
The methods proposed in this paper are the third way (iii) to create experimental distributions of K(φ; ω). In an abstract setting, they deliver surrogate or reduced models of the true DNS-derived distributions, and they do so very efficiently. While our simulations are not based on first principles, the geometries we generate are similar to those found by µ-CT imaging.

Transport from Pore-Scale to Darcy Scale
The Darcy coefficients φ and K are important for modeling the flow. For transport, other important quantities of interest at Darcy scale are the tortuosity T and breakthrough curves relevant for transport models.
First, we recall the transport model at pore-scale. Once the velocity u(x) is known from (4), we can simulate the advective transport, the first-order PDE The inlet value c D > 0 is constant, and in our experiments, we set c init ≡ 0.
To understand the impact of geometry of pore-scale D v on the transport at the Darcy scale, we can upscale the transport solution c by averaging over D to obtain but a lot of information is lost in this step. As is widely known, it is useful in this context to consider the notion of tortuosity discussed in detail, see, e.g., [1,53] for various definitions and uses. In particular, the tortuosity helps to distinguish channel-like pores from those of more intricate shapes. For the left-to-right flow considered here, we calculate hydraulic tortuosity T = l eh L , where L is the "straight-line" distance in D from left to right, and l eh is the effective path length taken by the fluid through D v . The path length l s eh for each streamline (s) is easy to find, and l eh is found by averaging over all streamlines. In addition, we recall the breakthrough curves which assess the average time it takes to travel through D To compare the breakthrough curves, we also define This α-breakthrough time T (α) is the time when the pore is at α * 100 percent of the total possible outflux.
The tortuosity and breakthrough curves provide information on the transport that is partly lost in upscaling. In particular, tortuosity is used to identify dispersion in the porous medium.
To find T and the breakthrough curves, we apply Lagrangian and Eulerian methods, respectively. In the Eulerian frame, we approximate the solutions to (7a), and calculate C n jl ≈ c(x jl , t n ) with an explicit upwind Finite Volume scheme known as Donor in Cell [54]. The method is stable if the time step τ satisfies the Courant-Friedrichs-Lewy (CFL) condition. In the Lagrangian frame, we first find N s streamlines for the flow by explicit time-stepping of dx dt = u starting from Γ in , and using the approximation to u interpolated on each cell from the solution to (4). Next, we propagate the constant value of c along the streamlines from Γ in to Γ out . In the Eulerian frame, we approximate B(t n ) by summing c n jl over the cells D jl adjacent to Γ out . In the Lagrangian frame, we track the position ξ n s of the front of each streamline s at time t n , and approximate B(t n ) by summing the number of streamlines that have reached the outlet before t n .
Summary. We have identified the important quantities that will be used to study average properties of porous media at Darcy scale from the knowledge of the pore-scale geometry of the void domain D v . In the sections below, we present methods (1A-1B-1C), respectively, which generate the realizations of D (i) v . We denote the upscaled quantities with subscript i. In addition to porosities φ (i) and permeabilities K (i) , we consider T (i) and T (α) i .

DNS for Single-Pore Geometry with Obstacles of the Pore-Coating Type
In this section, we generate multiple realizations of obstructions using the DNS of the process of biomass-nutrient (BN) dynamics. This process creates pore-coating obstructions, as noted in [6]. We consider a fixed D = (0, L x ) × (0, L y ) and a fixed D f . Our DNS produces the domains (D The model we use was first proposed in [6], recently refined in [30], and is further modified to reinforce the pore coating behavior.

Biomass-Nutrient Model
Let (B, N) represent biomass and nutrient concentrations, respectively, chosen in convenient non-dimensional units. Let B * = 1 denote the maximum possible biomass density, and B * = 0.9B * the threshold, such that the region x : B * ≤ B(x, t) ≤ B * represents the mature biofilm phase, which contains microbes as well as the biomass-produced extracellular polymer substance (EPS). This substance is beneficial to the microbes as it provides a safe environment for growth. For the purposes of this pape,r we therefore consider that obstructions are the same as the quite dense phase To model biomass growth and formation, we consider the model simplified from that in [6,30] Here, the nearly singular diffusivity is parametrized by γ = 2, motility d B,0 and B * = 1.01B * , and given by The Lagrange multiplier Λ(B) = ∂I (−∞,B * ] (B) enforces the constraint B(x, t) ≤ B * . (Here, ∂I S denotes the subgradient of the indicator function of a set S). A numerical solution of (12) is challenging due to the nearly singular behavior of d B , and the presence of ∂I (−∞,B * ] (B), which makes (12) a nonlinear parabolic variational inequality.
We use the linear diffusivity model for N where the nutrient diffusivity in D v d N,v = 6 [mm 2 /h], which is close to the molecular diffusivity of water, and d N, We also have the biomass growth, and the nutrient consumption rates are given as with the utilization rate κ in O([1/h]), the specific substrate uptake rate κ N ≈ 0.5, and the Monod half N 0 . In this paper, we focus on the nutrient rich case and assume that N D is large enough. The role of a(x) in (12j) is to "promote" the pore-coating behavior of biofilm. This is a well known feature of microbial growth: the microbes stay close to the wall, which provides protection to the community from detachment and other mechanical deformation due to large flow rates; they achieve this by excreting an adhesive substance which supports this mechanism. However, this pore-coating behavior is difficult to calibrate quantitatively, and we choose to do this by heuristically decreasing the reaction rates from Γ according to the following model, depending on a parameter A, as in (12j), where a o is the solution to the Poisson equation −∇ 2 a o = A with homogeneous Dirichlet boundary conditions on Γ and Neumann boundary conditions on ∂D ∩ ∂D f = Γ in ∪ Γ out ∪ Γ wall . For illustration, we provide the plot of a(x) in Figure 7. With A = 0, a ≡ 1, and the rate r B (B, N, x) = κB(x, t)m(N). However, with A = 1, we see that as the distance from the wall γ increases, a ↓ 1.
Finally, we close the system with initial conditions (12e) and (12f), and we choose D b to include only the voxels d jl immediately adjacent to the wall This choice is consistent with the usual assumption that the porous domain is inoculated with microbes which settle down next to the wall Γ. Figure 7. Contours of the attraction coefficient a(x) defined by (12j) in single-pore geometry and many-pore geometry.

Results of the DNS Model for (B, N) Dynamics
We start with a single-pore geometry D = (0, L x ) × (0, L y )[mm 2 ], in which we let B(x, t) evolve from initial state B(x, 0).
We use parameters listed in Table 2 Table 2. Parameters for the DNS model and simulations of (B, N) described in Section 3.1. Simulations cases (a,b,c,d) are presented.

Parameter
Description Value/Units Ref.

D b
Localization of initial biomass random The biomass in each case evolves first, and the biofilm phase becomes mature and impermeable when B * ≤ B(x, t) ≤ B * . Next, this phase continues growing by interface creating new obstructions D o (t), whose volume increases over time. The growth pattern depends on the initial state. With A = 0 and small V b , we see that the biofilm tends to grow spherically without strictly adhering to the walls. However, when we increase V b to cover Γ evenly, the biofilm sticks together but grows gradually away from walls.
On the other hand, with A = 1, the growth away from the walls is suppressed to promote growth near the walls. We see a dramatic difference as compared to the case when

Reduced CLPS Model: Constrained Localized Phase Separation
In this section, we describe our phase separation model. We work with a function ψ(x, t), x ∈ D f , t > 0 called the phase or order parameter, defined in the pore-scale domain D f and evolving in time t. The function ψ(x, t) is a solution of a certain partial differential equations called phase separation model. We approximate ψ(x, t) with a finite difference scheme on the grid (d jl ) jl covering D f ; see Appendix. We post-process the numerical solution to obtain D o .
The phase separation models in the literature include the Allen-Cahn (AC) and the Cahn-Hilliard (CH) Equation [56], often studied together as gradient flows of the functional The parameter has the physical meaning of the interface width, which arises due the competition between the diffusion and coarsening; it also affects the time scale of the evolution; see, e.g., [56,57]. The potential energy density W(ψ) part of (13) is typically a non-convex function with two stable minima ψ * and ψ * , and an unstable local maximum between these. One common choice is ψ * = −1, which corresponds to the solid phase, and ψ * = 1, which corresponds to the liquid phase, with the energy density given by the "double well" function ∼ (1 − ψ 2 ) 2 . Non-polynomial functions, such as in Ginzburg-Landau models, are also possible [56].
The gradient flow for J(ψ) in (13) is the Allen-Cahn (AC) equation which describes the evolution of ψ(x, t) towards one of many stationary solutions ψ(x, ∞) = lim t→∞ ψ(x, t); these are the local minima of J(ψ). The dynamics of the resulting semilinear parabolic PDE are well known: the stationary solutions feature patterns with aggregation in the regions with some diffuse interface region D separating D * and D * of width proportional to , due to the competition of diffusive and reaction terms associated with −∆ψ and f (ψ) = dW dψ , respectively. Various rigorous analyses are available; see, e.g., [58]. With bounded initial data in In this paper, we exploit the ability of AC model to produce patterns as the main design idea, but we modify the model to fit the interactions with the boundary of D, as described next.

Volume Constrained Phase Separation CLPS Model
We exploit the non-convexity of J(ψ) and the non-uniqueness of the stationary solutions to the AC model as the main design idea. We aim to generate patterns of region D o with obstructions with a certain given prescribed volume V(0), and we aim towards |D o | = V(0). We choose ψ * = 1 to denote cells with obstructions and ψ * = 0 (no obstructions). With these, we obtain that We pose the AC model, the gradient flow of (14) in where, for ψ * = 1, ψ * = 0, we use the double-well potential density W(ψ) and its derivative f (ψ) The solutions to (17) aggregate in the regions Typically, these regions D f , * (∞), D * f (∞) are non-empty and they stabilize as t ↑ ∞, and as the values of ψ(x, t) separate on the path to stationary solutions. However, the volume of obstructions is not preserved. This can be seen by integrating (17) over D f × [0, t] and applying Neumann boundary conditions.
As part of our algorithm, we aim to maintain We impose this equation as a constraint, and follow the general strategy outlined in [59]. We augment J(ψ) by the term −λ(V(t) − V(0)) with a Lagrange multiplier λ(t). The resulting evolution model for J(ψ) − λ(V(t) − V(0)) reads Independently, a similar strategy is considered in [60], framed as a non-local order conserving gradient flow model in which one substitutes λ(t) = 1 (20a), and skips (20d).

Controlling the Location of Obstructions in CPS with Localization in CLPS
Consider now a given pore-scale domain D, with D r and Γ known, as shown in Figure 4, and a simulation of the CPS model (20) in D f . Starting with some initial condition, the model will produce aggregated patterns of obstructions indicated in (18), but with the aggregates in D o somewhat randomly located in D f . In particular, these aggregated domains may be close or far from Γ, and therefore will not adhere to a pattern with qualitative properties consistent with the "pore-coating" or "pore-filling" behavior illustrated in Figure 4.
To rectify this, we modify the CPS model to promote the desired qualitative properties of ψ(x, ∞). In particular, we want to control the shape of D o and its distance 1.7em(D o , Γ) to Γ, as we promote the formation of obstruction with the "pore-coating" or the "pore-filling" behavior. In the pore-coating behavior, the obstructions aggregate adjacent to the solid matrix so that the distance 1.7em(D o , Γ) is small. In the pore-filling behavior, they aggregate away from the solid matrix and 1.7em(D o , Γ) is as large as possible. In other words, in the pore-filling behavior, D * f , is "repelled" from Γ, and in the pore-coating behavior it is "attracted" to Γ.
Such behavior can be realized by an application of some form of physically motivated "taxis" towards or away from Γ. In particular, one could set up the taxis with the gradient of an electrostatic or van der Waals potential to model the actual physical attractive or repulsive forces towards or away from Γ. In particular, such taxis can involve electrochemical interactions such as those included in the work on Poisson-Nernst-Planck models as in [61][62][63]. These complex models require a multitude of physical parameters that might be difficult to measure and require additional computations.

Properties of the Solution to CLPS
The dynamics of (21) are competition of coarsening due to f (·), diffusive action of −∆ψ, and the localization function g δ (x, ψ), which is only active when ψ(x) ∈ (0, 1). The solution to (21) maintains its total volume due to the Lagrange multiplier construction and Neumann boundary conditions. Additionally, the interplay of f (·) and g δ (·) is most interesting.
When θ = 0, we work with the dynamical system ψ t + f (ψ) = 0 to see that at each x, the solutions tend towards one of the two stable equilibria ψ * and ψ * , but may also remain at the unstable local maximum 1 2 The action of −∆ψ as in (17) provides spatial variability and allows the patterns of the sets D f , * and D * f to emerge and coalesce.
More generally, consider θ > 0 and the ODE ψ t + θq(ψ) = 0, with a stable equilibrium ψ * = 1, and an unstable one at ψ * = 0. Adding − f to the right hand side, as in the ODE ψ t + f (ψ) + θq(ψ) = 0, does not change the qualitative nature of this evolution towards the stable equilibrium ψ * = 1. In contrast, if θ < 0 and |θ| is large enough, the stable equilibrium ψ * = 0 is "promoted" in spite of the competition between f (ψ) and θq(ψ). The presence of −∆ψ leads to the formation of patterns.

CLPS Model to Generate Multiple Realizations of Obstructions in Pore-Scale Domain
Consider now simulations of (21) in a domain which is a realistic 2D slice of a porescale domain from µ-CT experiments from [6][7][8], with rock domain D r , shown in Figure 3. For this domain, we have a discretization parameter h = 0.05 relative to the domain size L x = 7.6, L y = 5.7; we provide an outline of the numerical scheme in Appendix A.
We must choose the desired volume of obstructions V(0), and the parameters δ, , and θ, which control the qualitative nature of the domain D o with obstructions. We recall that δ controls the width of the interface close to Γ. With θ > 0, we promote pore coating, and with θ < 0, we promote pore filling. The choice of and δ is guided by heuristics and intuition, so that D o has the desired characteristics when assessed visually. In the future, we plan to guide the choice of parameters by a systematic approach guided by data science. Typically, we choose = O(0.1diam(D f )) and δ = O(h). See Table 3 for a summary of all the parameters. Table 3. Parameters for CLPS model from Section 4.3 in single-pore and many-pore geometry. The model uses nondimensional rather than physical units.

Parameter
Description Value/Units Ref.
Model parameters (fixed) The Algorithm CLPS (Algorithm 2) has (A) a pre-processing stage in which we choose parameters, (B) simulation of ψ(x, t), and (C) post-processing step, in which we use ψ(x, T ∞ ) to determine D o . . Choose some random ψ init (x) and run the transient simulation of (21) until t = T ∞ , when the system seems to reach a stationary solution. In practice, we run the simulations until reaching large T ∞ , such that for t > T ∞ , the solution ψ(x, t) does not change significantly.

CLPS [C]. Post-process ψ(x, T ∞ ) and set
Now, as seen in (18), the stationary solution ψ(x, T ∞ ) is close to the characteristic function of D * f (t) ≈ D o (t). Thus, we assign the cell d jl to D o when ψ jl = ψ * = 1, and to D v whenever ψ jl = ψ * = 0. However, when ψ jl ∈ (ψ * , ψ * ), we must make a choice so as to honor (19) with the fixed V(0). We can do so by assigning d jl to D o if ψ jl ≥ψ, and to D v , otherwise, whereψ ∈ (0, 1) is some threshold value chosen depending on V(0): this strategy supports the interpretation that ψ jl represents the amount of obstruction material contained in cell d jl . Alternatively, we assign a certain number of cells to D o with the highest value of ψ jl selected randomly; this strategy is easier to implement than the former.
We can repeat the algorithm CLPS as many times as desired. Since we aim to obtain a large number of realizations, we consider a collection of "white noise" random initial conditions ψ (i) init (x) which satisfy (19), and which are generated by a uniform random generator on the voxel grid for D f . We seek a collection of stationary solutions ψ (i) (x, T ∞ ), where i = 1, . . ., and the corresponding post-processed obstruction region D Table 3. and the many-pore geometry shown in Figure 3, and set V(0) = 0.05. In Figure 10, we present two of the many stationary solutions obtained for each case, where we choose and apply algorithm CLPS. We plot the contours of the stationary solutions ψ (i) (x, T ∞ ) ∈ [0, 1], i = 1, 2, which satisfy V(T ∞ ) = V(0). These are later used in flow simulations with the post-processed results discussed in Section 6.

Pore coating obstructions
Pore filling obstructions

The LP Method from Statistical Mechanics to Generate Obstructions
In this section, we present a method called the LP method to generate geometries with obstructions. Unlike CLPS (21), which is a PDE for the minimizer(s) of the functional J(ψ) in (13) under constraints, the LP method does not use any differential equation, and is motivated by statistical mechanics models of phase transitions from the class of dynamic discrete lattice models from [64,65].
We recall the Markov Chain Monte Carlo (MCMC) method, and, in particular, the Metropolis algorithm [66][67][68], in a setting similar to the Ising model, as in [67]. The Ising model is widely known as a tool for generating arrangements of spins on a lattice, with applications in the study of phase transitions. In this paper, we apply it to the pore-scale geometries represented on a lattice, hence the name LP. We are motivated by the work in [69,70], which in [71], we adapted to multi-component adsorption. We modify the setting from [71] to the present case of pore-coating and pore-filling geometries from µ-CT. We explain the case of d = 2 lattice, but an extension to d = 3 is immediate.
We start with a given digitized rock image D from µ-CT, and a given V(0) = |D f |V o . The practical objective is to assign each voxel d jl of the voxel grid covering D f to either D v or D o . This is equivalent to setting an indicator variable x j,l ∈ {0, 1, 2} to one of three values indicating the void region, rock, and obstructions, as follows Now, each x j,l can be considered as a node of a lattice, and we collect these in the vector X = (x j,l ) j,l . For a collection of realizations, we denote the i'th realization by X (i) , with the corresponding D (i) o , i = 1, . . . I. We start from some initial configuration X (i) init , and we seek its stationary state, X (i) T ∞ , a local minimum of a certain non-convex Hamiltonian functional H(X), which determines the quality of the fit to our particular design of interest. We find these local minima using the Metropolis algorithm.
After the ensembles of geometries (X (i) T ∞ ) i are found, we post-process each one to get D (i) o and to calculate the porosity φ (i) , permeability K (i) and tortuosity T (i) .

LP Model for Minimization of a Hamiltonian Functional Defined on a Lattice
Next, we define the Hamiltonian energy functional H(X), which assigns some energy to each X. This Hamiltonian depends on various parameters which eventually lead to different cell aggregation patterns. In particular, we consider the aggregation patterns of the pore-coating or pore-filling character.

Hamiltonian for Generation of Obstructions
In the design of H(X), we include weights for the interaction at each type of interface between the lattice nodes in D r , D o and D v . The weights allow control over the qualitative character of the configuration of obstructions, such as in the pore-filling or pore-coating patterns.
Consider the node i, j and let Λ i,j = { i + 1, j , i − 1, j , i, j + 1 , i, j − 1 } denote its immediate neighbors. Let the next nearest neighbors be defined as We use periodic boundary conditions; thus, any of the nodes in Λ i,j and Λ 2 i,j are well defined.
Next, we set some weights w ro , w rv , and w vo ∈ R and associate these with the rockobstruction, rock-void, and void-obstruction interfaces, respectively. Now, we define We also let W ≥ be the weight of the next-nearest neighbor interactions, and define the Hamiltonian As usual, we use the Boltzmann distribution e −H(X))/k B Z where Z is the normalizing factor, k B is the Boltzmann constant, and β = 1 k B is the inverse temperature. We follow the literature closely [64,65,67,68].

Metropolis Algorithm for Pore-Coating and Pore-Filling Obstruction Patterns
The Metropolis algorithm for minimization of H(X) starts from some random initial configuration X init to produce the chain, which evolves towards one of the local minima of H(X). The algorithm is inherently stochastic. In each step, the method rearranges the assignment of the cells to D o to find the configurations which are more likely, i.e., have a lower value of energy represented by the Hamiltonian (Algorithm 3). The process is similar to fixing the net magnetic spin in the Ising model and flipping randomly chosen spins so as to find a configuration corresponding to finding one of the local minima of H(X).

Algorithm 3. Algorithm LP [A-B-C]. LP [A].
Generate the initial state X 0 = X init randomly for the given V o . Choose parameters β and W, w ro , w rv , w vo ; see Table 4.
LP [B]. and follow with iteration for updating X t → X t+1 for t = 0, 1, . . . as follows. Let x = X t .

1.
Select the pair of cells x i,j , x r,s uniformly from D f ; 2.
Propose a new state y, similar to x but with the values x i,j of x r,s swapped; 3.
Calculate the change ∆H = H(x) − H(y) to the Hamiltonian. (This step is very inexpensive); 4.
Accept the new state y with probability h = min{1, exp(β∆H)} (a) Select u ∼ U(0, 1); Continue the Markov chain t → t + 1 until satisfied that the so-called "burn-in" process is completed. Stop the iteration at some t = T ∞ .

LP [C].
Postprocess the final state X T ∞ which corresponds to one realization of D o pore geometry with an obstruction. Specifically, the nodes corresponding to x j,l = 2 identify the voxels in D o .
Notice that a different X

Results of LP Model for Single-Pore Geometry
We now illustrate the algorithm for D r as shown in Figure 4. The µ-CT data provide the resolution with M x = M y = 50. We also choose β = 1, W = 1/2. To obtain pore-filling and pore-coating obstructions, we set set the weights, as in Table 4, which avoid or promote the clinging to rock matrix and tend towards clumping.
In Figure 11 we show the progression of the LP model for generating the obstructions for pore-filling. We see that the obstruction cells quickly aggregate. However, it takes many steps for the process to stabilize with a single clump.  We next consider the generation of obstructions that favor adjacency to the rock matrix. The chosen weights also discourage the formation of fingers of obstruction. In Figure 12, we show the progress of the LP model. We can see that the obstruction cells form a thin layer along the rock matrix. The burn-in time for this process is much shorter than in our previous example.

Results and Comparison of DNS and Reduced Models
In this section, we compare the three methods of generating domains with obstructions. For the pore-coating scenario, we compare a method of generating obstructions from (1A) DNS (BN) from Section 3, (1B) CLPS from Section 4 and (1C) LP from Section 5. For the pore-filling scenario, we only consider (1B) and (1C).
The workflow in our experiments is as follows. We start with a porous domain D = D r ∪ D f ∪ Γ with some known D f , i.e., with some known initial porosity φ 0 . This domain may be either a many-pore model, such as from µ-CT imaging, as shown in Figure 3, or a smaller idealized single-pore model, such as in Figure 4. We choose the type of obstructions from either biofilm (pore coating) or hydrate (pore filling). The generation of obstructions is followed by simulations of flow and transport in D v and upscaling to Darcy-scale parameters φ, K, T , T (α) impacted by the presence of obstructions; these steps are as described in Section 2. Of interest is the relative change in the upscaled quantities from their unobstructed values denoted with subscript "0", e.g., φ 0 , K 0 or T 0 , corresponding to D o = ∅.
The first method (1A) produces I independent realizations D (i) o (t) with DNS for the pore-coating scenario, with the associated K (i) (t). The simulations start with initial conditions B init,(i) (x) localized in randomly chosen grid cells close to Γ, with a fixed initial amount given by parameter V b , which indicates the total initial biomass amount. We simulate each case up until the timepoint t at which the desired ratio V o in |D o (t)| = V o (t)|D f | is obtained, or until clogging, whichever occurs first.
The method (1B) CLPS produces I independent realizations D init . We run the evolution until the stationary solutions are found.
The first overall observation is that we find a qualitative agreement between the DNS model BN and the reduced models CLPS and LP both by visual inspection of geometries and through the probability distribution of upscaled properties. There are also some differences, which we discuss in detail. This is of primary importance for our purposes, as we aim towards the creation of reduced models.
We also present additional interesting features for the pore-coating and pore-filling obstructions, which we discuss in detail below. We first focus on flow properties in Section 6.1, and then on transport results in Section 6.2.

Permeability Distributions in the Idealized "Single-Pore" Geometry
We consider an idealized single-pore domain D, as in Figure 4. We generate I = 100 geometries with the DNS model BN, starting with a completely random initial microbial amount situated next to Γ, as shown in Figure 13. We also use this for initial BN D b with a V b = 0.1 identical to the regions produced by LP method. v are shown for the single-pore geometry for different methods in Sections 3-5. They are generally qualitatively similar to each other within the category of pore-coating or pore-filling, but also feature some differences, revealed by visual inspection. The solutions to BN models depend on how closely we enforce the adhesion of biofilm to the boundary Γ, i.e., on parameter A.

Pore-Coating Biofilm-Like Obstructions
We start by emphasizing the difference between (1A) and (1B-C). We recall from (3) that Thus, φ (i) (t) obtained with DNS model (1A) varies (continuously) in time t, but φ (i) are actually constant for each experiment with (1B) and (1C) when using selected volume fractions from (25). This is shown well in Figures 14 and 15, where we compare the scatter plots of K(V o ) with the results of BN vs. CLPS vs. LP model for the pore-coating scenario.
For BN discussed in Figure 14, the shape of geometries, and thus the relationship between K (i) (t) and V (i) o (t), depends on how these geometries were created, and, in particular, on the choice of A. We recall snapshots found with BN and shown in Figure 8, and now study the evolution of K (i) (t) in Figure 14. We see that K(t) decays close to linearly with V o , with more uniform pore-coating, and is promoted with A = 1 and large V b , but we see faster decay of K with A = 0 or large V b . In this context, an important natural question arises. Since all the geometries we create are random, it is not useful to compare them one by one. However, we want to know whether a BN initialized with one of geometries created by CLPS or LP with small V o = 0.1 will predict evolution towards large V o which is consistent, on average, with those produced by the reduced models. We confirmed this, and we saw a very close agreement between the distributions given by BN and CLPS and LP, as seen by comparison of plots in Figures 14 (bottom) and 15. This result is desirable but not unexpected, because the geometries corresponding to more tightly prescribed initial conditions tend to evolve in a more predictable way. Figure 15. Dependence of K on V o when D o are found with reduced models CLPS and LP discussed in Section 6.1.1.

Remark 1.
Furthermore, we want to know whether the trends in the decrease in permeabilities predicted by DNS and plotted in Figure 14 are realistic, and how they compare with the realistic data. For this, we recall the model for the dependence of K on V o based on Carman-Kozeny models for pore-coating and pore-filling scenarios ([27] Figure 4, p.8), given from algebraic expressions for the pore-coating strategies, and for pore-filling. In the range V o ∈ [0, 0.5] considered here, these formulas predict a mild (linear) decrease in K for pore coating, but a dramatic one (convex) for pore filling, in the shown range of V o .
Analysing Figures 14 and 15, we see that the case of strong pore coating corresponding to A = 1 is consistent with the mild linear decrease discussed in Remark 1.
Next, we discuss the histograms and other statistical information on the probability distribution of (K (i) /K 0 ). We compare the results for BN with A = 0 and A = 1 and the different initial biomass domain fraction V b , CLPS and LP. We consider the distributions shown in Figures 16-19 for DNS model BN, and CLPS and LP in Figures 20 and 21, respectively. For each, we present a cartoon of a sample geometry D Studying the sample snapshots of D o in these figures, we see that the obstructions aggregate in distinct "colonies" for DNS with A = 0 and for CLPS, but do not cover the entire rock-fluid interface. This feature is dissimilar to the case of DNS with A = 1, and to the results of the LP model in Figure 21. It seems that the LP model is tuned to promote adhesion to Γ more strongly than CLPS. The difference between CLPS and LP is similar to the use of A = 0 and A = 1 in Figure 16-19. Once again, we confirm the same trend for pore coating, consistent with the mild linear decrease discussed in Remark 1.
These observations correlate to the properties of probability distribution for permeability. The distribution for CLPS is less tight than for LP. The histograms of K (i) /K 0 for V o = 0.1 are symmetric, but less so than for the larger V o . For the higher values of V o , we see that, occasionally, the pore space is completely clogged, resulting in K (i) = 0 for some realizations. The trend in the mean of (K (i) /K 0 ) i depending on V o follows a similar trajectory to that for the LP model, but with a greater variance in (K (i) /K 0 ) i . These observations are also similar to those for the DNS model BN and the case A = 0 and A = 1, respectively.
Overall, we infer that pore-coating obstructions, such as in biofilm formation in porous media, may cause somewhat predictable changes in permeability for this range of V o , but that a larger V o corresponds to smaller predictability.

Hydrate-Like Obstructions
Next, we tune the CLPS and LP models to produce pore-filling (hydrate-line) obstructions. We use parameters in Tables 3 and 4. Results for CLPS are shown in Figure 22. The CLPS model predicts a drop in relative permeability of about one order of magnitude at V o = 0.1. For the LP model, the histograms of (K (i) /K 0 ) i show a distribution skewed towards the origin, and the change in the mean is sharp, dropping nearly 85% with V o = 0.1, but levels off as V o increases. The standard deviation of (K (i) /K 0 ) i is quite large, suggesting a very significant but less predictable impact on the sediment permeability when hydrate crystals form.
While the results of CLPS and LP are similar, we note that the CLPS model has a smaller spread with smaller standard deviation as compared to the LP model.
Finally, we see that the trend for pore filling is consistent with the dramatic decrease (sharp drop with a convex graph) discussed in Remark 1.

Disaggregated Obstructions
For additional interest, we here present the results for geometries D o without any aggregation, of the "free colloid" type. These can be generated by the LP model with parameters, as in Table 4. We show the results in Figure 23 for the distributions corresponding to V o ∈ {0.05, 0.1, 0.2}. We see that he impact on permeability appears much stronger than in the pore-coating or pore-filling case. We acknowledge that this effect may be somewhat overstated, due to our experiments being conducted in 2D and to the coarse grid resolution. Nevertheless, this is a striking result, which speaks for the need to simulate the aggregation and disaggregation of D o for realistic predictions of permeability.

Flow and Transport Properties in Many-Pore Geometries
In addition to studying the flow in obstructed geometries, it is important to study the transport time at which the outflow reaches some fixed ratio of the inflow of the transported species.
For these studies, we use a large REV with the many-pore geometry from Figure 3. We generate I = 100 obstructed geometries with CLPS and LP models for each pore-coating and pore-filling obstruction with V o ∈ {0.05, 0.1, 0.15}. For each geometry D (i) i , we solve for flow and upscale to the permeability K (i) . We follow up with transport solutions in the Eulerian frame of reference, and approximate the breakthrough curve B (i) (t) as in (9).
We study the scatter plots of K (i) /K 0 vs T (α) (i) /T (α) 0 to show that the latter quantity is not correlated to the former; thus, it is important to study both.
Example 3 (Many-pore with LP model). We use the parameters from Table 4 for many-pore geometries. We present the permeability results in Figure 24 and transport results for α = 0.8 and α = 0.95 in Figure 25, with a least squares fit to the data shown as a solid gray line. In Figure 24, we see good agreement with the single-pore results from Section 6.1. This includes a severe drop in permeability in the case of pore-filling obstructions and a less significant drop in the permeability in the case of pore-coating obstructions. Mean values and standard deviations for these cases are reported in Tables 5 and 6. For pore-coating obstructions, when V o = 0.1, there are 20 geometries that are impermeable and excluded from the mean; this explains the large standard deviation in relative permeability for that case.
For this example, we observe that, for many of the cases, we have T (α) (i) /T (α) 0 < 1. Compared to the intuition on a small or single-pore geometry, this result may appear surprising, since it suggests that transport in the presence of obstructions is happening more quickly. After reflection, in many-pore geometry, the result can be explained by the fact that our experiments use the same inflow value u in in a large REV, and the transport time substantially depends on a larger variety of flow paths than was possible in a singlepore geometry. Table 5. Sensitivity of the probability distribution of K (i) /K 0 to the simulation parameters of the CLPS model on the single-pore and many-pore geometries. The mean and standard deviation of K are presented for three selected vales of the relative volume V o = |D o | |D v | , as shown.  Table 6. Sensitivity of the probability distribution of K (i) /K 0 to the simulation parameters of the LP model on the single-pore and many-pore geometries. The mean and standard deviation of K are presented for several selected vales of the relative volume V o = |D o | |D v | , as shown. For the many-pore geometry and V o = 0.1, the method LP generates 20 out of I = 100 impermeable geometries. For V o = 0.15, most are impermeable. These facts result in a large standard deviation, and these results are not included. Next, for pore-coating obstructions, when α = 0.8, the trend line suggests a correlation between T (α) and K (i) with a slope about of 0.6. However, when α = 0.95, this trend diminishes. We also note an outlier where the permeability has dropped by nearly 60% while the breakthrough time has dropped by less than 20%; this is due to the clogging of certain throats in the REV. We next discuss the results obtained with the CLPS model. Since they are qualitatively similar for permeabilities, we present a different view of transport properties.
Example 4 (Many-pore with CLPS model). We use the parameters with (21) model as in Table 3 for the many-pore example. We present the results for permeability in Figure 26. Additionally, we study the tortuosity T as a function of obstruction volume, which we show in Figure 27. The case T (i) /T 0 > 1 or T (i) /T 0 < 1 corresponds to the increasing or decreasing average path lengths in the presence of obstructions.
We discuss these results here. We see good agreement between the single-pore, in Section 6.1), and many-pore cases. For pore-coating obstructions, there is a gradual decrease in permeability with V o , whereas, in the pore-filling case, the decrease in permeability is sharper. In the pore-filling case, the many-pore geometry with V o = 0.1 is, on average, more than twice as permeable as the single-pore case; see Table 5.
Next, from visual observations of geometries in Figure 27, we see D o for pore-coating featuring two instances of throat clogging and many small "colonies". In addition, we see that T (i) is usually greater than T o in the unobstructed geometry; with a greater obstruction volume V o , there is a greater increase in T . For pore-filling, we see tortuosity T decreasing as V o increases.

Sensitivity of BN and Reduced Models CLPS and LP to Parameters
For any computational model, it is important to study its sensitivity to data. For models whose solutions depend on the data, a formal framework is provided, e.g., through Sobol indices [72] based on multivariate polynomial expansions of the quantity of interest; see, e.g., their use in [22]. This type of analysis requires high regularity of the quantities of interest so that the sensitivities, i.e., "derivatives", can be assessed.
The quantities of interest in our models are (i) the geometries D o and (ii) the upscaled parameters such as K = K(V o ).
For the first, (i) the framework to study (the smoothness of) their dependence on data is largely unavailable, while it loosely connects to the abstract framework of shape analysis in high-dimensional spaces. In the practical setting of this paper, our experiments indicate that the most important qualitative characteristic is whether the geometries produced by our models are of pore-coating or of pore-filling character, but it is unclear whether these observations can be correlated with some metrics that are to be studied quantitatively. Further analyses and theoretical framework are needed; see future work discussed in Section 7.
For (ii), the Darcy scale quantities rely on an averaging process carried out in upscaling, and thus can be believed to be reasonably smooth. However, the dependence of K resulting from DNS models and CLPS and LP models on their parameters raises formal challenges. One challenge is as follows: the permeability values K for some realizations of geometry D o can be zero because the associated geometries may include obstructions critically placed in the flow path. This feature may occur without any apparent dependence on the parameters of the simulations, and is an example of the non-smooth behavior of even the upscaled quantities.
At this time, we are unsure whether the formal framework of sensitivity applies to the quantities (i) and (ii) from our simulations. We comment on the qualitative observations only briefly.

Sensitivity of BN Models
In the DNS model for the pore-coating behavior related to the BN model for biofilm growth, the solutions governed by (12) and the pattern of D o depend on all the data to the PDE model. In particular in the results discussed here, we keep boundary conditions as well as the diffusion and reaction rates fixed using parameters from the literature; see [6,30].
In this paper, we vary the initial conditions through parameter V b and the right-handside function depending on parameter A. As we show in the geometries in Figure 8, and in the plots of permeabilities in Figure 14, the initial condition parameter V b for the DNS simulation plays a significant role in the biofilm growth pattern. For large(r) V b , the growth in biofilm is more evenly distributed near the walls, which results in more pronounced pore-coating behavior. For small V b , the attraction parameter A, when set to A = 1, contributes significantly to the "pore coating" pattern, and eventually results in the dependence of the upscaled quantities consistent with this pattern, as discussed in Remark 1. More detailed data for this sensitivity study are provided in Table 7. Table 7. Sensitivity of the probability distribution of K (i) /K 0 to the simulation parameters V b and A for the geometries D 0 obtained with DNS for pore-coating scenarios. The values K presented here correspond to three selected vales of the relative volume

. Sensitivity of reduced models CLPS and LP
For the reduced models CLPS and LP the control parameters include those discussed in Sections 4 and 5. The CLPS model depends on parameters , δ, θ, as in Table 3, and the choice of is most delicate, as its small perturbations lead to different stationary solutions. The LP models solutions critically depend on w ro , w rv , w vo as shown, e.g., in Table 4. The CLPS model is far more sensitive to its parameters than the LP model.
The choice of parameters for CLPS and LP is guided by heuristics, and the qualitative character of geometries is assessed by visual inspection. Most significantly, the reduced models produce geometries of drastically different characters (pore-coating or pore-filling) due to a small change in the control parameters, i.e., they produce bifurcations. This feature seems similar to the sensitivity of the Ising-like models of phase transitions to their critical parameters, such as the inverse temperature β. We provide data for the sensitivity of the CLPS and LP models in Tables 5 and 6.
Lastly, we provide an aggregated set of figures, which, taken together, illustrate the sensitivity of the simulation results to the initial conditions and to the initial geometry (single-pore or many-pore). See Figure 28, including our commentary on the overall sensitivity and qualitative and quantitative comparison between DNS, CLPS, LP and algebraic models. Figure 28. Illustration of sensitivity of mean reduced permeability to the choice of V o and to the choice of initial pore geometry. The displayed results were obtained by DNS in various variants, and reduced models CLPS and LP for pore coating, and CLPS and LP for pore filling. The algebraic models guide the eye. The data are extracted from Tables 6 and 7, and the marker size is proportional to the standard deviation for the data. The illustrations confirm the observations made in text, which we paraphrase here. The reduced model results generally qualitatively agree with DNS and algebraic models from Remark 1, but simulation models predict a more severe reduction in K than predicted by algebraic models. They have larger uncertainty for pore-coating scenarios and pore-filling scenarios for small V o . Predictions for single-pore and many-pore geometries are similar to each other for small V o , but their mean is not reliable at large V o due to clogging, especially for many-pore geometry.
In the end, we acknowledge the sensitivity to parameters, but do not see it as useful to study these quantitatively within formal frameworks which use differentiation or otherwise rely on the regularity of the quantities of interest. In turn, we believe more theoretical work is needed to study the metrics, smoothness, and overall sensitivity of the solutions on the model parameters.

Computational Complexity of BN and Comparison with Reduced Models CLPS and LP
The overall cost of each method, per realization (i), is the cost of each of pre-processing + generation of D In the estimates below, we assume we have M grid cells, including those in D r . However, an efficient implementation may have lower complexity or be based on φ 0 M rather than M, depending on how we handle the keyed-out elements from D r . Regardless, we have M ≈ 10 4 for single-pore and M ≈ 2 × 10 4 many-pore geometries in d = 2. For large-case simulations in d = 3, which are not reported here, typically, M = O(10 6 ).
We report the cost of the algorithms, which create one realization of D o . The preprocessing cost involves the initialization of algorithms with some initial conditions, which use a random number generator called I times on a vector of size M.
Throughout, we denote the cost of solving a linear system with M unknowns by solver(M). We recall that the cost of direct methods such as LU and QR in the MATLAB environment for dense matrices is solver(M) = O(M η ) with η = 3. On the other hand, the discretizations we consider lead to sparse linear systems, and it is well known that appropriate direct and iterative solvers have a cost of O(pM), where p is band size. On the other hand, the solver in the post-processing environment HybGE-Flow3D has η ≈ 1.7, since it uses state-of-the-art preconditioners, and in d = 3, the band size p ≈ M 2/3 . Once (D (i) o ) i are found, we post-process. The cost of post-processing and upscaling to Darcy scale properties is about solver(2M) + 5M.
In general, we find that the cost of generation of (D (i) o ) i far outweighs that of postprocessing for DNS and CLPS, and is comparable to that for LP. The DNS of BN seems feasible for d = 2, but may be unfeasible for d = 3.
Moreover, we find that the cost of data storage and file IO is already considerable for I = 100, and may require creative approaches when I >> 100.

Complexity of BN
The solver BN for DNS requires the solution of two coupled nonlinear PDEs for (B, N), each of size M. Additionally, we keep track of and update the Lagrange multiplier, also of size M. In other words, we have 3M unknowns: (B, N, Λ). We also pre-compute a(x) prior to simulations.
We solve the problem over 100 to 1000 time-steps, depending on the case, and each requires, on average, about two iterations of the semi-implicit method. We require around 100M operations per iteration, plus the cost solver (3M). The former include computing diffusion coefficients, the Jacobian, residuals, and other operations. For the linear solver, the system is sparse but not symmetric.
We provide approximate times to indicate complexity when solving BN up until plugging. We report on wall clock time on a Linux multicore desktop with MATLAB. For single-pore geometry with V b = 0.07, this took about 30 [s] for A = 0 and 2.2 [m] for A = 1. For many-pore geometry with V b = 0.07, this took 5.3 [m] for A = 0 and 9.2 [m] for A = 1. These times simply indicate that it takes longer to simulate until complete plugging when the obstructions adhere more closely to the walls.

Complexity of CLPS
As compared to BN, the CLPS algorithm only solves one variable ψ, and uses linear symmetric non-degenerate diffusion and time-lagging, without any additional iterations. The localization functions are pre-computed. Thanks to these reductions, the complexity is considerable compared to BN.
The cost per time-step is about 5M operations, which icludes dealing with the volume constraint plus the cost of the linear solver solver(M + 1). We run time-steps until a stationary solution is found, which ranges from 100 to 10, 000.
Depending on the case, single-pore simulations run for about 10 [s] on average.

Complexity of LP
The LP algorithm might take as many as O(10 5 ) time-steps to reach the stationary solution, but each is very inexpensive, and requires only the cost of swaps and recalculation of Hamiltonian update ∆H(X). We approximate this cost to be 10 operations per time-step.
Depending on the case, LP simulations for single-pore on a laptop take less than 1[s] on average, and are about 6 times faster than CLPS.

Summary and Future Work
Imaging and other high-resolution data pose an opportunity and, simultaneously, a challenge to computational science. In this paper, our focus is on pore-scale data for flow and transport, and their upscaling to Darcy-scale properties of porous media, such as permeability and tortuosity. Specifically, we consider a medium which changes due to clogging, with two patterns of obstructions in pore-scale geometries, referred to as pore-coating and pore-filling. Both are motivated by their important applications to biofilm and hydrate crystal growth, respectively.
Interpreting the static imaging data as grids for flow simulations and upscaling is possible. However, imaging of pore-scale which changes in time due to clogging poses additional difficulties, and handling the potentially large collection of geometries is extremely computationally intensive, and requires calibration and tuning. Using DNS for time-dependent processes instead of imaging is quite complex and comes with considerable uncertainty.
In this paper, we presented an alternative reduced model, which can create off-line probability distributions of properties of interest without the DNS. In particular, we present two algorithms which generate plausible geometries with obstructions called D o , which clog the pores. These algorithms are not simulations of the actual physical process, but rather simulate the presence of D o which appear similar to those observed in experiments. The methods have roots in mathematical and statistical models of phase separation; they are easily tunable and efficient.
In the paper, we demonstrated that our reduced models produce comparable results to those of DNS, both by visual inspection of D o and when considering the upscaled quantities. Overall, the CLPS model is about one order of magnitude more efficient than DNS, and the LP model is about one further order of magnitude more efficient than CLPS. However, as we move from DNS to CLPS to LP models, we also move farther and farther away from the actual physical model of generating obstructions.
More work is needed to understand this trend and how to use and calibrate efficient reduced models which maintain close agreement with physical simulations. At the same time, validation efforts are needed to further calibrate the DNS models, which can serve as the basis for reduced models. In addition, theoretical analyses, and, specifically, new metrics, are needed to compare the images we generate with reduced models to the DNS and to the imaging data. In this paper, we used judgement made by the human eye, but we are currently exploring various emerging data science approaches which can be used to analyze the similarity between shapes.
Finally, we aim to continue the validation efforts of the probability distributions we generate with the experimentally obtained Darcy-scale quantities. In particular, in this paper, we considered only two representative geometries, called single-pore and manypore. Further testing and validation efforts should include testing the reduced models in a collection of pore geometries ranging from simple to complicated. Since the permeability depends very much on the local arrangement of the pores and throats, it depends not just on the size but also on the choice of an REV, and this randomness is even more pronounced when the geometry is modified by reactive transport starting from random initial conditions. These facts motivate this paper, since reduced models would produce more reliable results than, say, algebraic formulas, such as those in Remark 1; we observed this fact in [6].
We believe that it would be useful to create a library of probability distributions corresponding to the obstructed domains classified by type, say, channel-like, or large porosity (glass-beads), or sand-stone geometries. These efforts may be connected to some theoretical analyses of shapes, which we hope to address in the future.
Our current work also includes further calibration of the models to discover the optimal values of parameters for CLPS and LP starting from the available data. In this direction, we plan to explore data science approaches. There are also many interesting application-specific challenges, including the dependence of the results on I, the size and the shape of the rock domain D r in the selected REV, and the voxel resolution h. Example A1. We set = 0.005, θ = 10, and δ = 0.1 and approximate the solutions to (21). We let the solution evolve until the time T ∞ , after which the solutions remain essentially stationary. In Figure A1, we show the two stationary states corresponding to i = 1 and i = 2. In Figure A2, we illustrate the well-known impact of parameter , when the initial state ψ init (x) = ψ (2) init (x), θ = 1 and δ = 0.1.  Next, we study the impact of the parameter θ and of the localization functions.