# Interplay of biofilm growth, NAPL biodegradation and micro-scale heterogeneity in natural attenuation of aquifers delineated by pore-network modelling

Morteza Aminnaji, Fabrice Golfier, Vahid Niasar, Masoud Babaei

► To cite this version:

Morteza Aminnaji, Fabrice Golfier, Vahid Niasar, Masoud Babaei. Interplay of biofilm growth, NAPL biodegradation and micro-scale heterogeneity in natural attenuation of aquifers delineated by pore-network modelling. Advances in Water Resources, 2020, 145, pp.103750. 10.1016/j.advwatres.2020.103750. hal-02966298

HAL Id: hal-02966298

https://hal.univ-lorraine.fr/hal-02966298v1

Submitted on 26 Sep 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Interplay of biofilm growth, NAPL biodegradation and micro-scale heterogeneity in natural attenuation of aquifers delineated by pore-network modelling

Morteza Aminnaji<fnref n="1" />, Fabrice Golfier<fnref n="2" />, Vahid J. Niasar<fnref n="1" />, Masoud Babaei*<fnref n="1" />

## Key Points

* Bioremediation of NAPL pollutant was studied at pore scale using pore network modelling.

* NAPL dissolution and biodegradation in the natural attenuation process were studied in presence of pore-scale spatial correlation between pore radii.

* Spatially pore-size correlated networks showed lower NAPL dissolution and biodegradation.

* Uncorrelated pore networks with well-distributed NAPL through the entire pore network led to higher bio-growth resulting in more reduction in permeability and porosity.

## Abstract

The biofilm plays an important role in NAPL bioremediation, natural attenuation and the flow behaviour of contaminant in porous media. Therefore, NAPL biodegradation in the presence of any bacteria in the porous media needs to be considered for more accurate prediction of contaminant remediation from the soil and groundwater. The complex interaction between the NAPL biodegradation/dissolution and physiochemical bacteria (biofilm) growth — which depends on both time and space — results in difficulties in simulation of NAPL bioremediation in the porous media. In this work, the natural attenuation process is simulated using three-dimensional pore network models to investigate natural bioremediation and dissolution of NAPL in presence of different heterogeneities arising from uncorrelated and correlated pore radii distribution. The impact of an additional phase (biofilm) on transport and NAPL dissolution is dynamically incorporated in the pore network model with simplifying assumptions such as fully available NAPL for biodegradation and dissolution. Under such assumptions, the results indicate that as pore scale heterogeneity (in the form of correlation length of pore radii of pore networks) increases, the biofilm growth decreases resulting in a reduction of NAPL bioremediation and dissolution. As such it is critical to consider pore scale heterogeneity in predicting bioremediation efficiency. Additionally, we found that Kozeny-Carman equation fail to predict the permeability changes due to biofilm growth/extinction.

<footnote n="1">¹Department of Chemical Engineering and Analytical Science, the University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote n="2">²Université de Lorraine, CNRS, CREGU, GeoRessources Laboratory, 54518 Vandoeuvre-lès-Nancy Cedex, France</footnote>

<footnote>* Corresponding author: telephone: +44 (0)161 306 4554, email: masoud.babaei@manchester.ac.uk</footnote>
Nomenclature

<table><tr><td>A<sub>ij</sub></td><td>throat cross-section area between pore i<sup>th</sup> and j<sup>th</sup> (L<sup>2</sup>)</td><td>l</td><td>throat length (L)</td></tr><tr><td>C</td><td>concentration (ML<sup>-3</sup>)</td><td>M<sup>t</sup><sub>NAPL,i</sub></td><td>mass of NAPL in the pore i<sup>th</sup> at the times t (M)</td></tr><tr><td>C<sub>A</sub></td><td>concentration of species A (ML<sup>-3</sup>)</td><td>N<sub>i</sub></td><td>number of pores connected to pore i<sup>th</sup></td></tr><tr><td>C<sub>B</sub></td><td>concentration of species B (ML<sup>-3</sup>)</td><td>N<sub>n,i</sub></td><td>number of NAPL-filled pores connected to pore i<sup>th</sup></td></tr><tr><td>C<sup>t</sup><sub>ave</sub></td><td>average concentration of species at time t (ML<sup>-3</sup>)</td><td>N<sup>t</sup><sub>colony,i</sub></td><td>number of micro-colonies in pore i<sup>th</sup></td></tr><tr><td>C<sup>t</sup><sub>i</sub></td><td>NAPL concentration in pore i<sup>th</sup> at time t (ML<sup>-3</sup>)</td><td>P</td><td>pressure (ML<sup>-1</sup>T<sup>-2</sup>)</td></tr><tr><td>C<sup>*t</sup><sub>ij</sub></td><td>concentration of pore j<sup>th</sup> connected to pore i<sup>th</sup> if the flow direction is from pore j<sup>th</sup> to i<sup>th</sup>; and the concentration of pore i<sup>th</sup> if the direction is from pore i<sup>th</sup> to j<sup>th</sup></td><td>PVI</td><td>pore volume injected</td></tr><tr><td>ΔC<sup>t</sup><sub>ij</sub></td><td>concentration difference between pore i<sup>th</sup> and j<sup>th</sup></td><td>Q<sub>ij</sub></td><td>aqueous flow rate between pore i<sup>th</sup> and j<sup>th</sup> (L<sup>3</sup>T<sup>-1</sup>)</td></tr><tr><td>C<sub>s</sub></td><td>equilibrium concentration, (ML<sup>-3</sup>)</td><td>R<sub>A</sub></td><td>specific growth rate, nonlinear Haldane kinetic reaction term for species A (T<sup>-1</sup>)</td></tr><tr><td>D</td><td>diffusion coefficient (L<sup>2</sup>T<sup>-1</sup>)</td><td>r<sub>i</sub></td><td>radius of pore i<sup>th</sup> (L)</td></tr><tr><td>K<sub>abs</sub></td><td>absolute permeability (L<sup>2</sup>)</td><td>r<sub>ij</sub></td><td>radius of throat between pore i<sup>th</sup> and j<sup>th</sup> (L)</td></tr><tr><td>K<sub>effec, with bio</sub></td><td>effective permeability (L<sup>2</sup>) when biofilm is assumed as an impermeable/solid phase</td><td>t</td><td>time (T)</td></tr><tr><td>K<sub>effec, no bio</sub></td><td>effective permeability (L<sup>2</sup>) when biofilm is assumed as a very high permeable phase</td><td>V<sub>colony</sub></td><td>volume of each bacteria colony (L<sup>3</sup>)</td></tr><tr><td>K<sub>ext</sub></td><td>extinction coefficient of bacteria (T<sup>-1</sup>)</td><td>V<sub>i</sub></td><td>pore volume (L<sup>3</sup>)</td></tr><tr><td>K<sub>i</sub></td><td>inhibition coefficient (M L<sup>-3</sup>)</td><td>V<sup>t</sup><sub>bio,i</sub></td><td>Volume of bio in pore i<sup>th</sup> at time t (L<sup>3</sup>)</td></tr><tr><td>K<sub>m,A</sub></td><td>half-saturation constant for species A biodegradation rate (M L<sup>-3</sup>)</td><td>Y<sub>bio</sub></td><td>yield or stoichiometric coefficient (mg of biomass per mg of NAPL)</td></tr><tr><td>K<sub>m,B</sub></td><td>half-saturation constant for species B biodegradation rate (M L<sup>-3</sup>)</td><td>Δt</td><td>increment in time (T)</td></tr><tr><td>K<sub>part</sub></td><td>partitioning coefficient between the NAPL and the biofilm phase</td><td>μ</td><td>viscosity (M L<sup>-1</sup>T<sup>-1</sup>)</td></tr><tr><td></td><td></td><td>μ<sub>bio</sub></td><td>maximum reaction rate (specific degradation rate) for the species A (T<sup>-1</sup>)</td></tr><tr><td></td><td></td><td>ρ<sub>bio</sub></td><td>biomass density (wet mass per volume) (M L<sup>-3</sup>)</td></tr><tr><td></td><td></td><td>φ</td><td>porosity (fraction)</td></tr></table>
# 1 – Introduction

Contamination of soil and aquifer by non-aqueous-phase liquids (NAPLs) and other industrial waste are encountered because of their ubiquitous use, unintended release and, perhaps poor disposal. Several remediation/removal options have been used to remove NAPL from soil and groundwater including ex-situ and in-situ treatment (biological, chemical, and physical). In the ex-situ treatment, the source of contamination is taken out from the soil and then treated. Ex-situ treatment is not always an efficient method, so in-situ NAPL remediation methods have instead been proposed including chemical oxidation (Soga et al., 2004, Huling & Pivetz, 2006), surfactant-enhanced aquifer remediation (SEAR) (Fountain et al., 1996, Mulligan et al., 2001), bioremediation (Singh et al., 2006), vapour extraction (Heron et al., 2005), and natural attenuation (Wiedemeier et al., 1999, Lari et al., 2019). In the natural attenuation process, NAPL is naturally removed and remediated from the groundwater and no specific engineering technologies are employed except NAPL monitoring. Natural attenuation process includes various mechanisms such as biodegradation, volatilization, sorption, and natural dissolution.

Among the abovementioned technologies and options, bioremediation is one of the cost-effective and environmentally friendly technologies for pollutant remediation (Bhatnagar & Kumari, 2013). In the biodegradation process, the microorganism breakdown of the contaminants to smaller and non-toxic compounds occur (Van Dillewijn et al., 2007, Bhatnagar & Kumari, 2015). In the natural biodegradation process in aquifer, microorganisms are naturally present as a biofilm at the surface of the grains (Orgogozo et al., 2013). In some applications bio-surfactant or nutrients are injected to the porous media to produce microorganism for NAPL bioremediation or oil mobilisation in reservoir similar to processes involved in microbial-enhanced oil recovery or MEOR (Banat, 1995, Sivasankar & Kumar, 2014, Soudmand-asli et al., 2007, Lazar et al., 2007, Brown, 2010).

The biofilm growth in the porous media (mostly without NAPL phase) has been explored visually using x-ray micro-tomography or X-ray CMT (Davit et al., 2011, Larue et al., 2018, Ostvar et al., 2018). Vayenas et al., 2002, studied NAPL biodegradation through injection of dissolved NAPL into biofilm populated porous media constructed of etched glass. They used quasi-steady-state theoretical model of biodegradation to calculate mesoscopic biochemical rate and thickness of biofilm. However no separate
NAPL phase was considered. Significant research has been dedicated to model in-situ NAPL bioremediation. Transport and biodegradation of dissolved NAPL plume are coupled together with microbial growth and described by advection-dispersion-reaction equations at the continuum scale where Darcy's law is valid (Yang et al., 1995, Brauner & Widdowson, 2001, Clement et al., 2004, Becker & Seagren, 2009, Davit et al., 2010, Chambon et al., 2010, Orgogozo et al., 2010, Manoli et al., 2012, Orgogozo et al., 2013). However, research is lacking and limited to the mechanisms driving NAPL bioremediation at the pore-scale. That is, although pore scale and pore network models have been developed for implementation of biofilm growth by many researchers (Suchomel et al., 1998, Thullner et al., 2002, Thullner & Baveye, 2008, Ezeuko et al., 2011, Pintelon et al., 2012, Tang et al., 2013, Peszynska et al., 2016, Benioug et al., 2017), they have not considered the NAPL bioremediation itself.

Amongst the existing research studies, Qin and Hassanizadeh (2015), used Hagen-Poiseuille equation to calculate the steady-state water flux and studied the solute transport and the biofilm growth in a pore network model. They investigated the permeability reduction due to biofilm growth, biofilm clogging, and non-equilibrium mass exchange of solute between the water and biofilm phases. Von der Schulenburg et al. (2009) used three-dimensional pore-scale model of biofilm growth in porous media based on a lattice Boltzmann simulation for prediction of permeability as a function of biofilm accumulation. Moreover, a pore-scale description of NAPL transport and biodegradation in a porous medium containing biofilm was used and up-scaled for the relevant Darcy-scale equations using the method of volume averaging (Bahar et al., 2016, Golfier et al., 2009). Schmidt et al., 2018, conducted two-dimensional (2D) finite element simulations of fluid flow and solute transport with reaction to show the effect of bacterial distribution on the biodegradation. Recently, Benioug et al. (2019) derived an immersed boundary-lattice Boltzmann model to study the interaction between biofilm growth and NAPL dissolution at the pore-scale. However, this study was limited to 2D geometry and a small number of pores, most likely due to computational expenses in simulation.

In this work, particular focus will be given to fill the lack of knowledge related to the impacts of “heterogeneity in the form of spatial correlated pore radii distribution” on NAPL biodegradation,
dissolution, and biofilm growth. To the best of our knowledge, use of pore network for natural
bioremediation and dissolution has not been reported previously for micro-heterogeneous porous media.
Therefore, this study simulates the bioremediation using the unstructured three-dimensional pore
network model in which the effect of pore-scale heterogeneity (in the form of spatially correlated and
uncorrelated pore radii distribution) on NAPL dissolution/biodegradation, biofilm growth, and long
term natural remediation are investigated. In addition, the impact of the additional phase (biofilm
growth) on transport and NAPL dissolution is investigated in the pore network model.

## 2 – Methodology

### 2.1. NAPL biodegradation and dissolution modelling

Bacteria and microbes usually aggregate and adhere to the solid surfaces and form a cluster of microbes as a biofilm attached to the surface (Costerton et al., 1999). The impact of cell morphology on biofilm architecture and competition was studied by Smith et al., 2017. The biofilms exist in the pores and throats where a sufficient nutrient is available. As shown in **Figure 1**(a), the immobile biofilm phase is considered in the pore network that adheres to the pore and throat surfaces (the biofilm is considered as a uniform thickness in the pores and throats). Bacterial populations are considered to grow as biofilm (biomass phase) on the walls of throats and water-filled pores (on the rock surface). The biomass phase in the NAPL-filled pore is assumed to be made of micro-colonies (Baveye & Valocchi, 1989) which are sparsely distributed on the NAPL blob surfaces within NAPL-filled pores (in this work, NAPL is considered as trapped only in the pores). While planktonic bacteria are neglected in this work, surfactant-producing bacteria and inhibition effect on bacterial growth due to NAPL toxicity are considered.

There have been numerous research studies to find the mechanism involved in NAPL dissolution at the pore scale, and different models have been used for the pore network modelling. For example, Held & Celia, (2001), assumed the NAPL dissolution through the stagnant-layer (e.g. through the stagnant throat by considering a linear driving force for the mass transfer between NAPL-filled pores and water-filled pores). However, they did not consider the NAPL dissolution through the corner flow of wetting
fluid. This method has been used by other researchers including (Agaoglu et al., 2016, Sarikurt et al., 2017, Aminnaji et al., 2019). We would like to note that the assumption of stagnant-layer concept for NAPL dissolution could be over-simplifying as Sahloul et al., (2002) showed experimentally that the corner flow of water phase (wetting phase) controls the rate of NAPL dissolution. Additionally, there are more complicated and more inclusive modelling methods in literature that had accounted for physical processes such as corner flow, corner diffusion and convection diffusion, (Dillard & Blunt, 2000, Zhou et al., 2000, Zhao & Ioannidis, 2003, Zhao & Ioannidis, 2007b, Zhao & Ioannidis, 2007a, Zhao & Ioannidis, 2011, Huang et al., 2015, Khasi et al., 2020). Moreover, ganglion movement and snap-off during the NAPL dissolution and biodegradation are important mechanisms which need to be considered to have a more accurate pore network modelling for the prediction of contaminant remediation.

The concept of ganglion movement and snap-off requires the implementation of a dynamic pore network modelling (Chrysikopoulos & Vogler, 2006, Joekar-Niasar et al., 2010, Joekar-Niasar & Hassanizadeh, 2011, Joekar-Niasar & Hassanizadeh, 2012, Aghaei & Piri, 2015). In this work, however, we use the concept of NAPL dissolution through the stagnant throat in which the corner flow of aqueous phase in the NAPL-filled flow is ignored. We use a quasi-static pore network model with an immobile NAPL phase. Acknowledging the simplicity or over-simplicity of this modelling approach (with further implications on results that will be discussed in Section 3), we point out that although the results of each pore network modelling algorithms (either using concept of stagnant NAPL dissolution or concept of advection-corner NAPL dissolution) will be different, the incremental impact of biodegradation in the NAPL remediation in pore network structures *may* become most likely the same. Moreover, the particular focus of this work is to highlight implementation of NAPL biodegradation in the pore network modelling under the heterogeneity at pore-scale, and their interplay with dissolution, and biofilm growth, albeit with a simple pore network modelling. That is, in a comparison basis, the difference of biodegradation impact between two heterogeneous pore networks *may* not significantly change from one algorithm to another. However, further work in future is needed to fully delineate the impact of microscopic heterogeneity on biodegradation using dynamic pore network modelling with corner
diffusion incorporated. To this end, in this work it is assumed that the NAPL dissolution into water
occurs where the NAPL-water interface exists in the stagnant throats only and corner diffusion model
is ignored. In such throats the aqueous flow rate is zero and the advection regime does not take place.

A linear driving force is considered for mass transfer from NAPL-filled pore to the water-filled pore
and the NAPL mass transfer is described using the term $DA_{ij}(C_s - C_i^t)/l$ in the stagnant throats, where
$D$ is the diffusion coefficient, $A_{ij}$ is the throat cross-section area between pore $i^{th}$ and $j^{th}$, $C_s$ is the
solubility of NAPL in the water, and $l$ is the throat length. This linear mass transfer simplifies the
multistep partitioning and diffusion of NAPL species from the NAPL phase to the aqueous phase (three
steps of interphase mass transfer as shown in **Figure 1**(b) corresponding to (1) diffusion of NAPL
species in its own phase from $C_{\text{NAPL}}^n$ to $C_{\text{NAPL}}^{\text{int}}$ (NAPL concentration at the interface), (2) partitioning
through thermodynamic equilibrium in the interface from $C_{\text{NAPL}}^{\text{int}}$ to $C_s$ (solubility/equilibrium
concentration), and (3) diffusion of NAPL species in the aqueous phase from $C_s$ to the actual
concentration of NAPL in the aqueous phase, $C_i^t$. The approach has been previously employed
elsewhere (Jia et al., 1999, Held & Celia, 2001, Detwiler et al., 2001). It is assumed that the first two
steps are not rate-limited. When a biofilm exists on the NAPL surfaces, while NAPL can be degraded
within the biofilm, an extra mass transfer limitation will occur (Mulder et al., 1998). For the sake of
simplicity, this mass transfer resistance is discarded and a micro-colony description is considered on
the NAPL surface to prevent the total covering of NAPL blob by biofilm. Therefore, for the diffusion
of NAPL into the biofilm for the biodegradation, it is assumed no rate limited in diffusion of NAPL
species in its own phase. The partitioning of nutrient from the NAPL phase to the biofilm phase is
schematically shown in **Figure 1**(c). Although there is mass transfer limitation occurring within a
boundary layer between biofilm and the bulk liquid (Pizarro et al., 2001), we assume an instant
partitioning of NAPL into biofilm using a partitioning coefficient. However, a limitation for biomass
build up is considered based on the amount of NAPL exists and the space volume available in the pores
and throats

The diffusion and partitioning of NAPL species into the aqueous and biofilm phases, NAPL dissolution
and biodegradation (in the throats, water-filled pores, and NAPL-filled pores) are numerically simulated
within three-dimensional unstructured pore networks. In addition, biofilm could grow in each throat, water-filled pore, and NAPL-filled pore dynamically affecting the aqueous phase velocity. It is assumed that water exists in the stagnant throats which connect NAPL-filled pore and water filled pore, i.e., a possible layer of water could exit at the interface of NAPL-filled pore in which provide a possible growth of biofilm. The following sections describe the numerical modelling of the NAPL dissolution and biodegradation, and biofilm phase growth at the pore scale.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) Schematic of an aquifer contaminated by NAPL in the presence of biofilm together with (b) a schematic of NAPL species diffusion process to the aqueous phase and (c) NAPL partitioning into the biofilm phase.</figcaption></figure>

### 2.1.1. NAPL transport in the aqueous phase

The methodology presented in the following has been used elsewhere (Held & Celia, 2001, Agaoglu et al., 2016, Sarikurt et al., 2017, Aminnaji et al., 2019). The NAPL phase is assumed as an immobile phase in which the aqueous phase does not flow in the NAPL-filled pores, i.e., the water flow is restricted to the water-filled pores. The pressure distribution through the pore network is calculated using the volumetric flow balance equation for each of the pores:

$$ \sum_{j}^{N_i} Q_{ij} = 0 \qquad \text{Eq. 1} $$

where $Q_{ij}$ is the flow rate between pore $i^{th}$ and $j^{th}$, and $N_i$ is the number of pores which are connected to pore $i^{th}$. The throat flow rate which connects pore $i^{th}$ to $j^{th}$, $Q_{ij}$, is calculated using Hagen-Poiseuille equation, i.e., no pressure drop is considered in the pore bodies:
$$Q_{ij} = \frac{\pi r_{ij}^4}{8\mu l_{ij}} \Delta P_{ij} \qquad \text{Eq. 2}$$

where $r_{ij}$ is the throat radius, $l_{ij}$ is the throat length, $\mu$ is the viscosity of the aqueous phase, and $\Delta P_{ij}$ is a pressure difference between two sides of the throat. Combining Eq. 2 and Eq. 1 for each of the pores gives a linear system of equations for pressure:

$$\sum_{j}^{N} \frac{\Delta P_{ij} r_{ij}^{4}}{l_{ij}} \Delta P_{ij} = 0 \qquad \text{Eq. 3}$$

After simulation of pressure distribution and calculation of flow rates in the throats, the dissolved NAPL transport is simulated using the advection-diffusion equation:

$$V_i \frac{C_i^{t+1} - C_i^t}{\Delta t} = \sum_j^{N_i} Q_{ij} C_{ij}^{*t} + \sum_j^{N_i} \frac{DA_{ij}}{l_{ij}} \Delta C_{ij}^t + \sum_j^{N_{n,i}} \frac{DA_{ij}}{l_{ij}} (C_s - C_i^t) \qquad \text{Eq. 4}$$

where, $C_{ij}^{*t}$ is the concentration of pore $j^{th}$ which is connected to pore $i^{th}$ if the aqueous phase does flow from pore $j^{th}$ to $i^{th}$ and vice versa. The other parameters are defined in Nomenclature. The last term in the right-hand side of Eq. 4 describes the NAPL mass transfer from NAPL-filled pore to the water-filled pore.

### 2.1.2. Addition of biofilm on water-filled pores and throats

To implement NAPL biodegradation and biofilm growth in the pore network, the following assumptions are considered:

* Biofilm may occupy both pores and throats, that is in the throats and water-filled pores, biofilm is coating pore walls and is found as an immobile layer of uniform thickness which reduces the tube conductance accordingly. In the NAPL-filled pores, a micro-colony description (Baveye & Valocchi, 1989) is considered for NAPL biodegradation and biofilm-growth. This is to avoid the biofilm to totally cover the NAPL blob and to maintain the water-NAPL interface.

* Mono-species biofilm is considered where the kinetic of NAPL biodegradation is described using Haldane model which considers the effect of inhibitory substrate on biodegradation kinetics (Okaygun et al., 1992, Saravanan et al., 2008). Haldane model is defined as:

$$R_A = \mu_{\text{bio}} C_B^* \frac{C_A}{C_A + K_{m,A} + \frac{C_A^2}{K_i}} \qquad \text{Eq. 5}$$

$$C_B^* = \frac{C_B}{C_B + K_{m,B}} \qquad \text{Eq. 6}$$
where, $R_A$ [s⁻¹] is (specific) growth rate, $\mu_{bio}$ [s⁻¹] is the maximum specific degradation rate for the species A, $K_i$ [kg.m⁻³] is the inhibition coefficient, $K_{m,A}$ [kg.m⁻³] and $K_{m,B}$ [kg.m⁻³] are the half-saturation constants for the species A and B, respectively. This model considers a system with a single substrate (A: carbon and energy source) and a single electron acceptor (B). A large excess of electron acceptor (e.g. Oxygen) is assumed, so $C_B$ is constant and $C_B^* \approx 1$ (Golfier et al., 2009, Orgogozo et al., 2013, Bahar et al., 2016). Such conditions may occur, for instance, close to the water table for LNAPL aerobic biodegradation if re-aeration of the aquifer from the vadose zone is sufficient but it should be kept in mind that this assumption may lead to overestimating the impact of bacterial activity.

* It is assumed that the biofilm density ($\rho_{bio}$) is constant, also internal and external mass transfer limitations are discarded for biofilm (biofilm is fully penetrated).

* Bio-surfactant production is only considered when biofilm is in contact with NAPL blobs (i.e., in NAPL-filled pores). In addition, concentration of bio-surfactants is assumed to be uniform at the surface of NAPL blobs and transport of bio-surfactant from NAPL-filled pores to throats and water-filled pores is neglected.

Based on the above assumptions, the NAPL transport equation, Eq. 4, is modified for implementation of NAPL biodegradation as:

$$V_i \frac{c_i^{t+1} - c_i^t}{\Delta t} = \sum_j^{N_i} Q_{ij} c_{ij}^{*t} + \sum_j^{N_i} \frac{DA_{ij}}{l_{ij}} \Delta c_{ij}^t + \sum_j^{N_{n,i}} \left( \frac{DA_{ij}}{l_{ij}} (C_s - c_i^t) \right) - \sum_j^{N_i} \left( \frac{R_j(c_{ij}^t)}{Y_{bio}} \rho_{bio} V_{bio,j}^t \right) - \sum_j^{N_{n,i}} \left( \frac{R_j(C_s)}{Y_{bio}} \rho_{bio} V_{bio,i}^t \right) - \frac{R_i(c_i^t)}{Y_{bio}} \rho_{bio} V_{bio,i}^t \quad \text{Eq. 7}$$

$$R_j(C_{ij}^t) = \mu \frac{c_{ij}^t}{c_{ij}^t + K_m + \frac{c_{ij}^t}{K_i}} \quad \text{Eq. 8}$$

$$R_j(C_s) = \mu \frac{c_s}{c_s + K_m + \frac{c_s^2}{K_i}} \quad \text{Eq. 9}$$

The two first additional terms in the right-hand side of Eq. 7 represent biodegradation of dissolved NAPL in the throats which connect pore iᵗʰ to water-filled pores and NAPL-filled pores. The first additional term is only applied when pore iᵗʰ is in the downstream (flow direction is to the pore iᵗʰ). The last additional term describes biodegradation within the water-filled pores. In Eq. 7, $Y_{bio}$ is the yield or
stoichiometric coefficient (mg of biomass per mg of NAPL), and $V_{bio}^t$ is the volume of biofilm in the throat $j^{th}$ or in the water-filled pore $i^{th}$ which is described as:

$$V_{bio}^{t+1} = V_{bio}^{t} + (R - K_{ext})V_{bio}^{t}\Delta t \qquad \text{Eq. 10}$$

where, $K_{ext}$ [s⁻¹] is the extinction coefficient of bacteria and R is described by Eq. 8 or Eq. 9. However, the size of water-filled pores and throats are reduced due to biofilm growth, so their radii are updated as:

$$r_{\text{water-filled pore},i}^{t+1} = \left( (r_{\text{water-filled pore},i}^{t=0})^3 - \frac{V_{\text{bio},i}^{t+1}}{4\pi/3} \right)^{1/3} \qquad \text{Eq. 11}$$

$$r_{\text{throat},ij}^{t+1} = \left( (r_{\text{throat}}^{t=0})^2 - \frac{V_{\text{bio},ij}^{t+1}}{\pi l_{ij}} \right)^{1/2} \qquad \text{Eq. 12}$$

### 2.1.3. Effect of biofilm on NAPL-filled pores and NAPL reduction

Biodegradation does not occur within the NAPL phase, i.e., water and oxygen (or any other electron acceptor) which are necessary for biodegradation are not present in the NAPL phase. Although bacteria may only survive within the aqueous phase, it may grow at the NAPL interface where substrate concentrations are highest. This mechanism is described in Atlas, 1981, and used in another NAPL biodegradation model (De Blanc et al., 1996). The only limitation of biomass growth in our modelling is the amount of NAPL in NAPL-filled pores. We discuss about the implications of this limitation in Section 4.

Volume and mass of NAPL in pores are reduced due to NAPL dissolution and biodegradation until the NAPL in the pores are completely removed and the pore is filled with water (i.e., when a NAPL-filled pore is completely remediated). While the NAPL biodegradation is applied to all NAPL-filled pores, NAPL dissolution is only applied to those NAPL-filled pores which are connected to the water-filled pores as the corner flow was not considered in this modelling. At this stage, the aqueous fluid flow over the pore network should be updated as described in Section 2.1.1. Because of NAPL toxicity, bacteria are assumed to be sparsely distributed at the NAPL blob surface as micro-colonies. The NAPL blob reduction due to dissolution and biodegradation in the NAPL-filled pore $i^{th}$ is given by:
$$M_{NAPL,i}^{t+1} = M_{NAPL,i}^{t} - \sum_{j}^{K} D A_{ij} \frac{c_{s}-c_{ij}^{t}}{l} \Delta t + \frac{R_{i}(K_{\text{part}}C_{\text{NAPL}})}{Y_{\text{bio}}} \rho_{\text{bio}} N_{\text{colony},i}^{t} V_{\text{colony}} \Delta t \quad \text{Eq. 13}$$

where, $M_{NAPL,i}^t$ and $M_{NAPL,i}^{t+1}$ represent the mass of NAPL in the pore $i^{th}$ at the times $t$ and $t + 1$, $N_{colony,i}^t$ is the number of micro-colonies occupying the pore $i^{th}$, and $V_{colony}$ is the volume of each bacteria colony which is assumed constant. In Eq. 13, $K_{part}C_{NAPL}$ is the concentration in the thermodynamic equilibrium in the biofilm phase in which $C_{NAPL}$ is the mass concentration of pure NAPL and $K_{part}$ is the partitioning coefficient between the NAPL and the biofilm phase. Although $K_{part}C_{NAPL}$ is not necessarily equal to $C_s$ when bio-surfactant is produced, for sake of simplicity, $R_i(K_{part}C_{NAPL})$ is calculated using Eq. 9. The growth of micro-colonies is predicted using Eq. 14 and as the number of micro-colonies increases, the radius of NAPL-filled pores decreases and must be updated using Eq. 15. It is assumed that the reduction of NAPL-volume is roughly equal to the increase of biomass volume in the NAPL-filled pores.

$$N_{\text{colony},i}^{t+1} = N_{\text{colony},i}^{t} + (R_{i}(K_{\text{part}}C_{\text{NAPL}}) - K_{\text{ext}})N_{\text{colony},i}^{t} \Delta t \quad \text{Eq. 14}$$

$$r_{\text{NAPL-filled pore},i}^{t+1} = \left( (r_{\text{NAPL-filled pore},i}^{t=0})^3 - \frac{N_{\text{colony},i}^{t+1} V_{\text{colony}}}{4\pi/3} \right)^{1/3} \quad \text{Eq. 15}$$

## 2.2. Simulation parameters

### 2.2.1. Properties of pore networks

The pore networks are generated according to the algorithm used in (Babaei & Joekar-Niasar, 2016) with two realizations of uncorrelated and two realization of correlated networks (statistically identical realizations) as described briefly below. Also **Table 1** summarises the correlation lengths and geometric properties of four networks generated in this work. Note that all these networks have similar permeability but their pore size distribution differ and so does the initial NAPL distribution. These pore networks have $\mu_R = 16.65 \ \mu m$ (mean of pore radii), $\sigma_R = 9.57 \ \mu m$ (standard deviation of pore radii), and $R \le 50 \ \mu m$ (pore radii) with the average coordination number of 5.3, making the networks’
correlation lengths close to the Berea sandstones†. There are examples of contaminated sandstone aquifers across the world (Wealthall et al., 2001, Gooddy et al., 2002, Powell et al., 2003, Javanbakht & Goual, 2016). We can obtain the mean and variance of the natural logarithm of R ($\mu_Y$ and $\sigma_Y$) as functions of the mean and variance of R ($\mu_R$ and $\sigma_R$) by $\mu_Y = \ln(\mu_R^2/\sqrt{\mu_R^2 + \sigma_R^2})$ and $\sigma_Y^2 = \sqrt{\ln(\sigma_R^2/\mu_R^2 + 1)}$. The maximum pore radius in these networks is 50 µm with a minimum distance of 100 µm between each pair of pores (body centres) (i.e., the space is filled randomly with pores). The average distance between the centres of pore bodies of these networks is 119 µm (i.e., this distance is called as the characteristic length of the pore network (Bijeljic et al., 2004)). This average characteristic length is close to the Berea sandstone with a characteristic length of 131 µm and other modelling studies with the range of 100 to 150 µm (Øren & Bakke, 2003, Bijeljic et al., 2004, Mostaghimi et al., 2012). Additionally the correlation lengths reported in Babaei & Joekar-Niasar (2016), correspond to the typical pore scale correlation lengths observed by imaging for sandstones (Babaei & Joekar-Niasar, 2016 and reference therein). In this study, the radii of throats ($r_{th}$) were calculated using the radii of the two neighbouring pore bodies ($r_i$ and $r_j$) [Joekar-Niasar et al., 2008]:

$$r_{th} = r_i r_j (r_i^{1/0.9} + r_j^{1/0.9})^{-0.9} \qquad \text{Eq. 16}$$

<table><caption>**Table 1.** List of generated uncorrelated and correlated pore networks, their geometry, pore radii distribution correlation length and absolute permeability. The coordination number of all generated networks is 5.3. Flow direction is along the x-direction.</caption><thead><tr><th rowspan="2">Pore network</th><th rowspan="2">Number of pore bodies</th><th colspan="3">Domain size (µm)</th><th colspan="3">Dimensionless correlation lengths</th><th rowspan="2">Correlation type</th><th rowspan="2">Absolute permeability (mD)</th></tr><tr><th>L<sub>x</sub></th><th>L<sub>y</sub></th><th>L<sub>z</sub></th><th>x-direction (λ<sub>x</sub>/L<sub>x</sub>)</th><th>y-direction (λ<sub>y</sub>/L<sub>y</sub>)</th><th>z-direction (λ<sub>z</sub>/L<sub>z</sub>)</th></tr></thead><tbody><tr><td>UNC1</td><td>22204</td><td>7,000</td><td>4,600</td><td>1150</td><td>0.005</td><td>0.005</td><td>0.005</td><td>uncorrelated</td><td>95.3</td></tr><tr><td>UNC2</td><td>22166</td><td>7,000</td><td>4,600</td><td>1150</td><td>0.005</td><td>0.005</td><td>0.005</td><td>uncorrelated</td><td>96.8</td></tr><tr><td>COR1</td><td>22207</td><td>7,000</td><td>4,600</td><td>1150</td><td>0.05</td><td>0.05</td><td>0.05</td><td>correlated</td><td>86.7</td></tr><tr><td>COR2</td><td>22167</td><td>7,000</td><td>4,600</td><td>1150</td><td>0.05</td><td>0.05</td><td>0.05</td><td>correlated</td><td>107.6</td></tr></tbody></table>

<footnote>† The median coordination number for pores in the Berea sandstone is 6, compared to the median coordination number of 2 for the Fontainebleau sandstone (Thomson et al., 2020).</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The pore size distribution of pore networks of (a) UNC1, (b) UNC2, (c) COR1, and (d) COR2.</figcaption></figure>

### 2.2.2. Properties of fluid and bacteria

Toluene is selected as the NAPL phase in the pore network with a diffusion coefficient of $8.6 \times 10^{-6} \text{ cm}^2 \text{ sec}^{-1}$, solubility of $495.6 \text{ mg L}^{-1}$ in water, and density of $0.865 \text{ g cm}^{-3}$. In this work, it is assumed the pore network is water-wet, so the residual NAPL occupies the large pores. Most of the studies have shown that the residual non-wetting phase is trapped in the larger pore space (Powers, 1992, Chatzis et al., 1983, Wardlaw, 1982, Zhang & Guo, 2008). In this work 40% residual NAPL saturation was used for simulation, accordingly for this entrapped NAPL saturation/distribution, the largest pores are filled with toluene until the NAPL saturation reaches to 40%.

To simulate natural attenuation of toluene, it is assumed that toluene is biodegraded by *P. putida* as a bacterium which exists naturally in pore bodies and throats. There is a large range of variability in the parameter values of biodegradation model which are found in the literature. For instance, different values of $\mu_{max}$ measured for toluene biodegradation for different types of Pseudomonas are reported by Bordel et al., 2007, (0.78 h⁻¹), Oh et al., 1994, (1.56 h⁻¹), and Câmara et al., 2019, (1.72 h⁻¹).

In a similar way, the yield coefficient for biodegradation of toluene by Pseudomonas putida is 0.6 biomass mg/mg NAPL from Robledo-Ortiz et al., 2011, and 1.28 biomass mg/mg from Reardon et al., 2000. In this work, we used a biomass density as the total (wet) mass per volume. Classically, the dry weight represents between 20 and 60% of the respective wet mass of biomass depending on the amount
of intercellular water so that the yield coefficient (based on dry bio mass) should be scaled by a factor of 2 to 5. As a consequence, in this work, we have kept the value of 1.28 (wet) mg/mg for the yield coefficient (reported by Reardon et al., 2000), i.e., but using the value of 0.6 (dry biomass) mg/mg found in Robledo-Ortíz et al., 2011, with a conversion factor of 0.48 (ratio of dry mass over wet mass of bacterial cells).

In this work, the growth kinetic parameters (Haldane model) for biodegradation of toluene by *P. putida* together with other parameters used for bioremediation model are listed in **Table 2**. The growth kinetic parameters taken from Mathur and Majumder, 2010 (μ, Kᵢ, and Kₘ) are estimated in the range of 10-400 mg.L⁻¹ toluene concentration. Since the toluene concentration in the presence of bio-surfactant may increase up to 800-900 mg/L, it is assumed that these growth kinetic parameters can be extrapolated.

In this work, rhamnolipids is assumed as a produced bio-surfactant by *Pseudomonas* species/bacteria including *Pseudomonas putida* (Wittgens et al., 2011, Nanganuru & Korrapati, 2012, Meliani & Bensoltane, 2014, Kaskatepe & Yildiz, 2016). Rhamnolipid is the best-known glycolipid bio-surfactant with one or two molecules of β-hydroxydecanoic acid (Kaskatepe & Yildiz, 2016). In this work, the partitioning coefficient of toluene in the presence of rhamnolipids is approximately set to 1.5 (McCray et al., 2001). In addition, it is assumed that the initial thickness of biomass (which is attached to the rock surface) in the water-filled pores and throats is the 0.1% of their radii.

<table><caption>**Table 2.** List of parameters used in bioremediation model. * The yield coefficient found in Robledo-Ortíz et al. (2011) is expressed as a function of the dry mass of biomass. It was scaled by assuming a dry weight/wet weight ratio of 0.48 for bacterial cells in our simulations.</caption><thead><tr><th>Property</th><th></th><th>Value</th><th>Representative value from</th></tr></thead><tbody><tr><td>μ</td><td>maximum reaction rate (hr<sup>-1</sup>)</td><td>0.1722</td><td>(Mathur & Majumder, 2010)</td></tr><tr><td>K<sub>i</sub></td><td>inhibition coefficient (mg L<sup>-1</sup>)</td><td>380.62</td><td>(Mathur & Majumder, 2010)</td></tr><tr><td>K<sub>m</sub></td><td>half-saturation constant (mg L<sup>-1</sup>)</td><td>62.56</td><td>(Mathur & Majumder, 2010)</td></tr><tr><td>Y<sub>bio</sub></td><td>yield or stoichiometric coefficient</td><td>1.28</td><td>(Robledo-Ortíz et al., 2011)*</td></tr><tr><td>ρ<sub>bio</sub></td><td>bacterial density (g cm<sup>-3</sup>)</td><td>1</td><td>(Zhang & Bishop, 1994)</td></tr><tr><td>K<sub>ext</sub></td><td>extinction coefficient (h<sup>-1</sup>)</td><td>0.0086</td><td>5% of the maximum specific growth rate</td></tr><tr><td>K<sub>part</sub></td><td>partitioning coefficient</td><td>1.5</td><td>(McCray et al., 2001)</td></tr><tr><td>V<sub>colony</sub></td><td>volume of colony μm<sup>3</sup></td><td>1</td><td>Estimated</td></tr><tr><td>N<sub>colony,i</sub><sup>t=0</sup></td><td>initial number of colonies</td><td>10</td><td>estimated</td></tr></tbody></table>
## 3 – Results

To investigate the effect of pore scale heterogeneities on NAPL dissolution and bio-degradation in the process of natural attenuation for different pore network structures listed in **Table 1**, freshwater is flushed from the left-hand side of pore network in the x-direction. At field and laboratory conditions, the aqueous flow rate for groundwater is low (Imhoff et al., 2002, Soga et al., 2004, Farthing et al., 2012, Aminnaji et al., 2015). Therefore, a constant differential pressure (5 Pa) was applied through the pore network in the x-direction to achieve a low flow rate of ~1.3×10⁻⁵ cm³.min⁻¹ at initial condition. Varying heterogeneities and NAPL distribution in the pore networks and also NAPL dissolution and biofilm growth over time result in different flow rates during the simulation (as they modify available space for flow). Therefore, a dimensionless time that represents the pore volume injected (flow rate × time/pore volume) is used to compare the results and hereafter is referred to as PVI.

NAPL recovery (due to both NAPL dissolution and bio-degradation) as a function of dimensionless time is presented in **Figure 3** for different pore network structures with an initial NAPL saturation of 40%. The NAPL recovery due to only NAPL dissolution is also presented in the inset of **Figure 3** (it is assumed that there are not any bacteria for NAPL biodegradation). The results point out the significance of NAPL biodegradation when there is biofilm in the porous media. Firstly until 100 PVI, all non-bio cases (shown in blue) overlap each other showing that biofilm growth is more susceptible to the pore structure than non-bio only-dissolution cases. Secondly, a large discrepancy appears between the non-bio cases. This implies that, although the natural biodegradation process could be very slow, it has a significant effect compared to the only natural dissolution process. Chu et al., 2003, showed that there is a significant enhancement in the NAPL pool dissolution in the presence of biomass near to the water-NAPL interface when there is unlimited electron donor.

The results for bioremediation indicate that as NAPL is more distributed through the pore network (uncorrelated cases shown in green), the process of natural attenuation enhances NAPL removal. Two reasons are suggested as being the causes of the enhanced process of natural attenuation in the uncorrelated pore networks: (1) enhanced NAPL dissolution and (2) enhanced biodegradation. The NAPL concentration profiles shown in the inset of **Figure 3** indicate that while the dissolved NAPL in
the aqueous phase is distributed through the entire uncorrelated pore networks, the high concentration of NAPL is only observed in some areas in the correlated pore networks (i.e., NAPL is trapped as clusters in these areas). As NAPL is more distributed through the pore network, it provides more NAPL-water interface surface area which results in enhancement of dissolution.

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

**Figure 3.** NAPL recovery as a function of PVI for various pore network structures with and without bio implementation. Initial NAPL saturation is 40%. Insets: bio mass distribution through the networks at 60% NAPL recovery and dissolved NAPL concentration in the water phase at 50% NAPL recovery.

**Figure 4** shows the contribution of biodegradation (including biodegradation of dissolved NAPL and NAPL phase) and dissolved NAPL production in the NAPL recovery (remediation) as a function of time and pore volume injection for the correlated and uncorrelated pore networks. The results indicate the significant influence of biodegradation compared to the natural dissolution in NAPL remediation, as natural dissolution of NAPL into the water could be a very low process due to its low solubility and mass transfer coefficient. However, the contribution of NAPL biodegradation very much depends on the amount of biomass and its growth rate (through the parameters listed in **Table 2**).
<figure><img src="image_5.png" /><figcaption>**Figure 4.** NAPL recovery due to total biodegradation and production (NAPL which was flushed out from the system) as a function of time and PVI for (a) uncorrelated pore network, and (b) correlated pore network.</figcaption></figure>

Results are shown in **Figure 5** point out that the bio-growth in the uncorrelated pore networks was higher than in the correlated pore networks. The reason for this behaviour is that the dissolved NAPL is more distributed in the uncorrelated pore networks, providing the nutrition of bacteria for growth through the entire uncorrelated pore networks. For example, as shown in the inset of **Figure 3** for biomass distribution, while the map shows high biomass in some areas of correlated pore networks, it has grown through the entire uncorrelated pore networks. The biomass distribution agrees with dissolved NAPL distribution through the pore networks. The results of bio mass distribution agree with the Darcy numerical simulation for PCE biodegradation conducted by Chu et al., 2003, i.e., more biomass is expected to be produced in the vicinity of NAPL area. A higher biomass growth suggests more NAPL remediation by biodegradation. As shown in **Figure 5**, the mass of bio and the number of colonies for uncorrelated pore networks were higher than correlated pore networks at a specific PVI. Therefore, the NAPL biodegradation rate in the uncorrelated pore networks was higher than that of the correlated pore networks. **Figure 5** also shows that the rate of increase of biomass was high at the initial time of simulation and gradually decreased and finally the mass of bio started to decrease. This is because of the fact that as NAPL was removed over time, the amount of dissolved NAPL available for bio-growth (as the nutrient for bacteria) decreased.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Bio-growth: (a) mass of biofilm relative to its initial mass (b) number of colonies relative to its initial number as a function of PVI for different pore network structures with an initial NAPL saturation of 40%. The initial biomass and colonies refer to the amount of biomass and number of colonies respectively at the time PVI = 0 when the flow is started to inject.</figcaption></figure>

Additional insight can be gained from these pore-scale simulations by upscaling the averaged properties and identifying the heterogeneities effect on these properties. **Figure 6 (a)** and **Figure 6 (b)** show the average specific growth rate as a function of PVI and average resident concentration of NAPL in the aqueous phase respectively. This growth rate includes biodegradation in both aqueous phase (dissolved NAPL) and NAPL phase. It clearly indicates that the average specific growth rate for a specific average resident concentration of NAPL for uncorrelated pore network is higher than the correlated pore networks, albeit they are both lower than intrinsic specific growth rate which is defined by Eq. 5 with the value of parameters listed in **Table 2**. The reason for the deviation of average specific growth rate from intrinsic growth rate is the difference of NAPL concentration of the pore bodies and pore throats in the network. As NAPL is more distributed in the uncorrelated pore networks, the dissolved NAPL is dispersed in the aqueous phase through the entire network which results in a reduction in the variation of dispersion of concentration (the concentrations tend to be close to the mean). In contrast to the uncorrelated pore networks, the dispersivity of concentration in the correlated pore network is high (as shown in the inset of **Figure 3**) which results in more deviation of average specific growth rate from intrinsic biodegradation rate.

To put the above observation into perspective, as shown in **Figure 6 (b)**, the maximum specific growth rate (based on Haldane model) is at the concentration of ~150 mg L⁻¹, so the concentration below and above this threshold value results in the reduction of biodegradation rate. Therefore, for example for
the correlated pore networks, as the most of the NAPL concentrations in the pores are either close to solubility of toluene (495.6 mg L⁻¹) (for those pore bodies close to the NAPL-filled pores) or close to zero (for those pore bodies far from the NAPL-filled pores), the average specific growth rate is lower than the maximum and also lower than calculated intrinsic growth rate at the average resident concentration of NAPL in the pore network. This behaviour agrees with the results of pore scale simulation of NAPL biodegradation reported by Benioug et al. (2019). They pointed out that in areas of porous media which are too close or too far from the NAPL source (e.g. the NAPL distribution observed in correlated pore network structure), the bacterial mortality increases which results in reduction in total biomass growth rate.

The deviation of average specific growth rate for differently heterogeneous pore networks from the intrinsic biodegradation rate, as shown in **Figure 6** (b), points out that the parameters used in Eq. 5 (maximum reaction rate, inhibition coefficient, and half-saturation constant) should be upscaled from pore scale to be used for various heterogeneities of the continuum scale (Darcy scale) modelling. The necessity for scaling of these parameters and mass transfer coefficient has been demonstrated in other studies. Chu et al., 2004, pointed out that the inhibition coefficient which is measured in the laboratory could not be used in macroscopic model and it is challenging to be used for the complex NAPL distribution in the real field. On the other hand, Chu et al., 2007, showed that lump mass transfer coefficient depends on the scale and biological reaction. Therefore, the three terms for dissolved-NAPL biodegradation in the pore network modelling (last three terms in Eq. 7) could be lumped into one macroscopic term ($R(C_{ave}^t)\rho_{bio}V_{bio}^t/Y_{bio}$) to be used in the continuum scale single-porosity model (Aminnaji et al., 2019). **Table 3** lists the upscaled coefficients of Eq. 5 for different pore network structures (i.e., Eq.5 is fitted to data presented in **Figure 6** (b) to calculate the upscaled coefficients).

The results indicate that when the heterogeneity of the pore network is correlated, while upscaled maximum reaction rate ($\mu$) decreases, inhibition coefficient ($K_i$) and half-saturation constant ($K_m$) increase.
<table><caption>**Table 3.** The upscaled coefficients of Eq. 5 (maximum reaction rate, inhibition coefficient, and half-saturation constant) for various heterogeneities of pore networks.</caption><thead><tr><th rowspan="2">Coefficient</th><th colspan="4">Upscaled</th><th rowspan="2">Intrinsic value</th></tr><tr><th>UNC1</th><th>UNC2</th><th>COR1</th><th>COR2</th></tr></thead><tbody><tr><td>μ</td><td>0.1231</td><td>0.1300</td><td>0.0976</td><td>0.1001</td><td>0.1722</td></tr><tr><td>K<sub>m</sub></td><td>82.50</td><td>88.29</td><td>91.26</td><td>99.26</td><td>62.56</td></tr><tr><td>K<sub>i</sub></td><td>879.72</td><td>757.70</td><td>3103.33</td><td>2678.77</td><td>380.62</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Average specific growth rate over the entire pore network (with an initial NAPL saturation of 40%) as a function of (a) PVI, and (b) average resident concentration of NAPL in the aqueous phase.</figcaption></figure>

As biomass increases due to biofilm growth in the pores and throats, the porosity and permeability of pore networks change as the biofilm is considered as an impermeable/solid phase which adheres to the surface resulting in permeability and porosity reduction. **Figure 7** (a) shows the change of effective permeability relative to the initial absolute permeability due to the biofilm growth and NAPL remediation as a function of PVI for different pore network structures. The inset of **Figure 7** (a) shows the change of effective permeability relative to the initial absolute permeability as a function of PVI in the case if the biofilm does not clog the porous media and biofilm behaves as a highly permeable phase.

In this case, the biofilm growth does not affect the effective permeability and only NAPL dissolution results in the increase of effective permeability. These figures clearly show that any reduction in effective permeability of the pore network is related to the increase of biomass and any increase of effective permeability is associated with either NAPL removal/remediation or bio-extinction (as shown in **Figure 5**, biomass started to decrease at the end of the process).
Interestingly, the effective permeability started to increase when the porosity was still decreasing (e.g., 80 PVI for the UNC1 and UNC2). The reason for this behaviour is that while the porosity decreased due to biofilm growth (reducing the solid space), NAPL was dissolved and biodegraded resulting in increase of effective permeability (without change in the solid space). This interplay of reduction in porosity and increase of effective permeability is confirmed by the inset of Figure 7 (a). The figure shows the effective permeability (if biofilm assumed as a highly permeable phase) was non-decreasing and it significantly increased after ~70 PVI, indicating increase of NAPL biodegradation.

<figure><img src="image_8.png" /><figcaption>Figure 7. (a) Permeability with considering of biofilm (assumed as an impermeable/solid phase) relative to the initial absolute permeability, inset: permeability with considering the biofilm as a highly permeable phase relative to initial absolute permeability, and (b) porosity relative to the initial porosity as a function of PVI for different pore network structures. The initial permeability and porosity refer to the permeability and porosity at the time PVI=0 when the flow is started to inject.</figcaption></figure>

Figure 8 shows the change of absolute permeability relative to the initial absolute permeability as a function biomass/porosity change. As biomass increased, the porosity reduced resulting in absolute permeability reduction. Kozeny-Carman is one the simple model for the permeability-porosity relationship which has the general form of Eq. 17 (Rodriguez et al., 2004).

$$K_{abs} = C^* \frac{\phi^{n+1}}{(1-\phi)^n} \qquad \text{Eq. 17}$$

where n and C* are empirical parameters. The absolute permeability of pore network due to the change of porosity/biomass is predicted using Eq. 17 as shown in Figure 8 (assuming n = 2, and C* was calculated using the initial absolute permeability and porosity of pore network as listed in Table 1). As shown in Figure 8, the range of changes in porosity is small, so the Kozeny-Carman equation for
permeability prediction due to the biofilm growth could not be re-examined in this way. However, the highlighted area in **Figure 8** indicates that the Kozeny-Carman equation could not be applicable in the case of the permeability changes due to biofilm growth/extinction as with decrease of porosity the deviation from the Kozeny-Carman equation increases. The highlighted area shows that the absolute permeability increased after the reduction (indicating biofilm extinction). This highlighted area illustrates a hysteresis for the absolute permeability due to biofilm extinction which cannot be captured in Kozeny-Carman equation by a single set of parameters. Nevertheless, Chu et al., 2003, showed that biologging can significantly reduce the hydraulic conductivity and Chu et al., 2004, pointed out that the hydraulic conductivity of porous media could be described as function of biomass. However, the distribution and morphology of the bacteria have a significant impact on the relationship between conductivity and biomass (Dupin & McCarty, 2000, Dupin et al., 2001).

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Main figure: absolute permeability as a function of porosity, where biofilm is assumed as an impermeable /solid phase which adheres to the pore and throat surfaces, resulting in porosity reduction. Inset: absolute permeability as a function of biomass. KC stands for Kozeny-Carman. The initial permeability and porosity refer to the permeability and porosity at the time PVI=0 when the flow is started to inject.</figcaption></figure>

The results also indicate that the permeability of uncorrelated pore networks was reduced more than that of correlated pore networks. This behaviour follows the trend of the amount of grown bacteria in
the pore networks as shown in **Figure 5**. This also corresponds to the volume change shown in **Figure 9**. **Figure 9** shows the histogram of the percentage of pore/throat volume decrement changes (at the end of the simulation when 99% of NAPL was recovered) relative to their initial volumes for different pore network structures. It clearly points out that the number of pore/throats which underwent more volume reduction in the correlated pore networks are less than uncorrelated pore networks. For example, while ~56% of pores have a reduction of their volumes by <5% in the correlated pore networks, the reduction happened for less than ~15% of pores in the uncorrelated pore network. This indicates that the bio-growth was higher in the uncorrelated pore networks which results in more reduction in their porous space.

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Histogram of the percentage of volume change reduction (at the end of the simulation when 99% of NAPL was recovered) relative to the initial volume for (a) pore bodies, and (b) pore throats for different pore network structures with an initial NAPL saturation of 40%.</figcaption></figure>

## 4 – Discussions about assumptions, conclusions and future works

In this work, natural bioremediation and dissolution of NAPL were investigated using three-dimensional pore network modelling. Firstly an idealized condition was assumed where the pores were fully NAPL-filled and a large excess of electron acceptor (e.g. Oxygen) is available, leading to overestimation of biodegradation in our modelling. However, we ignored corner flow (allowing for flow of water along the corners of pores and throats containing trapped NAPL ganglia), and by doing so the biodegradation rate was adversely affected. In presence of these limiting assumptions, we found that structures of pore network (in the form of different pore radii correlation lengths) had a significant
effect on the natural attenuation process. The results showed enhancement of natural attenuation process
in the uncorrelated pore networks due to the NAPL being well-distributed through the pore network.
Two processes were suggested as being the root of the enhanced process of natural attenuation in the
uncorrelated pore networks —that due to simplifying assumptions behave similar to batch-type
biodegradation at idealized condition: (1) enhancement in NAPL dissolution and (2) enhancement in
biodegradation. While NAPL biodegradation in the uncorrelated pore network was higher than
correlated pore network, more reduction was observed in the physical properties of uncorrelated pore
network (e.g., permeability and porosity) as more biomass was formed.

Results also indicated that average specific bio-growth rate deviated from the intrinsic specific growth
rate in various pore network structures and decreased with pore radii correlation length. Therefore, in
the real field application, the parameters used for specific bio-growth rate (maximum reaction rate,
inhibition coefficient, and half-saturation constant) should be upscaled to be used in the continuum scale
(Darcy scale) model with respect to the various microscopic heterogeneities of the porous medium.

While in this work we showed how biofilm growth and NAPL dissolution were governed by
microscopic heterogeneity of pore network structures (mainly through the initial distribution of NAPL
in pores), this dependence will increase when more physically sound approaches for NAPL dissolution
such as corner flow diffusion or more realistic limitation for biodegradation (not all NAPL fully
available to be remediated) are considered. In such cases, the difference in the distribution pattern of
NAPL ganglia will be more important as the throats can provide pathways for water flow to the pores
that are filled with NAPL in clusters (correlated heterogeneity) rather than scattered blobs (uncorrelated
heterogeneity), and also the limitation in biodegradation rate will lead to different interplay with NAPL
dissolution for different structures. Therefore, the future recommendation for a better prediction of pore
scale biodegradation process is the use of dynamic pore network modelling and more realistic
biodegradation models.

Additionally the technical accuracy of the proposed model (or any future proposed models for water-NAPL-biofilm interplay) should be examined experimentally and more in-depth by using x-ray microtomography imaging of the process in microscopic set-ups specially designed in such a way that
heterogeneity at pore scale is representative. This will lead to more insights into how fluid flow, transport phenomena and biofilms interact within microscopically heterogeneous structures. The current state-of-the-art imaging of biofilm growth simultaneous to NAPL dissolution faces multiple challenges as biofilms are difficult to detect due to the lack of biofilm-specific biomarkers and methods for non-destructive imaging (Xu et al., 2020).

## Acknowledgement

The authors sincerely thank the anonymous reviewers for their constructive comments during the revision process.

## References

Agaoglu, B., Scheytt, T., & Copty, N. K. (2016). Impact of NAPL architecture on interphase mass transfer: A pore network study. *Advances in Water Resources, 95*, 138–151.

Aghaei, A., & Piri, M. (2015). Direct pore-to-core up-scaling of displacement processes: Dynamic pore network modeling and experimentation. *Journal of Hydrology, 522*, 488–509.

Aminnaji, M., Fazeli, H., Bahramian, A., Gerami, S., & Ghojavand, H. (2015). Wettability alteration of reservoir rocks from liquid wetting to gas wetting using nanofluid. *Transport in Porous Media, 109*(1), 201–216.

Aminnaji, M., Rabbani, A., Niasar, V. J., & Babaei, M. (2019). Effects of pore-scale heterogeneity on macroscopic NAPL dissolution efficiency: A two-scale numerical simulation study. *Water Resources Research*.

Atlas, R. M. (1981). Microbial degradation of petroleum hydrocarbons: an environmental perspective. *Microbiological Reviews, 45*(1), 180.

Babaei, M., & Joekar-Niasar, V. (2016). A transport phase diagram for pore-level correlated porous media. *Advances in Water Resources, 92*, 23–29.

Bahar, T., Golfier, F., Oltéan, C., & Benioug, M. (2016). An Upscaled Model for Bio-Enhanced NAPL Dissolution in Porous Media. *Transport in Porous Media, 113*(3), 653–693. https://doi.org/10.1007/s11242-016-0718-7

Banat, I. M. (1995). Biosurfactants production and possible uses in microbial enhanced oil recovery and oil pollution remediation: a review. *Bioresource Technology, 51*(1), 1–12.

Baveye, P., & Valocchi, A. (1989). An evaluation of mathematical models of the transport of biologically reacting solutes in saturated soils and aquifers. *Water Resources Research, 25*(6), 1413–1421.

Becker, J. G., & Seagren, E. A. (2009). Modeling the effects of microbial competition and hydrodynamics on the dissolution and detoxification of dense nonaqueous phase liquid contaminants. *Environmental Science & Technology, 43*(3), 870–877.

Benioug, M., Golfier, F., Oltéan, C., Buès, M. A., Bahar, T., & Cuny, J. (2017). An immersed boundary-lattice Boltzmann model for biofilm growth in porous media. *Advances in Water Resources, 107*, 65–82.

Benioug, M., Golfier, F., Fischer, P., Oltean, C., Buès, M. A., & Yang, X. (2019). Interaction between biofilm growth and NAPL remediation: A pore-scale study. *Advances in Water Resources, 125*, 82–97.

Bhatnagar, S., & Kumari, R. (2013). Bioremediation: a sustainable tool for environmental management--a review.
Annual Research & Review in Biology, 974-993.

Bhatnagar, S., & Kumari, R. (2015). Bioremediation: a sustainable tool for environmental management - A review. *Nature Environment and Pollution Technology, 14*(1), 181-190.

Bijeljic, B., Muggeridge, A. H., & Blunt, M. J. (2004). Pore-scale modeling of longitudinal dispersion. *Water Resources Research, 40*(11).

De Blanc, P. C., McKinney, D. C., & Speitel Jr, G. E. (1996). Modeling subsurface biodegradation of non-aqueous phase liquids. In *Advances in porous media* (Vol. 3, pp. 1-86). Elsevier.

Bordel, S., Muñoz, R., Dv'liaz, L. F., & Villaverde, S. (2007). New insights on toluene biodegradation by Pseudomonas putida F1: influence of pollutant concentration and excreted metabolites. *Applied Microbiology and Biotechnology*, 74(4), 857-866.

Brauner, J. S., & Widdowson, M. A. (2001). Numerical simulation of a natural attenuation experiment with a petroleum hydrocarbon NAPL source. *Groundwater*, 39(6), 939-952.

Brown, L. R. (2010). Microbial enhanced oil recovery (MEOR). *Current Opinion in Microbiology*, 13(3), 316-320.

Câmara, J. M. D. A., Sousa, M., & B Neto, E. L. (2019). KINETICS OF THE BIODEGRADATION OF MONOAROMATICS BY Pseudomonas aeruginosa. *Brazilian Journal of Chemical Engineering*, 36(1), 65-72.

Chambon, J. C., Broholm, M. M., Binning, P. J., & Bjerg, P. L. (2010). Modeling multi-component transport and enhanced anaerobic dechlorination processes in a single fracture--clay matrix system. *Journal of Contaminant Hydrology*, 112(1-4), 77-90.

Chatzis, I., Morrow, N. R., & Lim, H. T. (1983). Magnitude and detailed structure of residual oil saturation. *Society of Petroleum Engineers Journal*, 23(02), 311-326.

Chrysikopoulos, C. V., & Vogler, E. T. (2006). Acoustically enhanced ganglia dissolution and mobilization in a monolayer of glass beads. *Transport in Porous Media*, 64(1), 103-121.

Chu, M., Kitanidis, P. K., & McCarty, P. L. (2003). Effects of biomass accumulation on microbially enhanced dissolution of a PCE pool: a numerical simulation. *Journal of Contaminant Hydrology*, 65(1-2), 79-100.

Chu, M., Kitanidis, P. K., & McCarty, P. L. (2004). Possible factors controlling the effectiveness of bioenhanced dissolution of non-aqueous phase tetrachloroethene. *Advances in Water Resources*, 27(6), 601-615.

Chu, M., Kitanidis, P. K., & McCarty, P. L. (2007). Dependence of lumped mass transfer coefficient on scale and reactions kinetics for biologically enhanced NAPL dissolution. *Advances in Water Resources*, 30(6-7), 1618-1629.

Clement, T., Gautam, T., Lee, K., Truex, M., & Davis, G. (2004). Modeling of DNAPL-dissolution, rate-limited sorption and biodegradation reactions in groundwater systems. *Bioremediation Journal*, 8(1-2), 47-64.

Costerton, J. W., Stewart, P. S., & Greenberg, E. P. (1999). Bacterial biofilms: a common cause of persistent infections. *Science*, 284(5418), 1318-1322.

Davit, Y., Debenest, G., Wood, B. D., & Quintard, M. (2010). Modeling non-equilibrium mass transport in biologically reactive porous media. *Advances in Water Resources*, 33(9), 1075-1093.

Davit, Y., Iltis, G., Debenest, G., VERAN-TISSOIRES, S., Wildenschild, D., Gérino, M., & Quintard, M. (2011). Imaging biofilm in porous media using X-ray computed microtomography. *Journal of Microscopy*, 242(1), 15-25.

Detwiler, R. L., Rajaram, H., & Glass, R. J. (2001). Nonaqueous-phase-liquid dissolution in variable-aperture fractures: Development of a depth-averaged computational model with comparison to a physical experiment. *Water Resources Research*, 37(12), 3115-3129.

Dillard, L. A., & Blunt, M. J. (2000). Development of a pore network simulation model to study nonaqueous phase liquid dissolution. *Water Resources Research*, 36(2), 439-454.

Van Dillewijn, P., Caballero, A., Paz, J. A., González-Pérez, M. M., Oliva, J. M., & Ramos, J. L. (2007).
Bioremediation of 2,4,6-trinitrotoluene under field conditions. *Environmental Science and Technology*, 41(4), 1378-1383. https://doi.org/10.1021/es062165z

Dupin, H. J., & McCarty, P. L. (2000). Impact of colony morphologies and disinfection on biological clogging in porous media. *Environmental Science & Technology*, 34(8), 1513-1520.

Dupin, H. J., Kitanidis, P. K., & McCarty, P. L. (2001). Simulations of two-dimensional modeling of biomass aggregate growth in network models. *Water Resources Research*, 37(12), 2981-2994.

Ezeuko, C. C., Sen, A., Grigoryan, A., & Gates, I. D. (2011). Pore-network modeling of biofilm evolution in porous media. *Biotechnology and Bioengineering*, 108(10), 2413-2423.

Farthing, M. W., Seyedabbasi, M. A., Imhoff, P. T., & Miller, C. T. (2012). Influence of porous media heterogeneity on nonaqueous phase liquid dissolution fingering and upscaled mass transfer. *Water Resources Research*, 48(8).

Fountain, J. C., Starr, R. C., Middleton, T., Beikirch, M., Taylor, C., & Hodge, D. (1996). A controlled field test of surfactant-enhanced aquifer remediation. *Groundwater*, 34(5), 910-916.

Golfier, F., Wood, B. D., Orgogozo, L., Quintard, M., & Buès, M. (2009). Biofilms in porous media: Development of macroscopic transport equations via volume averaging with closure for local mass equilibrium conditions. *Advances in Water Resources*, 32(3), 463-485. https://doi.org/10.1016/j.advwatres.2008.11.012

Gooddy, D. C., Bloomfield, J. P., Harrold, G., & Leharne, S. A. (2002). Towards a better understanding of tetrachloroethene entry pressure in the matrix of Permo-Triassic sandstones. *Journal of Contaminant Hydrology*, 59(3-4), 247-265.

Held, R. J., & Celia, M. A. (2001). Pore-scale modeling and upscaling of nonaqueous phase liquid mass transfer. *Water Resources Research*, 37(3), 539-549.

Heron, G., Carroll, S., & Nielsen, S. G. (2005). Full-scale removal of DNAPL constituents using steam-enhanced extraction and electrical resistance heating. *Groundwater Monitoring & Remediation*, 25(4), 92-107.

Huang, J., Christ, J. A., Goltz, M. N., & Demond, A. H. (2015). Modeling NAPL dissolution from pendular rings in idealized porous media. *Water Resources Research*, 51(10), 8182-8197.

Huling, S. G., & Pivetz, B. E. (2006). *In-situ chemical oxidation*.

Imhoff, P. T., Farthing, M. W., Gleyzer, S. N., & Miller, C. T. (2002). Evolving interface between clean and nonaqueous phase liquid (NAPL)-contaminated regions in two-dimensional porous media. *Water Resources Research*, 38(6), 21-29.

Javanbakht, G., & Goual, L. (2016). Mobilization and micellar solubilization of NAPL contaminants in aquifer rocks. *Journal of Contaminant Hydrology*, 185, 61-73.

Jia, C., Shing, K., & Yortsos, Y. C. (1999). Visualization and simulation of non-aqueous phase liquids solubilization in pore networks. *Journal of Contaminant Hydrology*, 35(4), 363-387.

Joekar-Niasar, V, & Hassanizadeh, S. M. (2012). Analysis of fundamentals of two-phase flow in porous media using dynamic pore-network models: A review. *Critical Reviews in Environmental Science and Technology*, 42(18), 1895-1976.

Joekar-Niasar, V, Hassanizadeh, S. M., & Dahle, H. K. (2010). Non-equilibrium effects in capillarity and interfacial area in two-phase flow: dynamic pore-network modelling. *Journal of Fluid Mechanics*, 655, 38-71.

Joekar-Niasar, Vahid, & Hassanizadeh, S. M. (2011). Effect of fluids properties on non-equilibrium capillarity effects: Dynamic pore-network modeling. *International Journal of Multiphase Flow*, 37(2), 198-214.

Kaskatepe, B., & Yildiz, S. (2016). Rhamnolipid biosurfactants produced by Pseudomonas species. *Brazilian Archives of Biology and Technology*, 59.

Khasi, S., Ramezanzadeh, M., & Ghazanfari, M. H. (2020). Experimentally based pore network modeling of NAPL dissolution process in heterogeneous porous media. *Journal of Contaminant Hydrology*, 228, 103565.
Lari, K. S., Davis, G. B., Rayner, J. L., Bastow, T. P., & Puzon, G. J. (2019). Natural source zone depletion of LNAPL: A critical review supporting modelling approaches. *Water Research*.

Larue, A. E., Swider, P., Duru, P., Daviaud, D., Quintard, M., & Davit, Y. (2018). Quantitative 3D comparison of biofilm imaged by X-ray micro-tomography and two-photon laser scanning microscopy. *Journal of Microscopy*, **271**(3), 302–314.

Lazar, I., Petrisor, I. G., & Yen, T. F. (2007). Microbial enhanced oil recovery (MEOR). *Petroleum Science and Technology*, **25**(11), 1353–1366. https://doi.org/10.1080/10916460701287714

Manoli, G., Chambon, J. C., Bjerg, P. L., Scheutz, C., Binning, P. J., & Broholm, M. M. (2012). A remediation performance model for enhanced metabolic reductive dechlorination of chloroethenes in fractured clay till. *Journal of Contaminant Hydrology*, **131**(1–4), 64–78.

Mathur, A. K., & Majumder, C. B. (2010). Kinetics modelling of the biodegradation of benzene, toluene and phenol as single substrate and mixed substrate by using Pseudomonas putida. *Chemical and Biochemical Engineering Quarterly*, **24**(1), 101–109.

McCray, J. E., Bai, G., Maier, R. M., & Brusseau, M. L. (2001). Biosurfactant-enhanced solubilization of NAPL mixtures. *Journal of Contaminant Hydrology*, **48**(1–2), 45–68.

Meliani, A., & Bensoltane, A. (2014). Enhancement of hydrocarbons degradation by use of Pseudomonas biosurfactants and biofilms. *Journal of Petroleum & Environmental Biotechnology*, **5**(1), 1.

Mostaghimi, P., Bijeljic, B., Blunt, M., & others. (2012). Simulation of flow and dispersion on pore-space images. *SPE Journal*, **17**(04), 1–131.

Mulder, H., Breure, A. M., Van Honschooten, D., Grotenhuis, J. T. C., Van Andel, J. G., & Rulkens, W. H. (1998). Effect of biofilm formation by Pseudomonas 8909N on the bioavailability of solid naphthalene. *Applied Microbiology and Biotechnology*, **50**(2), 277–283.

Mulligan, C. N., Yong, R. N., & Gibbs, B. F. (2001). Surfactant-enhanced remediation of contaminated soil: a review. *Engineering Geology*, **60**(1–4), 371–380.

Nanganuru, H. Y., & Korrapati, N. (2012). Studies on biosorption of cadmium by Pseudomonas putida. *Int J Eng Res Appl*, **2**(3), 2217–2219.

Oh, Y.-S., Shareefdeen, Z., Baltzis, B. C., & Bartha, R. (1994). Interactions between benzene, toluene, and p-xylene (BTX) during their biodegradation. *Biotechnology and Bioengineering*, **44**(4), 533–538.

Okaygun, M. S., Green, L. A., & Akgerman, A. (1992). Effect of Consecutive Pulsing of an Inhibitory Substrate on Biodegradation Kinetics. *Environmental Science and Technology*, **26**(9), 1746–1752. https://doi.org/10.1021/es00033a005

Øren, P.-E., & Bakke, S. (2003). Reconstruction of Berea sandstone and pore-scale modelling of wettability effects. *Journal of Petroleum Science and Engineering*, **39**(3–4), 177–199.

Orgogozo, L., Golfier, F., Buès, M., & Quintard, M. (2010). Upscaling of transport processes in porous media with biofilms in non-equilibrium conditions. *Advances in Water Resources*, **33**(5), 585–600.

Orgogozo, L., Golfier, F., Buès, M. A., Quintard, M., & Koné, T. (2013). A dual-porosity theory for solute transport in biofilm-coated porous media. *Advances in Water Resources*, **62**, 266–279. https://doi.org/10.1016/j.advwatres.2013.09.011

Ostvar, S., Iltis, G., Davit, Y., Schlüter, S., Andersson, L., Wood, B. D., & Wildenschild, D. (2018). Investigating the influence of flow rate on biofilm growth in three dimensions using microimaging. *Advances in Water Resources*, **117**, 1–13.

Peszynska, M., Trykozko, A., Iltis, G., Schlueter, S., & Wildenschild, D. (2016). Biofilm growth in porous media: experiments, computational modeling at the porescale, and upscaling. *Advances in Water Resources*, **95**, 288–301.

Pintelon, T. R. R., Picioreanu, C., van Loosdrecht, M. C. M., & Johns, M. L. (2012). The effect of biofilm permeability on bio-clogging of porous media. *Biotechnology and Bioengineering*, **109**(4), 1031–1042.

Pizarro, G., Griffeath, D., & Noguera, D. R. (2001). Quantitative cellular automaton model for biofilms. *Journal*
of Environmental Engineering, 127(9), 782-789.

Powell, K. L., Taylor, R. G., Cronin, A. A., Barrett, M. H., Pedley, S., Sellwood, J., et al. (2003). Microbial contamination of two urban sandstone aquifers in the UK. *Water Research*, 37(2), 339-352.

Powers, S. E. (1992). *Dissolution of nonaqueous phase liquids in saturated subsurface systems*.

Qin, C.-Z., & Hassanizadeh, S. M. (2015). Pore-network modeling of solute transport and biofilm growth in porous media. *Transport in Porous Media*, 110(3), 345-367.

Reardon, K. F., Mosteller, D. C., & Bull Rogers, J. D. (2000). Biodegradation kinetics of benzene, toluene, and phenol as single and mixed substrates for Pseudomonas putida F1. *Biotechnology and Bioengineering*, 69(4), 385-400.

Robledo-Ort\'iz, J. R., Ram\'irez-Arreola, D. E., Pérez-Fonseca, A. A., Gómez, C., González-Reynoso, O., Ramos-Quirarte, J., & González-Núñez, R. (2011). Benzene, toluene, and o-xylene degradation by free and immobilized P. putida F1 of postconsumer agave-fiber/polymer foamed composites. *International Biodeterioration & Biodegradation*, 65(3), 539-546.

Rodriguez, E., Giacomelli, F., & Vazquez, A. (2004). Permeability-porosity relationship in RTM for different fiberglass and natural reinforcements. *Journal of Composite Materials*, 38(3), 259-268.

Sahloul, N. A., Ioannidis, M. A., & Chatzis, I. (2002). Dissolution of residual non-aqueous phase liquids in porous media: pore-scale mechanisms and mass transfer rates. *Advances in Water Resources*, 25(1), 33-49.

Saravanan, P., Pakshirajan, K., & Saha, P. (2008). Growth kinetics of an indigenous mixed microbial consortium during phenol degradation in a batch reactor. *Bioresource Technology*, 99(1), 205-209.
https://doi.org/10.1016/j.biortech.2006.11.045

Sarikurt, D. A., Gokdemir, C., & Copty, N. K. (2017). Sherwood correlation for dissolution of pooled NAPL in porous media. *Journal of Contaminant Hydrology*, 206, 67-74.

Schmidt, S. I., Kreft, J.-U., Mackay, R., Picioreanu, C., & Thullner, M. (2018). Elucidating the impact of micro-scale heterogeneous bacterial distribution on biodegradation. *Advances in Water Resources*, 116, 67-76.

von der Schulenburg, D. A. G., Pintelon, T. R. R., Picioreanu, C., Van Loosdrecht, M. C. M., & Johns, M. L. (2009). Three-dimensional simulations of biofilm growth in porous media. *AIChE Journal*, 55(2), 494-504.

Singh, R., Paul, D., & Jain, R. K. (2006). Biofilms: implications in bioremediation. *Trends in Microbiology*, 14(9), 389-397.

Sivasankar, P., & Kumar, G. S. (2014). Numerical modelling of enhanced oil recovery by microbial flooding under non-isothermal conditions. *Journal of Petroleum Science and Engineering*, 124, 161-172.
https://doi.org/10.1016/j.petrol.2014.10.008

Smith, W. P. J., Davit, Y., Osborne, J. M., Kim, W., Foster, K. R., & Pitt-Francis, J. M. (2017). Cell morphology drives spatial patterning in microbial communities. *Proceedings of the National Academy of Sciences*, 114(3), E280-E286.

Soga, K., Page, J. W. E., & Illangasekare, T. H. (2004). A review of NAPL source zone remediation efficiency and the mass flux approach. *Journal of Hazardous Materials*, 110(1-3), 13-27.

Soudmand-asli, A., Ayatollahi, S. S., Mohabatkar, H., Zareie, M., & Shariatpanahi, S. F. (2007). The in situ microbial enhanced oil recovery in fractured porous media. *Journal of Petroleum Science and Engineering*, 58(1-2), 161-172. https://doi.org/10.1016/j.petrol.2006.12.004

Suchomel, B. J., Chen, B. M., & Allen, M. B. (1998). Network model of flow, transport and biofilm effects in porous media. *Transport in Porous Media*, 30(1), 1-23.

Tang, Y., Valocchi, A. J., Werth, C. J., & Liu, H. (2013). An improved pore-scale biofilm model and comparison with a microfluidic flow cell experiment. *Water Resources Research*, 49(12), 8370-8382.

Thomson, P. R., Ellis, R., Chiarella, D., & Hier-Majumder, S. (2020). Microstructural Analysis From X-Ray CT Images of the Brae Formation Sandstone, North Sea. Front. *Earth Sci.*, 8, 246.

Thullner, M., & Baveye, P. (2008). Computational pore network modeling of the influence of biofilm permeability
on bioclogging in porous media. *Biotechnology and Bioengineering*, **99**(6), 1337–1351.

Thullner, M., Zeyer, J., & Kinzelbach, W. (2002). Influence of microbial growth on hydraulic properties of pore networks. *Transport in Porous Media*, **49**(1), 99–122.

Vayenas, D. V., Michalopoulou, E., Constantinides, G. N., Pavlou, S., & Payatakes, A. C. (2002). Visualization experiments of biodegradation in porous media and calculation of the biodegradation rate. *Advances in Water Resources*, **25**(2), 203–219.

Wardlaw, N. C. (1982). The effects of geometry, wettability, viscosity and interfacial tension on trapping in single pore-throat pairs. *Journal of Canadian Petroleum Technology*, **21**(03).

Wealthall, G. P., Steele, A., Bloomfield, J. P., Moss, R. H., & Lerner, D. N. (2001). Sediment filled fractures in the Permo-Triassic sandstones of the Cheshire Basin: observations and implications for pollutant transport. *Journal of Contaminant Hydrology*, **50**(1–2), 41–51.

Wiedemeier, T. H., Rifai, H. S., Newell, C. J., & Wilson, J. T. (1999). *Natural attenuation of fuels and chlorinated solvents in the subsurface*. John Wiley & Sons.

Wittgens, A., Tiso, T., Arndt, T. T., Wenk, P., Hemmerich, J., Müller, C., et al. (2011). Growth independent rhamnolipid production from glucose using the non-pathogenic Pseudomonas putida KT2440. *Microbial Cell Factories*, **10**(1), 80.

Yang, X., Erickson, L. E., & Fan, L. T. (1995). A study of the dissolution rate-limited bioremediation of soils contaminated by residual hydrocarbons. *Journal of Hazardous Materials*, **41**(2–3), 299–313.

Zhang, L., & Guo, F. (2008). Micro-mechanisms of residual oil mobilization by viscoelastic fluids. *Petroleum Science*, **5**(1), 56–61.

Zhang, T. C., & Bishop, P. L. (1994). Density, porosity, and pore structure of biofilms. *Water Research*, **28**(11), 2267–2277.

Zhao, W., & Ioannidis, M. A. (2003). Pore network simulation of the dissolution of a single-component wetting nonaqueous phase liquid. *Water Resources Research*, **39**(10).

Zhao, Weishu, & Ioannidis, M. A. (2007a). Convective mass transfer across fluid interfaces in straight angular pores. *Transport in Porous Media*, **66**(3), 495–509.

Zhao, Weishu, & Ioannidis, M. A. (2007b). Effect of NAPL film stability on the dissolution of residual wetting NAPL in porous media: A pore-scale modeling study. *Advances in Water Resources*, **30**(2), 171–181.

Zhao, Weishu, & Ioannidis, M. A. (2011). Gas exsolution and flow during supersaturated water injection in porous media: I. Pore network modeling. *Advances in Water Resources*, **34**(1), 2–14.

Zhou, D., Dillard, L. A., & Blunt, M. J. (2000). A physically based model of dissolution of nonaqueous phase liquids in the saturated zone. *Transport in Porous Media*, **39**(2), 227–255.

Xu, Y., Dhaouadi, Y., Stoodley, P., & Ren, D. (2020). Sensing the unreachable: challenges and opportunities in biofilm detection. *Current Opinion in Biotechnology*, **64**, 79-84.