Multirate mass transfer and biofilm growth modeling in porous media

DOCTORAL THESIS Multirate mass transfer and biofilm growth modeling in porous media Author:Jingjing Wang Supervisors: Prof. Jesús Carrera Ramírez Prof. Maarten Willem Saaltink Dr. Cristina Valhondo A thesis submitted for the degree of Doctor of Philosophy in Geotechnical Engineering, Department of Civil and Environmental Engineering, Universitat Politecnica de Catalunya (UPC) June, 2021 Barcelona, Spain The author would like to acknowledge the financial support of the project NITREM (NITrogen REMoval from waste rock) initiated and funded by EIT (European Institute of Innovation and Technology), and the Water-JPI European Union project MARadentro. i Acknowledgements The PhD thesis has been a challenging experience for me. Not only my academic research in Geoscience, but also my life abroad. During the four years, I would like to express my thanks for the people mentioned below, without them I won’t be able to finish my thesis. First of all, I would like to express my heartfelt gratitude to my supervisor Jesús Carrera. I would like to thank Jesús for providing me the opportunity to achieve my PhD in Spain, for trusting me in the research topic, giving me unconditional support and encouragement during the whole process. Without his support, I won’t have the chance to achieve my PhD career. For me he is not only my supervisor, but also my friend and family. I would like to acknowledge the role of my advisor Maarten Saaltink. Thanks for his responsible guidance and participate in each discussion. For this, I respect and love him lots. I would also like to thank my co-supervisor Cristian Valhondo. It’s her independence, sense of responsibility and excellent ability that inspires me to perform well in the face of challenges. I would like to thank the people in Groundwater Hydrology Group (UPC, IDAEA-CSIC). They are excellent researchers and kind friends. I would like to specially thank the former and current members of the group that have helped me in one way or another. Thanks the people in office 1428, Alba, Quim, Rotman, Sonia, Nacho, Pedro. Thanks my friends Andrea and Tybaud for going through the PhD experience together, we witness the growth of each other, and we all achieve the PhD. Thanks the friends I have known in CSIC Alessandro, Vivien, Alec, Lazaro, Tomas, Silvia, Kevin, Miguel, Carolina, Robert, Francesco, Jordi, Julia, Berta, Yoar, Laura, S., Laura M., Lidia, Lucia, Ashkan, Ahmad, Max, Jiaqi, Joan, Sara, Carlos, Marco, Juan, Enric, Jordi, Josep. Thanks the people in UPC Michela, Laura V. A., Sonia, Oriol, Arnau, Guillem, Paula, Lurdes, Xavier, Daniel, Albert. I would especially like to thank Tere for her help in my living and working in Barcelona. I would also like to thank the experiences through my PhD. Thanks the training courses I have attended, including the 4th Cargèse summer school: flow and transport in porous and fractured media, the short course of applied reactive transport modelling (Bordeaux, France), and the short course of reactive transport in natural and engineered systems (Barcelona, Spain). Thanks the international conferences I have attended, including the EGU general assembly (Vienna, Austria), the Goldschmidt (Barcelona, Spain) and InterPore 2021 online (Edinburgh, United Kingdom). These experiences promote me to be more professional on my research topic. Finally, I would like to thank my family and friends in China for their love, encouragement and support. Thanks for every video and voice call that we have had on WeChat. A special thanks to myself for my persistence, brave and endeavor on the way of my PhD. These impressive experience promote my maturity and growth both in my life and research work. iii To my family v Abstract Reactive transport modeling is a methodological tool to study the coupled physical, chemical and biological processes in Earth system. It is complex not only because of the nature of the equations, but also because of the effects of the porous medium heterogeneity on reactive transport. This thesis aims to deepen the understanding of reactive transport processes in order to explain the biochemical degradation process in porous media, with special emphasis on the role of biofilm and its growth. First, we propose a general and efficient numerical solution of reactive transport in multicontinuum media using Multirate Mass Transfer (MRMT) approach. To overcome the non-linearity of the problem, induced by non-linear kinetics, we use the Newton-Raphson method to get the global solution. We solve the system of equations in block form, which allow us to reduce the unknowns to those of mobile zones and to, thus improving efficiency. The solution is validated by comparison with analytical solution for linear kinetics. The code is developed in Object-oriented way, which enables the code reusability and data polymorphism. Second, we investigate the conditions for chemical localization (i.e., the occurrence of reactions that would not be possible in single continuum media). To this end, we write the multicontinuum transport equations in dimensionless form to find that reactive transport in multicontinuum media is governed by three characteristic times: the distribution of residence times in immobile zones, and the characteristic reaction and transport times. To study the interplay between these three characteristic times, we simulate three chemical systems: conservative, single reaction and sequential reaction. Results demonstrate that reactions driven by species that result from previous reactions will localize in immobile zones whose residence time is comparable to reaction times. Furthermore, immobile zones with residence times much smaller than those for transport can be lumped together (assuming that very fast reactions are assumed in equilibrium), which greatly reduces computations. Third, we perform simulations of reactive transport incorporating biochemical reactions that not only oxidize organic carbon, but also produce biomass, thus causing biofilm growth. Biofilm growth is known to cause clogging (i.e., reduction of permeability), which has concentrated most research on the topic. But it also causes a significant change in the pore space geometry and connectivity, which leads to not only an overall increase in mean residence time in immobile regions, but also on its distribution. As discussed above, this is critical to (bio)chemical localization, especially considering that microbial mediated reactions tend to concentrate in biofilms. We propose a model for the evolution of residence time distribution in immobile zones in response to biofilm growth. We test this model by comparison with laboratory experiments extracted from the literature, where tracer tests have been performed at various stages of growth. Results show that the dynamic MRMT model is capable of reproducing the salient features of these experiments. vii Resumen El modelado de transporte reactivo es una herramienta metodológica para estudiar los procesos físicos, químicos y biológicos acoplados en el sistema terrestre. Es complejo no solo por la naturaleza de las ecuaciones, sino también por los efectos de la heterogeneidad del medio poroso sobre el transporte reactivo. Esta tesis tiene como objetivo profundizar en el conocimiento de los procesos de transporte reactivo para explicar el proceso de degradación bioquímica en medios porosos, con especial énfasis en el papel del biofilm y su crecimiento. En primer lugar, proponemos una solución numérica general y eficiente de transporte reactivo en medios multicontinuum utilizando el enfoque de Transferencia de Masa Multivelocidad (MRMT). Para superar la no linealidad del problema, inducida por una cinética no lineal, utilizamos el método de Newton-Raphson para obtener la solución global. Resolvemos el sistema de ecuaciones en forma de bloque, lo que nos permite reducir las incógnitas a las de zonas móviles y a, mejorando así la eficiencia. La solución se valida por comparación con la solución analítica para cinética lineal. El código se desarrolla de forma orientada a objetos, lo que permite la reutilización del código y el polimorfismo de los datos. En segundo lugar, investigamos las condiciones para la localización química (es decir, la ocurrencia de reacciones que no serían posibles en un medio continuo único). Con este fin, escribimos las ecuaciones de transporte multicontinuum en forma adimensional para encontrar que el transporte reactivo en medios multicontinuum se rige por tres tiempos característicos: la distribución de los tiempos de residencia en zonas inmóviles, y los tiempos de reacción y transporte característicos. Para estudiar la interacción entre estos tres tiempos característicos, simulamos tres sistemas químicos: conservador, reacción única y reacción secuencial. Los resultados demuestran que las reacciones impulsadas por especies que resultan de reacciones previas se localizarán en zonas inmóviles cuyo tiempo de residencia es comparable a los tiempos de reacción. Además, las zonas inmóviles con tiempos de residencia mucho más pequeños que los de transporte pueden agruparse (asumiendo que se asumen reacciones muy rápidas en equilibrio), lo que reduce en gran medida los cálculos. En tercer lugar, realizamos simulaciones de transporte reactivo incorporando reacciones bioquímicas que no solo oxidan el carbono orgánico, sino que también producen biomasa, lo que provoca el crecimiento de biopelículas. Se sabe que el crecimiento de biopelículas causa obstrucciones (es decir, reducción de la permeabilidad), lo que ha concentrado la mayor parte de la investigación sobre el tema. Pero también provoca un cambio significativo en la geometría del espacio poroso y la conectividad, lo que conduce no solo a un aumento general del tiempo medio de residencia en las regiones inmóviles, sino también a su distribución. Como se discutió anteriormente, esto es crítico para la localización (bio) química, especialmente considerando que las reacciones mediadas por microbios tienden a concentrarse en biopelículas. Proponemos un modelo para la evolución de la distribución del tiempo de residencia en zonas inmóviles en respuesta al crecimiento de biopelículas. Probamos este modelo comparándolo con experimentos de laboratorio extraídos de la literatura, donde se han realizado pruebas de trazadores en varias etapas de crecimiento. Los resultados muestran que el modelo MRMT dinámico es capaz de reproducir las características más destacadas de estos experimentos. CONTENTS xiv 3.4 Results ................................................................................................................................... 23 3.4.1 No reaction system ........................................................................................................ 23 3.4.2 Single reaction system .................................................................................................. 24 3.4.3 Sequential reaction system ............................................................................................ 26 3.5 Conclusions ........................................................................................................................... 28 4 A conceptual model for reactive transport with dynamic biofilm growth in multicontinuum media 31 4.1 Introduction ........................................................................................................................... 32 4.2 Conceptual models ................................................................................................................ 33 4.2.1 Reactive transport with MRMT .................................................................................... 33 4.2.2 Biofilm growth .............................................................................................................. 34 4.3 Simulation methods .............................................................................................................. 35 4.3.1 Reactive transport ......................................................................................................... 35 4.3.2 Model update ................................................................................................................ 35 4.3.3 Algorithms .................................................................................................................... 38 4.4 Applications .......................................................................................................................... 38 4.4.1 Experimental data of Kone (2014) ................................................................................ 38 4.4.2 Experimental data of Seifert (2007) .............................................................................. 40 4.5 Conclusions ........................................................................................................................... 43 5 Conclusions ................................................................................................................................... 45 A. Numerical discretization of governing equations ......................................................................... 47 B. Analytical solution of reactive transport in multicontinuum media for first-order kinetics in the Laplace Domain .................................................................................................................................... 49 C. Breakthrough curves of species A, B and C for sequential reactions ........................................... 51 D. The relationship between the mass transfer rates and the volume fraction of immobile zones .... 55 E. Publications and Conference presentations ................................................................................... 57 Publications in scientific journals ................................................................................................. 57 Presentations in Conferences ........................................................................................................ 57 Presentations in Seminar ............................................................................................................... 57 Bibliography ......................................................................................................................................... 59 xv List of figures Figure 2.1 (a) Illustration of MRMT model in porous media, white areas bounded by black curves represent mobile zone, the black areas represent solid matrix, and the gray areas represent immobile zones, dark gray corresponding to a lower possibility of visit, while light gray denotes a higher possibility to visit (Gouze et al., 2008b). (b) Numerical discretization of mobile and immobile zones, each circle (labelled m for mobile, and im for immobile) is a node. ...................................................... 4 Figure 2.2 Comparison between the results of the non-reactive transport model calculated by the proposed method and the method of Silva et al. (2009). ........................................................................ 9 Figure 2.3 Distribution of (left) components 𝑢𝑚 and 𝑢𝑖𝑚, (right) reaction rates 𝑟𝑚 and 𝑟𝑖𝑚 versus time at distance 𝑥=100.0m. The blue color represents state variables in the mobile zone, and the red color represents state variables in immobile zones. The solid lines are the analytical solutions and the circle dot lines are the proposed numerical solutions. .......................................................................... 10 Figure 2.4 Distribution of (left) components 𝑢𝑚 and 𝑢𝑖𝑚, (right) reaction rates 𝑟𝑚 and 𝑟𝑖𝑚 versus time at distance 𝑥=100.0m. The blue color represents state variables in the mobile zone, and the red color represents state variables in immobile zones. The solid lines are the analytical solutions and the circle dot lines are the proposed numerical solutions. .......................................................................... 10 Figure 2.5 CPU time comparison between proposed method and full DSA for different number of immobile zones. The chemical systems are identical for all cases that is an irreversible bimolecular reaction 𝐴+𝐵→𝐶, with a second-order kinetics 𝑟𝑘=𝜅𝑐𝐴𝑐𝐵. ......................................................... 11 Figure 3.1 Pdf of exchange rates, 𝛼 (Gamma distibution), and residence times, 𝜏𝛼, plotted in log-log scale for several values of 𝜏𝛼0 and 𝛽=12. Note that the maximum of latter is around 𝜏𝛼0 (the mode is 𝜏𝛼0/𝛽+1 and the expected value 𝜏𝛼0(1−𝛽)). Its large 𝜏𝛼 slope is 𝛽+1. ................................. 17 Figure 3.2 Steady state concentrations of species A (blue color), B (red color), and C (yellow color) as a function of residence time in immobile zones, 𝜏𝛼[T], assuming that mobile concentrations of B and C are zero. First rate constant, 𝑘1 [T-1], for A→B transformation equals 0.1 (solid lines), 1 (dot lines), and 10 (dash lines). The rate constant, 𝑘2 [T-1], for B→C reaction is 0.01. ......................................... 22 Figure 3.3 The ratios between the cumulative concentration in immobile zones and mobile zone for reactant species A (left) in immobile zones characterized by the residence times and (right) at different reaction times. The solid lines represent the analytical solutions, and the squares represent the numerical solutions. ............................................................................................................................................... 23 Figure 3.4 Mobile (top row) and average immobile (bottom row) concentrations versus time at 𝑥𝐷= 100𝐿𝑐 in response to a pulse input (left in arithmetic scale, center in log scale) and to a continuous injection (right) for transport in a medium with a distribution residence times in immobile zones with 𝛽=12, 𝜏𝛼0𝐷=102, 101, 100, and 10−1. The log-scale enhances the BTC tail (note vertical scale), which is only relevant in this case for 𝜏𝛼0𝐷=102. When exchange with the immobile zones is slow (i.e, large 𝜏𝛼0𝐷), the BTC peaks at the advective time in the mobile zone (𝑡𝐷=100, indicated by blue arrows). The peak occurs at the advective time for the full porosity (𝑡𝐷=400, red arrows) when axchange is fast, which also reduces the time it takes for immobile zones to equilibrate with inflow water. ..................................................................................................................................................... 24 Figure 3.5 Cumulative immobile scaled concentration profiles versus dimensionless distance and dimensionless residence times (pdf with 𝛽=12 and 𝜏𝛼0𝐷=10−1) at 𝑡𝐷=300. Note that immobile concentrations are identical to mobile concentrations for 𝜏𝛼𝐷<10−1. .......................... 24 Figure 3.6 Mobile (top row) and average immobile (bottom row) concentrations of species A (first and third columns) and B (second and fourth columns) versus time at 𝑥𝐷=100𝐿𝑐 in response to a pulse input (left two columns) and to a continuous injection (right columns) for transport in a medium with a distribution residence times in immobile zones with 𝛽=12, and 𝜏𝛼0𝐷=101, and reaction times are LIST OF FIGURES xvi uniformly distributed in immobile zones with 𝜏𝑟𝐷= ∞ , 104, 103 and 102. Species B is absent when no reaction occurs. ................................................................................................................................ 25 Figure 3.7 Cumulative mobile and averaged immobile concentration profiles versus dimensionless time and dimensionless distance for species A (left two columns) and B (right two columns). Species transport in a medium characterized by the immobile zones in which the residence time follows gamma distribution with 𝛽=12, 𝜏𝛼0𝐷=101. The reaction times are uniform in immobile zones with 𝜏𝑟𝐷= ∞ , 104, 103 and 102 from top to bottom. ........................................................................................... 25 Figure 3.8 Cumulative concentration profiles versus dimensionless distance and dimensionless residence times for species A (left) and B (right) at 𝑡𝐷=300 (above), and 1×104 (below) from top to bottom after continuous injection of A in a medium with immobile zones distributed with 𝛽=12, 𝜏𝛼0𝐷=101. Reaction A→B occurs with a characteristic time 𝜏𝑟𝐷=102. ..................................... 26 Figure 3.9 Breakthrough curves of species A, B and C in different cases for 𝜏𝑟1𝐷=𝜏𝑟2𝐷=103. The black lines represent results only with mobile zone. Case 1: 𝜏𝛼𝐷=100, case 2: 𝜏𝛼𝐷=100,101, case 3: 𝜏𝛼𝐷=100,101,102, case 4: 𝜏𝛼𝐷=100,101,102,103, case 5: 𝜏𝛼𝐷=100,101,102,103,104, case 6: 𝜏𝛼𝐷=100,101,102,103,104,105,106, case 7: 𝜏𝛼𝐷=102,103 and case 8: 𝜏𝛼𝐷=101,102,103. .............................................................................................................................................................. 27 Figure 3.10 Cumulative concentration profiles versus dimensionless distance and dimensionless residence times for species (left) A, (middle) B and (right) C for 𝜏𝑟1𝐷=𝜏𝑟2𝐷=103 at 𝑡𝐷=50, 200, 500, 2000 and 1×104 from top to bottom. Species transport in a medium characterized by the immobile zones in which the residence time follows gamma distribution with 𝛽=12, 𝜏𝛼0𝐷=101. .............................................................................................................................................................. 28 Figure 4.1 Proposed stages in ideal biofilm growth. Stage I (upper left) represent the initial situation, where a few isolated colonies (red points) tend to grow laterally (as indicated by black arrows at points L). Lateral growth continues during stage II (upper right), but spherical growth (arrows at points S) dominates in some portions. Stage III (lower left) emerges as the biofilm covers all the grains, and growth occurs primarily by thickening (red arrows at points S), although spherical growth may still occur in some places. Eventually, biofilm tends to clog the system (stage IV, lower right), so that flow lines tend to crowd the few open paths. This causes a dramatic reduction in permeability and increase in diffusion times into immobile zones, which may starve (purple zones). Changes in flow path geometries are moderate during stages I through III, although permeability may be reduced as a part of the pores is occupied by the biofilm, thus reducing its size. ................................................................. 37 Figure 4.2 Relationship between residence time and the fraction of porous medium occupied by biofilm. .............................................................................................................................................................. 37 Figure 4.3 The evolution and distribution of the volume fraction of biofilm in porous media. ........... 39 Figure 4.4 Biofilm growth curve plotted as a function of time at x = 9 cm, (left) plotted on linear scale and (right) semi-log scale. ..................................................................................................................... 40 Figure 4.5 Breakthrough curves of solute at x = 9cm for continuous injection at t =29 days. ............. 40 Figure 4.6 The evolution of residence times and probabilities in different portions of immobile zones during biofilm growth at t = 29 days..................................................................................................... 40 Figure 4.7 The evolution and distribution of the volume fraction of biofilm in porous media. ........... 41 Figure 4.8 Breakthrough curves of solute at x = 5cm for continuous injection at t=13, 62 and 113 days. .............................................................................................................................................................. 42 Figure 4.9 The evolution of residence times and probabilities in different portion of immobile zones during biofilm growth at x = 5cm and t = 13, 62 and 113 days. ........................................................... 42 Figure 4.10 The distribution of reaction rates in space and immobile zones at t = 13, 62 and 113 days from top to bottom. ............................................................................................................................... 43 Figure C.1 Breakthrough curves of species A, B and C for 𝑞𝑚=1.0 [LT-1], thus 𝑡𝑐=0.1 [T1], 𝜏𝑟1𝐷=𝜏𝑟2𝐷=102, 103 and 104 from top to bottom. ................................................................ 51 Figure C.2. Breakthrough curves of species A, B and C for 𝑞𝑚=0.1 [LT-1], thus 𝑡𝑐=1.0 [T1], 𝜏𝑟1𝐷=𝜏𝑟2𝐷=101, 102 and 103 from top to bottom. ................................................................ 52 xvii Figure C.3. Breakthrough curves of species A, B and C for 𝑞𝑚=0.01 [LT-1], thus 𝑡𝑐=10.0 [T1], 𝜏𝑟1𝐷=𝜏𝑟2𝐷=100, 101 and 102 from top to bottom. ................................................................ 53 xix List of tables Table 4-1 Model setup for transport and biofilm growth ...................................................................... 39 Table 4-2 Model setup for transport and biofilm growth ...................................................................... 41 1 1 Introduction This thesis is motivated by projects dealing with water renationalization by Soil Aquifer treatment (SAT) and degradation processes in the subsurface. The hydrogeology research group has worked recently on a methodology to enhance water quality improvement processes during soil passage (Valhondo et al., 2020). The approach consists of installing a reactive barrier at the bottom of infiltration basins in SAT system to favor the adsorption of organic pollutants and their degradation. To this end, the reactive barrier consists of a mixture of sand (to ensure high permeability and to guarantee structural strength), organic matter (woodchips or vegetable compost, to provide sorption sites for neutral compounds and as a source of organic carbon to promote reducing conditions), clay (to provide sorption sites for cationic compounds). Ideally, having a broad range of sorption sites and redox conditions should ensure the removal of most organic contaminants. Many of them hardly degrade under aerobic conditions, but require reducing conditions (Christensen et al., 2001). Another type of application is the removal of nitrates from mine tailing's seepage. An efficient passive technology to reduce nitrate concentration consists of circulating the nitrate loaded water through a tank containing woodchips, which ensures denitrification conditions (Nordström et al., 2021; Nordström & Herbert, 2019; Nordström & Herbert, 2018; Nordström & Herbert, 2017). Both approaches have been successfully tested in the field (Nordström et al., 2021; Valhondo et al., 2020; Valhondo, Martínez-Landa, et al., 2020; Nordström & Herbert, 2018; Valhondo et al., 2018). However, quantitative interpretation of results is difficult because degradation reactions occur in the biofilm, which evolves overtime. Microbial communities grow, which causes degradation rates to take place more efficiently, but also risks clogging. Quantitative interpretation of results is required for proper understanding and for upscaling observation from pilot to large scale and from relatively short (2 years) to long term operation. Quantitative interpretation is made by means of numerical models, which is challenging for several reasons. First, degradation reactions take place within the biofilm, which hosts the vast majority of microorganisms (Flemming et al., 2007; Morales et al., 2010; Flemming & Wingender, 2010; Satpathy et al., 2016; Flemming & Wuertz, 2019; Wu et al., 2019), but where water does not flow (Costerton et al., 1995; Flemming et al., 2007; Flemming & Wingender, 2010). Obviously, mobile-immobile representations of porous media are required. Many such representations are available (Le Borgne et al., 2008). But such representations are primarily motivated to reproduce the impact of heterogeneity, which leads to a broad range of residence times because of the variability of water velocity in different flow paths, but not necessarily because of the time it takes to reach actual immobile zones. Second, as a result of the above, the experience about simulating reactive transport in this kind of model representations is limited (Donado et al., 2009; Willmann et al., 2010). Worse, the computational demand increases dramatically because the number of unknowns is multiplied by the number of immobile zones used for representing the broad range of residence times. This may lead to computationally unfeasible problems when the chemical system is complex. And, third, these immobile regions evolve in time as a result of biofilm growth. Many studies exists on Its effect on hydrological parameters, such as permeability, porosity and retention curve (e.g., Taylor et al., 1990; Morales et al., 2010; Carles Brangarí et al., 2017; Lopez-Peña et al., 2019). Its effects on the parameters controlling mass exchange between mobile and immobile zones have been studied and modelled much less. CHAPTER 1 2 To address these challenges, I have opted for a multi-rate-mass-transfer (MRMT) approach (Haggerty & Gorelick, 1995; ), because it localizes chemical data, thus making it possible to use conventional chemical calculation tools, which I have also developed as part of the thesis. In this context, I have made several contributions. Each of them is described in one chapter, which will hopefully be published. This implies that every chapter can be read independently, but also that a portion of the mathematical description is repeated in every chapter. Chapter 2 contains a new method to solve MRMT problems. The method benefits from the fact that immobile zones are connected to independent mobile nodes. This leads to a highly sparse system matrix, which facilitates solution. In fact, the final system only contains the mobile zone concentrations because immobile zones concentrations are expressed as a function of mobile zone concentrations. The resulting algorithm is very fast. Chapter 3 analyses the interplay of immobile zones residence times, transport time and reaction times with two goals: first, to minimize the number of immobile zones, and second, to identify the conditions for chemical localization (i.e., the occurrence of reactions that would not take place if all the water was mobile). Finally, Chapter 4 contains a proposal for biofilm growth. To date, biofilm growth models emphasize clogging (i.e., the reduction of permeability as the pore space becomes filled by biofilm. Here, I emphasize that the residence time distribution changes in response to biofilm growth. This change is relevant for proper representation of the evolution of redox states that occur both during artificial recharge across a reactive barrier and during flow through a woodchips porous medium, which motivated this thesis. 3 2 A general and efficient numerical solution of reactive transport with multirate mass transfer  The presence of low permeability regions within porous media impacts solute transport and the distribution of species concentrations. Therefore, (bio)chemical reactions are equally affected. Multirate Mass Transfer (MRMT) models can be used to represent this anomalous transport process. MRMT conceptualizes the medium as a set of multiple continua: one mobile zone and multiple immobile zones. It simulates species transport in mobile and immobile zones simultaneously, which are related by firstorder mass exchange. Numerical modeling of reactive transport in this kind of multicontinua media is complex and demanding because of the high dimensionality of the problem. In this paper, we establish the governing equations of reactive transport in multicontinuum media incorporating chemical kinetics into the governing equations. We propose a general numerical solution of reactive transport with MRMT by applying direct substitution approach (DSA) based on Newton-Raphson method. The efficiency of the proposed algorithm benefits of the block structure of the system, which allows us to eliminate immobile zones equations and leads to significant savings in CPU time. We test the validity of the developed solution by comparison with other numerical and analytical solutions.  This chapter is based on the paper Wang et al., 2021. A general and efficient numerical solution of reactive transport with multirate mass transfer, submitted to Computers & Geosciences, under 2-nd review. CHAPTER 2 10 𝒓𝑚(𝑥,𝑡)=−𝜕𝒄𝑚,2(𝑥,𝑡) 𝜕𝒖𝑚(𝑥,𝑡)[𝜙𝑖𝑚 𝜙𝑚∫ 𝑓(𝛼) ∞ 0𝜕𝒖𝑖𝑚(𝑥,𝛼,𝑡) 𝜕𝑡 𝑑𝛼] +𝜙𝑖𝑚 𝜙𝑚∫ 𝑓(𝛼) ∞ 0𝜕𝒄𝑖𝑚,2(𝑥,𝛼,𝑡) 𝜕𝒖𝑖𝑚(𝑥,𝛼,𝑡)𝜕𝒖𝑖𝑚(𝑥,𝛼,𝑡) 𝜕𝑡 𝑑𝛼 −𝜕2𝒄𝑚,2(𝑥,𝑡) 𝜕𝒖𝑚 2(𝑥,𝑡)∇𝑇𝒖𝑚(𝑥,𝑡)𝑫∇𝒖𝑚(𝑥,𝑡)−𝜙𝑖𝑚 𝜙𝑚𝒓𝑖𝑚(𝑥,𝑡) (2.17) Figure 2.3 Distribution of (left) components 𝑢𝑚 and 𝑢𝑖𝑚, (right) reaction rates 𝑟𝑚 and 𝑟𝑖𝑚 versus time at distance 𝑥=100.0m. The blue color represents state variables in the mobile zone, and the red color represents state variables in immobile zones. The solid lines are the analytical solutions and the circle dot lines are the proposed numerical solutions. 2.3.1.3 Kinetic reaction To test the performance of our proposed solution in chemical kinetics, we choose a simple first-order kinetics. In this case, the analytical solution can be found in the Laplace domain (the solution is deduced in Appendix B). Then we simulate the distribution of components and kinetic rates both in mobile and immobile zones over time with the component degradation rate 𝜅=0.01. As we can see the proposed numerical solutions agree with the analytical solution as displayed in Figure 2.4. Comparing the evolution of components in chemical equilibrium (Figure 2.4 left), the presence of chemical reactions decreases the concentrations of the components both in mobile and immobile zones. Figure 2.4 Distribution of (left) components 𝑢𝑚 and 𝑢𝑖𝑚, (right) reaction rates 𝑟𝑚 and 𝑟𝑖𝑚 versus time at distance 𝑥=100.0m. The blue color represents state variables in the mobile zone, and the red color 11 represents state variables in immobile zones. The solid lines are the analytical solutions and the circle dot lines are the proposed numerical solutions. 2.3.2 Efficiency Analysis The efficiency of the algorithm depends on the problem size 𝑛, the total number of unknowns of system equation (2.13) that equals to the number of mesh nodes times the number of components times the number of immobile zones plus one. In our problem, the matrix operation is the most time consuming. Both the computational cost of matrix multiplication and matrix inversion are 𝑂(𝑛3). Instead of solving system equation (2.13) globally, we solve it in blocks representing mobile and immobile zones separately, using equations (2.14) and (2.15). Since these two approaches are mathematically identical, the two will converge to the same solution within the same iterations. The advantage of the proposed block solution is that it reduces the size of the system to be the number of mesh nodes times the number of components. We simulate an irreversible bimolecular reaction 𝐴+𝐵→𝐶 for the number of immobile zones equal to 3, 10, 30, 50, 70 and 100, with a second-order kinetics 𝑟𝑘=𝜅𝑐𝐴𝑐𝐵. In Figure 2.5, we compare the CPU time of the proposed method (i.e., block solver) and full DSA on a log-log scale. As we can see, the CPU time of the block solver increases linearly with the number of immobile zones (𝑁), while that of the full DSA increases much faster, approximating to 𝑁3. Clearly, the block solver runs faster than full DSA for a higher number of immobile zones. For a small number of immobile zones, the full DSA runs faster because of the costs of building the more complicated structure of the block solver. Figure 2.5 CPU time comparison between proposed method and full DSA for different number of immobile zones. The chemical systems are identical for all cases that is an irreversible bimolecular reaction 𝐴+𝐵→𝐶, with a second-order kinetics 𝑟𝑘=𝜅𝑐𝐴𝑐𝐵. 2.4 Conclusions The proposed method is effective and efficient for reactive transport modeling capable of accounting for numerous immobile zones. We formulate the general numerical solution of reactive transport with MRMT based on the Newton-Raphson method, which enables us to simulate complex chemical kinetics. For chemical systems, whatever it is in equilibrium or in kinetics, the convergent solution is efficiently solved within several iterations. CHAPTER 2 12 The proposed numerical solution is verified in case of passive solute transport, chemical equilibrium and chemical kinetics. In all cases the simulations agree very well with available analytical solutions and other numerical solutions. The full consistency between the proposed numerical solutions and available analytical solutions indicate that our proposed method is capable of reproducing the anomalous transport of reactive transport with MRMT. The computational efficiency of the proposed algorithm is improved by solving system equations in block instead of full DSA, which eliminates the impact of the number of immobile zones on the computational complexity and decreases the size of the linear system to be the number of mesh nodes times the number of components. This may become particularly important for more complicated chemical models in combination with MRMT. In that case the full DSA can become prohibitively expensive in CPU and the block solver may be the only feasible method. The advantage of the method increases with the number of immobile zones. Babey et al. (2015) concluded that five immobile zones (compared to one hundred in our models) are sufficient to address anomalous transport of non-reactive solutes. However, this may be different for reactive transport because of the broad range of residence times (Haggerty et al., 2000) and chemical localization (SolerSagarra et al., 2016). The number of immobile zones needed for proper reproduction of reactive processes requires further research. 2.5 Computer Code Availability Name of code: Reactive Transport with Multirate Mass Transfer Developer: Jingjing Wang Contact detail: Department of Civil and Environmental Engineering, Universitat Politécnica de Catalunya (UPC), Jordi Girona 1-3, 08034 Barcelona, Spain Email: [email protected] Year first available: October, 2020 Hardware required: No specific hardware is required. However, the code has been developed on a computer with processor Intel® Core™ i5-6500 CPU @ 3.20GHz, 16GB Installed memory (RAM), 64-bit Operating System, x64-based processor. Software development platform: Microsoft Visual Studio 2015 Program language: object-oriented programming FORTRAN 2003 Compiler: Intel Parallel Studio XE 2017 Cluster Edition for Windows* Library: Intel® Math Kernel Library. To solve the system equation (2.14), we call routine dgbtrf() to compute the LU factorization of the left hand side matrix of system equation (2.14), then we call routine dgbtrs() to solve the linear system with the LU-factored square coefficient matrix returned by routine dgbtrf(). Program size: 4.39 MB Details on how to access the open-source code: the source code can be freely download from GitHub on the public repository https://github.com/Jingjingwangxiang/RT_MRMT_DSA. Our codes are developed in object-oriented instead of procedural-oriented which lacks flexibility and extensibility (Meysman et al., 2003a ; Meysman et al., 2003b). The object-oriented programming allows the code reusability and facilitates the implementation of reactive transport modeling. To simulate the reactions, two main modules are developed. The biochemical system module simulates the localized chemical reactions occurring in a (bio)chemical system, it contains procedures that are capable of constructing the stoichiometric matrix and component matrix. The local biochemistry module captures 13 (bio)chemical state variables at each mesh node, such as components, concentrations, equilibrium reaction rates, kinetics, as well as the derivatives of kinetics with respect to concentrations, etc. It encapsulates procedures of chemical calculations that are capable of computing the state variables at each mesh node by using the biochemistry system module. These two main developed modules are coupled with the transport equations in the reactive transport module, which enables the modeling of reactive transport. 15 3 On the localization of chemical reactions in multicontinuum media  The objective of this work is to study the localization of reactions in different portions of the porous medium. For this we use the non-local MRMT model, which views the porous media as being composed of one mobile zone and many immobile zones. The localized physical and chemical heterogeneity are reflected by a distribution of residence times and reaction times in the different immobile zones, which can be incorporated in the governing equations through memory functions and reaction rate kernels, respectively. We show that the dimensionless form of governing equations is characterized by reaction times, transport times and distribution of residence times. To analyze the interplay between them, we simulated three cases: one without reactions, one with a single reaction and one with sequential reactions. The results indicate that reactions driven by species that are not present in the inflowing water but are the result of previous reactions will take place in immobile zones, whose residence time is comparable to reaction times. Furthermore, mobile zones with residence times much smaller than those for transport can be lumped together (assuming that very fast reactions are assumed in equilibrium), which greatly reduces computations.  This chapter is based on the manuscript Wang et al., 2021. On the localization of chemical reactions in multicontinuum media, submitted to Water, under review. CHAPTER 3 16 3.1 Introduction Characterization of reactive transport in heterogeneous media is important for environmental science and engineering problems, such as managed aquifer recharge, seawater intrusion, CO2 geological storage, and in situ bioremediation. Physical heterogeneity may induce chemical heterogeneity or localization of chemical reactions by which different reactions occur in different portions of the porous medium (Soler-Sagarra et al., 2016; Babaei & Islam, 2018). Dentz et al. (2011) concluded that transport under physical and chemical heterogeneity cannot be upscaled separately: upscaling reactions depends on physical heterogeneity and upscaling transport is affected by chemical heterogeneity. In addition, we acknowledge that the behavior of macroscale transport in the mobile zone of the porous medium is determined by the microscale transport process in immobile zones (Gouze et al., 2008). Yet, it is not so clear how the localized physical and chemical heterogeneity affects the distribution and evolution of reactive species in different portion of immobile zones and how these microscale processes in immobile zones govern the reactive transport in mobile zone? It is well known that transport in heterogeneous porous media displays anomalous (non-Fickian) behavior (Kitanidis, 1988) both at field (Adams & Gelhar, 1992) and laboratory scales (Valocchi, 1985; Levy & Berkowitz, 2003; Berkowitz & Scher, 2009). Hence, the advection-dispersion-reaction equation (ADRE) is no longer appropriate. Many non-local methods have been developed to represent this anomalous transport, including continuous time random walks (CTRW) (Berkowitz & Scher, 1998; Dentz et al., 2004; Berkowitz et al., 2006; Berkowitz et al., 2006; Dentz et al., 2015), fractional advection-dispersion equations (fADE) (Benson et al., 2000; Schumer et al., 2003; Marseguerra & Zoia, 2008), multirate mass transfer (MRMT) (Haggerty & Gorelick, 1995; Wang et al., 2005; Salamon et al., 2006; Benson & Meerschaert, 2009; Dentz et al., 2011; Fernàndez-Garcia & Sanchez-Vila, 2015; De Dreuzy & Carrera, 2016), memory functions (Carrera et al., 1998; Haggerty et al., 2000; Willmann et al., 2008; Gouze et al., 2008) and so forth. Although these methods use different approaches, essentially they are equivalent (Dentz & Berkowitz, 2003; Silva et al., 2009; Neuman & Tartakovsky, 2009). Many studies focus on the behavior of breakthrough curves (BTCs) and conservative transport in heterogeneous porous media (Haggerty et al., 2000; Schumer et al., 2003; Willmann et al., 2008; Berkowitz & Scher, 2009; Dentz et al., 2015). Only few deals with reactive transport (Willmann et al., 2010; Dentz et al., 2011) probably due the complexity of reactive transport in heterogeneous porous media. For non-linear kinetic reactions, no analytical solutions exist in the Laplace domain which further limits the study of the problem. Among the non-local methods, the MRMT formulation allows localized concentrations which simplifies the simulations of reactive transport in physical and chemical heterogeneous media (Dentz et al., 2011; Soler-Sagarra et al., 2016; Babaei & Islam, 2018). Moreover, the MRMT formulation for conservative transport can be extended to reactive transport, including nonlinear reactions (Willmann et al., 2010). The objective of this work is to study the effects of localized physical and chemical heterogeneity on reactive transport based on the MRMT formulation. To do so, first, we establish the governing equations of reactive transport in multicontinuum media and deduce a dimensionless form of these equations. Then we present three models one without reactions, one with a single reaction, and one with two sequential reactions in section 3.3, whose results are analyzed in section 3.4. The last section is dedicated to the conclusions. 3.2 Methodology We model reactive transport in heterogeneous porous medium by using the MRMT approach (Haggerty & Gorelick, 1995). Every point (representative elementary volume) is viewed as consisting of a mobile zone and a distribution of immobile zones characterized by their residence time. Each of these 17 exchanges solute mass with the mobile portion proportionally to the mass transfer rate, 𝛼[T−1]. Therefore, it is natural to characterize this distribution by a probability density function (pdf) of mass transfer rates, 𝑓(𝛼) (see Figure 3.1 left). To facilitate comparison to other non-local methods, we also characterize immobile zones by the distribution, 𝑃(𝜏𝛼), of residence times, 𝜏𝛼≡1 𝛼 ⁄ [T] (Haggerty et al., 2004) (see Figure 3.1 right). 𝑃(𝜏𝛼) is given by (Dentz & Berkowitz, 2003; Dentz et al., 2011; Dentz et al., 2015) 𝑃(𝜏𝛼)=𝜏𝛼 −2𝑓(𝛼) (3.1) The residence time probability 𝑃(𝜏𝛼)𝑑𝜏𝛼 is the frequency of the immobile zone, characterized by the residence time in the interval [𝜏𝛼,𝜏𝛼+𝑑𝜏𝛼]. Obviously, it satisfies the condition ∫𝑃(𝜏𝛼) ∞ 0𝑑𝜏𝛼=1. Figure 3.1 Pdf of exchange rates, 𝛼 (Gamma distibution), and residence times, 𝜏𝛼, plotted in log-log scale for several values of 𝜏𝛼0 and 𝛽=1 2 ⁄. Note that the maximum of latter is around 𝜏𝛼0 (the mode is 𝜏𝛼0/(𝛽+1) and the expected value 𝜏𝛼0 (1−𝛽) ⁄). Its large 𝜏𝛼 slope is 𝛽+1. 3.2.1 Governing equations The transport of any reactive species in a medium that consists of a mobile continuum and multiple immobile continua can be expressed as (Donado et al., 2009; Willmann et al., 2010; Dentz et al., 2011) 𝜙𝑚𝜕𝑐𝑚 𝜕𝑡 =𝐿𝑡[𝑐𝑚]−𝜙𝑖𝑚∫1 𝜏𝛼𝑃(𝜏𝛼)[𝑐𝑚(𝑥,𝑡)−𝑐𝑖𝑚(𝑥,𝜏𝛼,𝑡)] ∞ 0𝑑𝜏𝛼+𝜙𝑚𝑟𝑚 (3.2) where 𝐿𝑡[𝑐𝑚]=−[𝑞𝑚𝛻𝑐𝑚]+𝛻∙[𝜙𝑚𝐷𝑚𝛻𝑐𝑚] is the transport operator that describes advection and dispersion, 𝑞𝑚 [LT-1] is Darcy flux, 𝐷𝑚 [L2T-1] is the hydrodynamic dispersion tensor; 𝜙𝑚 [-] and 𝜙𝑖𝑚 [-] are the porosities of mobile and immobile zones, which denote the volume fraction of the mobile pore zone and immobile pore zones over the bulk volume, respectively; 𝑐𝑚 [ML-3] and 𝑐𝑖𝑚 [ML-3] are concentrations in the mobile and immobile zone, which are expressed as mass per unit volume of mobile zone and mass per unit volume of immobile zone; 𝑟𝑚 [ML-3T-1] is the sink-source term due to chemical reactions in mobile zone and corresponds to the mass removed by reactions in mobile zone per unit volume of mobile water per unit time. Mass balance in each immobile zone is given by 𝜕𝑐𝑖𝑚(𝑥,𝜏𝛼,𝑡) 𝜕𝑡 =1 𝜏𝛼[𝑐𝑚(𝑥,𝑡)−𝑐𝑖𝑚(𝑥,𝜏𝛼,𝑡)]+𝑟𝑖𝑚(𝑥,𝜏𝛼,𝑡) (3.3) where 𝑟𝑖𝑚 [ML-3T-1] is the reactions sink-source term in the immobile zone with residence time 𝜏𝛼, which corresponds to the mass removed by reactions in the immobile zone per unit volume of immobile water per unit time. The mass exchange between mobile and immobile zone is modeled by a first-order mass transfer mechanism, represented by the continuous variable 𝜏𝛼 and characterized by the distribution density function 𝑃(𝜏𝛼) for the immobile zones. Thus, the total mass exchange is the weighted sum over all CHAPTER 3 18 immobile zones as expressed in the integral term of equation (3.2). Integrating equation (3.3) with weight 𝑃(𝜏𝛼) in terms of 𝜏𝛼, multiplying by 𝜙𝑖𝑚, and substituting it into equation (3.2), yields the total solute mass balance, 𝜙𝑚𝜕𝑐𝑚(𝑥,𝑡) 𝜕𝑡 +𝜙𝑖𝑚∫ 𝑃(𝜏𝛼)𝜕𝑐𝑖𝑚(𝑥,𝜏𝛼,𝑡) 𝜕𝑡 ∞ 0𝑑𝜏𝛼=𝐿𝑡[𝑐𝑚(𝑥,𝑡)]+𝑟(𝑥,𝑡) (3.4) in which, 𝑟 is the total reaction rate (now per unit volume of bulk porous medium) that integrates reactions in both mobile and immobile zones, 𝑟(𝑥,𝑡)=𝜙𝑚𝑟𝑚(𝑥,𝑡)+𝜙𝑖𝑚∫ 𝑃(𝜏𝛼)𝑟𝑖𝑚(𝑥,𝜏𝛼,𝑡) ∞ 0𝑑𝜏𝛼 (3.5) If the reaction follows first-order kinetics, the reaction rate in mobile zone will be 𝑟𝑚(𝑥,𝑡)=−𝑘𝑚𝑐𝑚(𝑥,𝑡) (3.6) where 𝑘𝑚 is the reaction rate constant. Similarly, the local reaction rate in the 𝜏𝛼 immobile zone becomes 𝑟𝑖𝑚(𝑥,𝜏𝛼,𝑡)=−𝑘𝑖𝑚(𝜏𝛼)𝑐𝑖𝑚(𝑥,𝜏𝛼,𝑡) (3.7) where 𝑘𝑖𝑚(𝜏𝛼) [T-1] is the local reaction rate constant. The reaction time in the immobile zone of 𝜏𝛼 is 𝜏𝑟≡1 𝑘𝑖𝑚(𝜏𝛼)⁄ , which is a variable dependent on the immobile zone. By inserting equation (3.7) into (3.3), then solving equation (3.3), we obtain the concentration in the immobile zone of 𝜏𝛼 as a function of mobile concentration history 𝑐𝑖𝑚(𝑥,𝜏𝛼,𝑡)=𝑐𝑖𝑚(𝑥,𝜏𝛼,𝑡=0)𝑒−[𝛼+𝑘𝑖𝑚(𝜏𝛼)]𝑡+∫𝜑(𝑡−𝑡′)𝑐𝑚(𝑥,𝑡′) 𝑡 0𝑑𝑡′ (3.8) in which, the definition of 𝜑(𝑡) is 𝜑(𝑡)=1 𝜏𝛼𝑒−[𝛼+𝑘𝑖𝑚(𝜏𝛼)]𝑡 (3.9) Assuming that the initial concentration in the immobile zones is zero (i.e., 𝑐𝑖𝑚(𝑥,𝜏𝛼,𝑡=0)=0), then by substituting equation (3.8) into equation (3.2) and rearranging terms, we obtain the total solute mass balance, as a sole function of 𝑐𝑚 𝜙𝑚𝜕𝑐𝑚(𝑥,𝑡) 𝜕𝑡 +𝜙𝑖𝑚 𝜕 𝜕𝑡∫𝑔(𝑡−𝑡′)𝑐𝑚(𝑥,𝑡′)𝑑𝑡′ 𝑡 0 =𝐿𝑡[𝑐𝑚(𝑥,𝑡)]−∫𝜅(𝑡−𝑡′)𝑐𝑚(𝑥,𝑡′)𝑑𝑡′ 𝑡 0 (3.10) where the memory function (Carrera et al., 1998; Haggerty et al., 2000) and reaction rate kernel (Dentz et al., 2011) are given by 𝑔(𝑡)=∫ 𝑃(𝜏𝛼)𝜑(𝑡)𝑑𝜏𝛼 ∞ 0 (3.11) and 𝜅(𝑡)=𝜙𝑚(𝑥)𝑘𝑚𝛿(𝑡)+𝜙𝑖𝑚(𝑥)∫1 𝜏𝑟𝑃(𝜏𝛼)𝜑(𝑡)𝑑𝜏𝛼 ∞ 0 (3.12) respectively. In which, 𝛿(𝑡) is the Dirac delta. The memory function can be viewed as the rate of change of concentration in the immobile zone that is caused by a unit change of concentration in the mobile zone at initial time 𝑡=0. In the presence of reactions, the memory function not only incorporates the 19 distribution of local mass transfer rates but also the distribution of local reaction rates. Similarly, the reaction rate kernel represents the local reactions in immobile zones, it integrates both the reactions in mobile and immobile zones, and accounts for the mass transfer between mobile and immobile zones. Although the memory function is mainly controlled by the distribution of residence times in immobile zones, it is also affected by the local reaction times in immobile zones. Likewise, the reaction rate kernel is mainly controlled by the chemical reactions, but it is also affected by the distribution of residence times in immobile zones (Dentz et al., 2011). Clearly, the governing equations of reactive transport under physical and chemical heterogeneity (3.10), (3.11) and (3.12) are controlled by the distribution of residence times and reaction times in immobile zones simultaneously. The physical and chemical heterogeneities interact with each other, and together they govern reactive transport (Dentz et al., 2011). 3.2.2 Dimensionless formulations To characterize the governing equation (3.10), (3.11) and (3.12), we define the characteristic length as follows 𝐿𝑐=𝜙𝑚𝐷𝑚 𝑞𝑚 (3.13) and the characteristic transport time as follows 𝑡𝑐=𝐿𝑐 𝑣𝑚=𝜙𝑚 2𝐷𝑚 𝑞𝑚 2 (3.14) Note, that in 1D, the characteristic length equals the longitudinal dispersivity 𝛼𝐿 due to the definition of dispersion 𝐷𝑚=𝛼𝐿𝑣𝑚, where 𝑣𝑚=𝑞𝑚𝜙𝑚 ⁄ is the mean fluid velocity in porous media. Introducing these definitions of characteristic length and characteristic transport time into equation (3.10), we obtain the dimensionless form of the governing equations, that is 𝜕𝑐𝑚𝐷 𝜕𝑡𝐷+𝜂 𝜕 𝜕𝑡𝐷∫ 𝑔(𝑡𝐷−𝑡𝐷 ′)𝑐𝑚𝐷(𝑡𝐷 ′)𝑑𝑡𝐷 ′ 𝑡𝐷 0=−𝛻𝐷𝑐𝑚𝐷+𝛻𝐷∙𝛻𝐷𝑐𝑚𝐷−∫ 𝜅(𝑡𝐷−𝑡𝐷 ′)𝑐𝑚𝐷(𝑡𝐷 ′)𝑑𝑡𝐷 ′ 𝑡𝐷 0 (3.15) with the following dimensional variables 𝑐𝑚𝐷=𝑐𝑚 𝑐𝑐, 𝑡𝐷=𝑡 𝑡𝑐, 𝑥𝐷=𝑥 𝐿𝑐, 𝜂=𝜙𝑖𝑚 𝜙𝑚, 𝜏𝛼𝐷=𝜏𝛼 𝑡𝑐, 𝜏𝑟𝐷=𝜏𝑟 𝑡𝑐 (3.16) In the dimensionless formulations of governing equation (3.15), the physical and chemical heterogeneity are represented by the distribution of residence times and reaction times, respectively, and reflected simultaneously in the memory function and reaction rate kernel. Obviously, in the situation where the characteristic length 𝐿𝑐 and the dimensionless ratio of porosity of immobile zone to porosity of mobile zone 𝜂 are invariables, the governing equation (3.15) in dimensionless form are totally governed by three characteristic times, that is the characteristic transport time 𝑡𝑐, the dimensionless residence time 𝜏𝛼𝐷 and the dimensionless reaction time 𝜏𝑟𝐷 in immobile zones. The solution of governing equation (3.15) in Laplace domain is expressed as ℒ{𝑐𝑚𝐷}(𝑠)=𝑒𝑥𝑝[𝑥𝐷 2(1−√1+4[𝑠[1+𝜂ℒ{𝑔}(𝑠)]+ℒ{𝜅}(𝑠)])] (3.17) where ℒ{∙}(𝑠) represents the Laplace transform of a function. The Laplace transform of the memory function and reaction rate kernel are defined as CHAPTER 3 26 distribution with 𝛽=1 2 ⁄, 𝜏𝛼0𝐷=101. The reaction times are uniform in immobile zones with 𝜏𝑟𝐷= ∞, 104, 103 and 102 from top to bottom. Figure 3.8 Cumulative concentration profiles versus dimensionless distance and dimensionless residence times for species A (left) and B (right) at 𝑡𝐷=300 (above), and 1×104 (below) from top to bottom after continuous injection of A in a medium with immobile zones distributed with 𝛽=1 2 ⁄, 𝜏𝛼0𝐷=101. Reaction A→B occurs with a characteristic time 𝜏𝑟𝐷=102. 3.4.3 Sequential reaction system Results for the sequential reactions case (A→B and B→C) are summarized in Figure 3.9 and Figure 3.10. Figure 3.9 makes it clear that BTCs for mobile-immobile regions with short residence times are identical to BTCs where all the porosity is mobile. That is, immobile regions with short residence time (compared to transport time) can be lumped in the mobile region, if reactivity is the same in both types of regions. Otherwise (i.e., when reactions occur primarily in the immobile region), all fast reactions can be lumped into a simple immobile region (a complete set of BTCs for several 𝜏𝑟𝐷 and 𝜏𝛼0𝐷 is presented in Appendix C). The most significant observation derives from Figure 3.10, which displays immobile regions concentrations as a function of space and residence time. This figure makes it clear that C concentrations are only relevant in immobile regions with residence times comparable to characteristic reaction times, regardless of the transport time. Therefore, the three sets of times (transport time, reaction time, and residence times) are relevant when deciding the appropriate discretization of residence times. 27 Figure 3.9 Breakthrough curves of species A, B and C in different cases for 𝜏𝑟1𝐷=𝜏𝑟2𝐷=103. The black lines represent results only with mobile zone. Case 1: 𝜏𝛼𝐷=[100], case 2: 𝜏𝛼𝐷=[100,101], case 3: 𝜏𝛼𝐷=[100,101,102], case 4: 𝜏𝛼𝐷=[100,101,102,103], case 5: 𝜏𝛼𝐷= [100,101,102,103,104], case 6: 𝜏𝛼𝐷=[100,101,102,103,104,105,106], case 7: 𝜏𝛼𝐷=[102,103] and case 8: 𝜏𝛼𝐷=[101,102,103]. CHAPTER 3 28 Figure 3.10 Cumulative concentration profiles versus dimensionless distance and dimensionless residence times for species (left) A, (middle) B and (right) C for 𝜏𝑟1𝐷=𝜏𝑟2𝐷=103 at 𝑡𝐷=50, 200, 500, 2000 and 1×104 from top to bottom. Species transport in a medium characterized by the immobile zones in which the residence time follows gamma distribution with 𝛽=1 2 ⁄, 𝜏𝛼0𝐷=101. 3.5 Conclusions Several conclusions can be drawn from the modeling exercise presented in this work 1) Some reactions may occur in immobile regions that would not occur if the entire medium was mobile. This is not surprising when the mineral composition is different in each immobile zone, as was discussed by Soler-Segarra et al (2016). What we show here is that, independently of the reactivity of mobile and immobile regions, the interplay between residence times and reaction times may cause some reactions to take place solely in the immobile region. This will occur whenever residence time is comparable to reaction time, and both are much longer than the transport time. It goes without saying that immobile region reactions will be enhanced if reactivity in immobile regions is higher than in the mobile region. This is the case for biofilms, where microbial mediated degradation reactions concentrate. 2) Regardless of reaction rates, concentrations in the fast immobile regions (i.e., regions where 𝜏𝛼≪ 𝑡𝑡=𝐿𝜙/𝑞) will tend to equilibrate with mobile concentrations. Therefore, little is gained by representing them explicitly in a model. All fast exchange immobile regions can be lumped into a zone with residence time slightly smaller than transport time. 29 3) Similarly, very slow immobile zones (𝜏𝛼≫ 𝑡𝑡) can be ignored because little mass will diffuse into them. Note that this conclusion is more relevant for laboratory experiments or short tracer tests, where transport forms are moderate. Under natural conditions, transport time can be very long (many years), so that most immobile regions are indeed accessible. 31 4 A conceptual model for reactive transport with dynamic biofilm growth in multicontinuum media  Biofilm growth in porous media changes the hydrodynamic properties of the medium: porosity and permeability are reduced, and dispersivity increases. However, the first arrival of breakthrough curves (BTCs) is more reduced than derived from the reduction in porosity, and the BTC tail becomes heavier. These observations suggest the need of multicontinuum models (Multirate-Mass-Transfer, MRMT) that evolve dynamically with the biofilm. The MRMT model is capable of representing reactive transport in heterogeneous porous media which facilitates the simulation of localized reactions often observed within biofilms. In this work, we present a conceptual model of reactive transport with dynamic biofilm growth based on MRMT formulations. It incorporates the microbial growth according to the stoichiometry and kinetic rate laws of biological reactions. The physical, including not only porosity, but also the distribution of residence times in immobile zones, and chemical properties are updated after the reactive transport simulations at each time. This model has been tested on laboratory data.  This chapter is based on the manuscript Wang et al., 2021. A conceptual model for reactive transport with dynamic biofilm growth in multicontinuum media, in preparation for Water Resource Research. CHAPTER 4 32 4.1 Introduction Biofilm growth in porous media is important for many natural and engineered environmental applications such as wastewater treatment (Nicolella et al., 2000; Lewandowski & Boltz, 2011), aquifer recharge (Kim et al., 2010), enhanced oil recovery (Raiders et al., 1989; Van Hamme et al., 2003) and in situ bioremediation of soil and aquifer contamination (Cunningham et al., 2003; Sethi & Di Molfetta, 2019). Biofilm consists of living and reproducing microbial communities. It forms when microbes adhering to the surface of grains excrete extracellular polymeric substance (EPS). EPS is a fibrous geltype matrix composed of polysaccharides, lipids, proteins, and DNA materials (Flemming et al., 2007; Flemming & Wingender, 2010). Biofilm mainly consists of water. Microbial cells, which account for less than 10% of the dry mass, and the EPS matrix, which can account for over 90% (Satpathy et al., 2016). Microorganisms in biofilms catalyze chemical reactions (Flemming & Wingender, 2010). The metabolism of microorganisms leads to the modification of biofilm surface structure and the deposition of microbes (Esperanza Cortés et al., 2011). The complex biofilm structure and metabolism gives the biofilm a function similar to tissues in higher organisms (Costerton et al., 1995). Models usually take microbes as an independent phase to account for biofilm growth (Chen-Charpentier, 1999; Von Der Schulenburg et al., 2009; Shafahi & Vafai, 2009). These models assume that biochemical reactions take place in the bulk mobile pore water. In reality, however, these reactions are catalyzed by microbes residing within the immobile biofilm. Moreover, the biofilm structure is very heterogeneous in space and evolves with biomass growth (Picioreanu et al., 1998; Eberl et al., 2000; Picioreanu et al., 2004). Biofilm growth is known to change the hydrodynamic properties in porous media such as reduction of porosity and permeability (Taylor et al., 1990; Cunningham et al., 1991; Von Der Schulenburg et al., 2009; Cunningham et al., 2003; Thullner, 2010), retention curve (Carles Brangarí et al., 2017), as well as dispersivity (Taylor & Jaffé, 1990b; Kone et al., 2014). As a result, it impacts the reactive transport processes, which in turn affects biofilm growth. Experimental results demonstrate that biofilm growth induces heterogeneities that affect the transport in porous media (Seifert & Engesgaard, 2007; Kone et al., 2014). To properly simulate biofilm growth and its effects on reactive transport, a model is needed that accounts for these heterogeneities. However, biofilm growth modeling is complex, because it requires coupling it with reactive transport. The challenges include how to model reactive transport in heterogeneous media, how to model the biofilm growth, and how to couple these two processes reasonably. As the hydraulic conductivity of biofilm is very low (Deng et al., 2013), exchange of solutes between pore water and biofilm is controlled by diffusion. Models of biofilm growth on a flat surface in constantly stirred tanks take this into account (Picioreanu et al., 2004). In porous media this has to be combined with the heterogeneity of the pore and biofilm sizes. This makes the MRMT a promising tool to model the dynamical interaction between solute transport and biofilm growth. In fact, diffusion models of mass transfer between mobile and immobile zones can also be described by a set of firstorder mass transfer models (Haggerty & Gorelick, 1995). The non-local MRMT formulation has been widely used to represent transport in heterogeneous media (Haggerty & Gorelick, 1995; Wang et al., 2005; Salamon et al., 2006; Benson & Meerschaert, 2009; Dentz et al., 2011; Fernàndez-Garcia & Sanchez-Vila, 2015; De Dreuzy & Carrera, 2016). In addition, Willmann et al. (2010) obtained quite accurate agreement between 2D heterogeneity and 1D MRMT, which demonstrates that the non-local MRMT formulation of conservative transport can be extended to reactive transport. The MRMT model allows the localized modeling of reactive transport in physical and chemical heterogeneous porous media (Soler-Sagarra et al., 2016; Babaei & Islam, 2018). The physical and chemical heterogeneity are characterized by a distribution of mass transfer rates and reactions rates (or residence times and reaction times) in immobile zones. 33 The major challenge is that biofilm growth causes the heterogeneity of porous media to change. In fact, it is clear that velocity distributions in the medium evolve as the biofilm grows. The velocity distribution becomes broad, which suggests an increase in channeling and immobile zones (Maxence Carrel et al., 2018; M. Carrel et al., 2018). Hence, the evolution of biomass distribution should be represented by the change of mass transfer rates and reaction rates in biofilm. The objective of this work is to build a concept model of biofilm growth in the framework of MRMT. We focus on the study of biofilm growth induces effects on transport. 4.2 Conceptual models In our concept model, we assume a biofilm consists of water and biomass i.e., EPS & microbial cells (Seifert & Engesgaard, 2007). The biofilm is represented by the various immobile zones of the MRMT approach. Solutes diffuse into biofilm and provide nutrients for the metabolism of microbes residing in the biofilm. Reactions mediated by microbes only occur in the biofilm, not in mobile water. Here, we explain the reactive transport formulation of the abiotic chemical species, and the incorporation of biofilm. 4.2.1 Reactive transport with MRMT In MRMT model, reactive transport in mobile zone is expressed as 𝜙𝑚𝜕𝒄𝑚 𝜕𝑡 =𝐿𝑡[𝒄𝑚]−∑𝜙𝑖𝑚,𝑗 𝑁 𝑗=1 𝛼𝑗(𝒄𝑚−𝒄𝑖𝑚,𝑗)+𝜙𝑚𝒇𝑄𝑚 (4.1) in which, 𝜙𝑚 [-], porosity of mobile water and 𝜙𝑖𝑚 [-], porosity of immobile water , denote the volume fraction of the mobile and immobile pore water over the total volume, respectively; 𝒄𝑚 [ML-3] and 𝒄𝑖𝑚,𝑗 [ML-3], arrays of concentrations of all aqueous species in the mobile and 𝑗th immobile zone, denote the mass of a species per unit volume of mobile water and immobile water. 𝐿𝑡[𝒄𝑚]=−𝒒𝑚𝛻𝒄𝑚+𝛻∙ (𝜙𝑚𝑫𝑚𝛻𝒄𝑚) is the transport operator that describe the advection and dispersion processes, 𝒒𝑚 [LT-1] is Darcy flux, 𝑫𝑚 [L2T-1] is the hydrodynamic dispersion tensor that includes the molecular diffusion and the mechanical dispersion. The second term on the right-hand side of equation (4.1) describes the mass transfer between mobile water and immobile water. 𝛼𝑗 [T-1] is a first-order mass transfer rate between the mobile and 𝑗th immobile zone. 𝜙𝑖𝑚,𝑗=𝜙𝑖𝑚𝑝𝑗, where 𝑝𝑗 is the probability of the 𝑗th immobile zone. The sum of 𝑝𝑗 of all immobile zones should be equal to one, that is, ∑𝑝𝑗 𝑁 𝑗=1 =1.0. 𝒇𝑄𝑚 [ML-3T-1] is the sink-source term due to reactions in the mobile zone which is the mass removed 𝒇𝑄𝑚<0 or added 𝒇𝑄𝑚>0 by reactions per unit volume of mobile water per unit of time. In the immobile zone, the mass balance is described by 𝜕𝒄𝑖𝑚,𝑗 𝜕𝑡 =𝛼𝑗(𝒄𝑚−𝒄𝑖𝑚,𝑗)+𝒇𝑄𝑖𝑚,𝑗 𝑗=1,⋯,𝑁 (4.2) in which, 𝒇𝑄𝑖𝑚,𝑗 [ML-3T-1] is the sink-source term of the 𝑗th immobile zone and corresponds to the mass removed or added by biochemical reactions per unit volume of immobile water and per unit of time, and 𝑁 is the number of immobile zones. Equations (4.1) and (4.2) need to be complemented with the mass balance of non-mobile species (biomass, minerals, sorbed species). These additional mass balances are identical to those of aqueous species, but neglecting transport terms (advections, diffusion, and exchange) The basic equations for biochemical sink/sources can be expressed as a function of reaction rates (Steefel & MacQuarrie, 1996; Saaltink et al., 1998) CHAPTER 4 34 𝒓=𝑺𝑒𝑇𝒓𝑒+𝑺𝑘 𝑇𝒓𝑘 (4.3) where 𝑺𝑒 is the stoichiometric matrix for equilibrium reactions, and 𝑺𝑘 is the stoichiometric matrix for kinetic reactions, 𝒓𝑒 and 𝒓𝑘 represent the vectors of reaction rates for equilibrium and kinetic reactions, respectively. Note that these matrices contain the full description of the chemical system, as every row contains the stoichiometric coefficients of all species (columns) participating in the reaction. In the case of biofilm growth, these species may include biomass, whose level of detail depends on the specific model, as well as electron acceptors and donors. Solving the system of equations (4.1), (4.2) and (4.3), including the mass balance of non-mobile species, require complementing these equations with the Mass Action Law for equilibrium reactions and empirical expressions. Solution can be complex, but it is greatly simplified by introducing the component matrix 𝑼 to eliminate the equilibrium reactions and reduce the number of unknowns in the governing equations (4.1) and (4.2). It is the kernel of 𝑺𝑒𝑇, defined as 𝑼𝑺𝑒𝑇=𝟎 (Saaltink et al., 1998; Molins et al., 2004). Multiplying equation (4.1) and (4.2) by the component matrix 𝑼, we obtain the governing equations of the components, 𝒖=𝑼𝒄 that is 𝜙𝑚𝜕𝒖𝑚 𝜕𝑡 =𝐿𝑡[𝒖𝑚]−∑𝜙𝑖𝑚,𝑗 𝑁 𝑗=1 𝛼𝑗(𝒖𝑚−𝒖𝑖𝑚,𝑗)+𝜙𝑚𝑼𝑺𝑘 𝑇𝒓𝑘,𝑚 (4.4) 𝜕𝒖𝑖𝑚,𝑗 𝜕𝑡 =𝛼𝑗(𝒖𝑚−𝒖𝑖𝑚,𝑗)+𝑼𝑺𝑘 𝑇𝒓𝑘,𝑖𝑚,𝑗 𝑗=1,⋯,𝑁 (4.5) 4.2.2 Biofilm growth In this study we simplify the biochemistry by two kinetic reactions: the growth of microbes or biomass (b) through a redox reaction involving an electron donor (D) and acceptor (A) and the death of biomass. Other chemical species, such as inorganic carbon, are not considered. Then, the stoichiometric matrix for kinetic reactions becomes: 𝑺𝑘=[−𝑆𝐷−𝑆𝐴] (4.6) where SD and SA are stoichiometric coefficients. We do not consider equilibrium reactions, which means that matrix U equals the identity matrix. The mass balances of non-mobile species. The metabolism of microbes in biofilm induces the growth and decay of biofilm which is given by 𝜌𝑏𝜕𝜙𝑏,𝑖𝑚,𝑗 𝜕𝑡 =𝑟𝑏,𝑖𝑚,𝑗𝜙𝑏,𝑖𝑚,𝑗 (4.7) where 𝜌𝑏 [ML-3] is the molar density of biomass (biomass per unit volume of biofilm), 𝜙𝑏,𝑖𝑚,𝑗= 𝜙𝑖𝑚,𝑗 𝜃𝑏 ⁄ [-] is the volume of biofilm in the 𝑗th immobile zone per unit volume of porous medium, 𝜃𝑏 is the volume ratio of immobile water in the biofilm, and 𝑟𝑏,𝑖𝑚,𝑗 [ML-3T-1] is the growth rate of microbes corresponding to the mass growth of microbes per unit volume of biofilm per unit of time that can be expressed as (Cirpka et al., 1999; Rodríguez-Escales et al., 2016) 𝑟𝑏,𝑖𝑚,𝑗=𝑌𝑟𝑘,𝑖𝑚,𝑗−𝑑𝜌𝑏 (4.8) where 𝑌 [-] is the yield coefficient which denotes the production of microbes per unit mass of substrate (electron donor), 𝑑 [T-1] is the death rate, 𝑟𝑘,𝑖𝑚,𝑗 [ML-3T-1] is the degradation rate of substrate (electron acceptor) catalyzed by microbes, which can be represented by the Monod kinetics equation (RodríguezEscales et al., 2016) 35 𝑟𝑘,𝑖𝑚,𝑗=𝜇 𝑐𝐷,𝑖𝑚,𝑗 𝐾𝐷+𝑐𝐷,𝑖𝑚,𝑗 𝑐𝐴,𝑖𝑚,𝑗 𝐾𝐴+𝑐𝐴,𝑖𝑚,𝑗𝜌𝑏 (4.9) where 𝜇 [T-1] is the maximum growth rate, 𝑐𝐷,𝑖𝑚,𝑗 [ML-3] and 𝑐𝐴,𝑖𝑚,𝑗 [ML-3] are the concentrations of electron donor and electron acceptor in biofilm, 𝐾𝐷 [ML-3] and 𝐾𝐴 [ML-3] are the half saturation constants for species of electron donor and electron acceptor, respectively. The solution of equation (4.7) is 𝜙𝑏,𝑖𝑚,𝑗(𝑡)=𝜙𝑏,𝑖𝑚,𝑗(𝑡=0)𝑒𝑟𝑏,𝑖𝑚,𝑗 𝜌𝑏𝑡 (4.10) It is an exponential function with its exponent depending on the microbial growth rate and molar density that equals to 𝑟𝑏,𝑖𝑚,𝑗 𝜌𝑏=𝑌𝜇 𝑐𝐷,𝑖𝑚,𝑗 𝐾𝐷+𝑐𝐷,𝑖𝑚,𝑗 𝑐𝐴,𝑖𝑚,𝑗 𝐾𝐴+𝑐𝐴,𝑖𝑚,𝑗−𝑑 (4.11) 4.3 Simulation methods The simulation follows a time marching scheme. At each time step, it mainly consists of two steps: reactive transport simulation and model update caused by biofilm growth. 4.3.1 Reactive transport We use the finite element method (FEM) to discretize the governing equations (4.4) and (4.5). Then we apply the Newton-Raphson method to solve the global non-linear equations (Wang et al., 2021, Chapter2). The reactive transport simulation is implemented by the code RT_MRMT_DSA developed by Wang et al. (2021, Chapter 2), which is coupled to biofilm growth by updating the hydrodynamic properties of porous media and the local residence time in the different portions of the biofilm. These two simulations are run sequentially. The reactive transport simulation solves the state variables of 𝑐𝑚, 𝑐𝑖𝑚,𝑗, 𝑟𝑘,𝑚, 𝑟𝑘,𝑖𝑚,𝑗 and 𝑟𝑏,𝑖𝑚,𝑗 at each mesh node. 4.3.2 Model update At each time step, the growth of the volume fraction of biofilm is updated according to equation (4.7), that is ∆𝜙𝑏,𝑖𝑚,𝑗=𝜙𝑏,𝑖𝑚,𝑗(𝑡)𝑟𝑏,𝑖𝑚,𝑗(𝑡) 𝜌𝑏∆𝑡 𝑗=1,⋯,𝑁 (4.12) The increase of the volume fraction of the 𝑗th immobile water is the contribution of all microbial species in which, ∆𝜙𝑖𝑚,𝑗=𝜃𝑏∑∆𝜙𝑏,𝑖𝑚,𝑗 𝑁𝑏 𝑏=1 𝑗=1,⋯,𝑁 (4.13) where 𝑁𝑏 is the number of microbial species. At each mesh node, the increase of the total volume fraction of immobile water is the sum of all the fractions of the volume of immobile water, which is ∆𝜙𝑖𝑚=∑∆𝜙𝑖𝑚,𝑗 𝑁 𝑗=1 (4.14) The total volume fraction of the mobile water and biofilm in porous media is a constant which equals to CHAPTER 4 42 Figure 4.8 Breakthrough curves of solute at x = 5cm for continuous injection at t=13, 62 and 113 days. Figure 4.9 The evolution of residence times and probabilities in different portion of immobile zones during biofilm growth at x = 5cm and t = 13, 62 and 113 days. 43 Figure 4.10 The distribution of reaction rates in space and immobile zones at t = 13, 62 and 113 days from top to bottom. 4.5 Conclusions We have developed a conceptual model for biofilm growth based on the non-local MRMT formulation, which takes into account the changes in transport characteristics. Clearly, biofilm growth affects transport. This is shown by both the conceptual model and the experiments we used to validate the models. Although the modeled results do not agree perfectly the experiments, the model can reproduce the change in transport behavior due to biofilm growth. We modeled two experiments, which revealed some interesting differences between the two. In the experiment of Kone et al. (2014) the abundance of electron acceptors and donors made the biofilm growth independent of transport. This caused an equal change for all immobile zone, that is, the probability density function 𝑝 only shifts, but does not change its form. On the other hand, in the experiment of Seifert and Engesgaard (2007) biofilm growth is dependent on transport behavior. This causes a change probability density function 𝑝 with a larger effect on the immobile zones with smaller residence time and a tendency to decrease the heterogeneity. For both models we used a constant shape parameter 𝛾=2/3. This leads to an exponential biofilm growth with time. This can be seen clearly for the experiment of Kone et al. (2014), but is also true for that of Seifert and Engesgaard (2007). Obviously, the shape parameter will increase, when biofilm volume fraction approach values equal to the porosity during clogging. In that case, models should consider changes in shape parameter. 45 5 Conclusions Each chapter of this thesis is independent and contains the relevant conclusions. Therefore, here I simply summarize the most salient conclusions. In chapter 2, I presented a novel method to solve MRMT reactive transport models. The method takes advantage of the sparse nature of the system matrix that results in MRMT problems. This allows solving immobile concentrations as a function of mobile concentrations, so that the final system is only a function of mobile concentrations. The method turns out to be extremely efficient. In chapter 3, We analyzed the interplay between characteristic reaction times, transport times and distribution of residence times. We conclude that reactions driven by species that are not present in the inflow water but are the result of previous reactions will take place in immobile zones, whose residence time is comparable to reaction times. Furthermore, immobile zones with residence times much smaller than those for transport can be lumped together (assuming that very fast reactions are assumed in equilibrium), which greatly reduces computations. Finally, a biofilm growth model was presented in Chapter 4. The singularity of this model lies in its ability to represent the changes in residence time distributions as a result of biofilm growth. The model has been successful in reproducing tracer test experiments extracted from the literature performed in the laboratory at different stages of biofilm growth. 47 A. Numerical discretization of governing equations In practice, MRMT model is substituted by a finite number of immobile zones. In this way, the continuous governing equation (2.7) and (2.8) are discretized as follows, 𝜙𝑚𝜕𝒖𝑚(𝑥,𝑡) 𝜕𝑡 =𝐿𝑡[𝒖𝑚(𝑥,𝑡)]−∑𝑭𝑗 𝑁 𝑗=1 +𝜙𝑚𝑹𝑚(𝑥,𝑡) (A.1) 𝑭𝑗=𝜙𝑖𝑚,𝑗𝛼𝑗[𝒖𝑚(𝑥,𝑡)−𝒖𝑖𝑚,𝑗(𝑥,𝑡)] (A.2) 𝑹𝑚(𝑥,𝑡)=𝑼𝑺𝑘 𝑇𝒓𝑘,𝑚(𝒖𝑚(𝑥,𝑡)) (A.3) 𝜕𝒖𝑖𝑚,𝑗(𝑥,𝑡) 𝜕𝑡 =𝛼𝑗[𝒖𝑚(𝑥,𝑡)−𝒖𝑖𝑚,𝑗(𝑥,𝑡)]+𝑹𝑖𝑚,𝑗(𝑥,𝑡), 𝑗=1,⋯,𝑁 (A.4) 𝑹𝑖𝑚,𝑗(𝑥,𝑡)=𝑼𝑺𝑘 𝑇𝒓𝑘,𝑖𝑚,𝑗(𝒖𝑖𝑚,𝑗(𝑥,𝑡)) (A.5) where 𝑭𝑗 is the mass exchange between mobile and 𝑗th immobile zone, 𝜙𝑖𝑚,𝑗=𝜙𝑖𝑚𝑝𝑗 is the weighted porosity in 𝑗th immobile zone, 𝑝𝑗 is the probability for mass exchange rate 𝛼𝑗 accounting in the 𝑗th immobile zone, 𝑁 is the number of immobile zones. To obtain the numerical discretization of the system equations, a forward finite difference method is used to discretize the first derivative in time, and the finite element method is applied to discretize the governing partial difference equations. Meanwhile, state variables are evaluated at some time between time step 𝑘 and 𝑘+1. The resulting discretized system equations are given below, (𝒈𝑚)𝑖=(𝑮 ∆𝑡+𝜃𝑡𝑬)(𝒖𝑚 𝑘+1)𝑖−(𝑮 ∆𝑡−(1−𝜃𝑡)𝑬)𝒖𝑚 𝑘+𝑮 𝜙𝑚∑(𝑭𝑗𝑘+𝜃𝑡)𝑖 𝑁 𝑗=1 −𝑮(𝑹𝑚 𝑘+𝜃𝑟)𝑖 =0 (A.6) (𝒈𝑖𝑚,𝑗)𝑖=(𝒖𝑖𝑚,𝑗 𝑘+1)𝑖−𝒖𝑖𝑚,𝑗 𝑘 ∆𝑡 −𝛼𝑗[(𝒖𝑚 𝑘+𝜃𝑡)𝑖−(𝒖𝑖𝑚,𝑗 𝑘+𝜃𝑡)𝑖]−(𝑹𝑖𝑚,𝑗 𝑘+𝜃𝑟)𝑖=0, 𝑗=1,⋯,𝑁 (A.7) in which, 𝑬 is the global matrix accounting for advection and dispersion, 𝑮 is the global matrix that assembles porosity in the mobile zone, 𝑖 is the iteration number at each time step, 𝜃𝑡∈[0,1] is a temporal weight factor for transport, and 𝜃𝑟∈[0,1] is a weight factor for kinetics. The entries of Jacobian matrix are formulated as follows, (𝜕𝒈𝑚 𝜕𝒖𝑚 𝑘+1)𝑖=(𝑮 ∆𝑡+𝑬𝜃)⊗𝑰𝑁𝑢+(𝑮/𝜙𝑚)∑𝜙𝑖𝑚,𝑗𝛼𝑗 𝑁 𝑗=1 𝜃𝑡𝑰𝑁𝑛⊗𝑰𝑁𝑢−𝑮𝜃𝑟(𝜕𝑹𝑚 𝑘+1 𝜕𝒖𝑚 𝑘+1)𝑖 (A.8) (𝜕𝒈𝑚 𝜕𝒖𝑖𝑚,𝑗 𝑘+1)𝑖=−𝜙𝑖𝑚,𝑗𝛼𝑗𝜃𝑡(𝑮/𝜙𝑚)𝑰𝑁𝑛⊗𝑰𝑁𝑢 (A.9) (𝜕𝒈𝑖𝑚,𝑗 𝜕𝒖𝑚 𝑘+1)𝑖=−𝛼𝑗 𝜃𝑡𝑰𝑁𝑛⊗𝑰𝑁𝑢 (A.10) APPENDIX A 48 (𝜕𝒈𝑖𝑚,𝑗 𝜕𝒖𝑖𝑚,𝑗 𝑘+1)𝑖=(1 ∆𝑡+𝛼𝑗𝜃𝑡)𝑰𝑁𝑛⊗𝑰𝑁𝑢−𝜃𝑟(𝜕𝑹𝑖𝑚,𝑗 𝑘+1 𝜕𝒖𝑖𝑚,𝑗 𝑘+1)𝑖 (A.11) in which, 𝑁𝑛 is the number of mesh nodes, 𝑁𝑢 is the number of components, ⊗ represents Kronecker product, 𝑰𝑁𝑛 and 𝑰𝑁𝑢 are identity matrix with dimensions equal to 𝑁𝑛×𝑁𝑛 and 𝑁𝑢×𝑁𝑢, respectively. The derivatives of reaction rates with respect to components at each node are calculated according to the chain rule, that is (𝜕𝑹𝑚 𝑘+1 𝜕𝒖𝑚 𝑘+1)𝑖=𝑼𝑺𝐾 𝑇(𝜕𝒓𝑘,𝑚 𝑘+1 𝜕𝒖𝑚 𝑘+1)𝑖=𝑼𝑺𝐾 𝑇(𝜕𝒓𝑘,𝑚 𝑘+1 𝜕𝒄𝟏𝑚 𝑘+1)𝑖[(𝜕𝒖𝑚 𝑘+1 𝜕𝒄𝟏𝑚 𝑘+1)𝑖]−1 (A.12) (𝜕𝑹𝑖𝑚,𝑗 𝑘+1 𝜕𝒖𝑖𝑚,𝑗 𝑘+1)𝑖=𝑼𝑺𝐾 𝑇(𝜕𝒓𝑘,𝑖𝑚,𝑗 𝑘+1 𝜕𝒖𝑖𝑚,𝑗 𝑘+1)𝑖=𝑼𝑺𝐾 𝑇(𝜕𝒓𝑘,𝑖𝑚,𝑗 𝑘+1 𝜕𝒄𝟏𝑖𝑚,𝑗 𝑘+1)𝑖[(𝜕𝒖𝑖𝑚,𝑗 𝑘+1 𝜕𝒄𝟏𝑖𝑚,𝑗 𝑘+1)𝑖]−1 (A.13) both are matrices of size 𝑁𝑢×𝑁𝑢. 49 B. Analytical solution of reactive transport in multicontinuum media for first-order kinetics in the Laplace Domain For a chemical system satisfies 𝑼𝑺𝑘 𝑇𝒓𝑘,𝑚(𝑥,𝑡)=−𝒓𝑘,𝑚(𝑥,𝑡), 𝑼𝑺𝑘 𝑇𝒓𝑘,𝑖𝑚(𝑥,𝛼,𝑡)=−𝒓𝑘,𝑖𝑚(𝑥,𝛼,𝑡), and the chemical kinetics follows first-order decay, 𝒓𝑘,𝑚(𝑥,𝑡)=𝜅𝒖𝑚(𝑥,𝑡), 𝒓𝑘,𝑖𝑚(𝑥,𝛼,𝑡)= 𝜅𝒖𝑖𝑚(𝑥,𝛼,𝑡), then the governing equation (2.7) and (2.8) in mobile and immobile domains simplify to be linear, that is, 𝜙𝑚𝜕𝒖𝑚(𝑥,𝑡) 𝜕𝑡 =𝐿𝑡[𝒖𝑚(𝑥,𝑡)]−𝜙𝑖𝑚∫ 𝛼[𝒖𝑚(𝑥,𝑡)−𝒖𝑖𝑚(𝑥,𝛼,𝑡)]𝑓(𝛼) ∞ 0𝑑𝛼 −𝜙𝑚𝜅𝒖𝑚(𝑥,𝑡) (B.1) 𝜕𝒖𝑖𝑚(𝑥,𝛼,𝑡) 𝜕𝑡 =𝛼[𝒖𝑚(𝑥,𝑡)−𝒖𝑖𝑚(𝑥,𝛼,𝑡)]−𝜅𝒖𝑖𝑚(𝑥,𝛼,𝑡) (B.2) Inserting equation (B.2) into (B.1) leads to the total governing equation, written as 𝜙𝑚𝜕𝒖𝑚(𝑥,𝑡) 𝜕𝑡 +𝜙𝑖𝑚∫ 𝑓(𝛼)𝜕𝒖𝑖𝑚(𝑥,𝛼,𝑡) 𝜕𝑡 ∞ 0𝑑𝛼=𝐿𝑡[𝒖𝑚(𝑥,𝑡)]−𝒓𝑘(𝑥,𝑡) (B.3) in which, 𝒓𝑘(𝑥,𝑡) is the total reaction rate that integrates reactions in both mobile and immobile zones, 𝒓𝑘(𝑥,𝑡)=𝜙𝑚𝜅𝒖𝑚(𝑥,𝑡)+𝜙𝑖𝑚∫ 𝑓(𝛼)𝜅𝒖𝑖𝑚(𝑥,𝛼,𝑡) ∞ 0𝑑𝛼 (B.4) Solving equation (B.2), we get the solution of components in the immobile zone, it is given by 𝒖𝑖𝑚(𝑥,𝛼,𝑡)=𝒖𝑖𝑚 0𝑒−𝜆𝑡+∫𝛼𝑒−𝜆(𝑡−𝜏)𝒖𝑚(𝑥,𝜏) 𝑡 0𝑑𝜏 =𝒖𝑖𝑚 0𝑒−𝜆𝑡+𝛼𝑒−𝜆𝑡∗𝒖𝑚(𝑥,𝑡) (B.5) in which, 𝒖𝑖𝑚 0 is the initial condition in immobile zones, 𝜆=𝛼+𝜅 is the decay rate that accounts for both mass exchange rate and kinetic rate. Assuming 𝒖𝑖𝑚 0=0, then substituting equation (B.5) into the total governing equation (B.3), we obtain the total governing equation with respect to component only in mobile zone, that is 𝜙𝑚𝜕𝒖𝑚(𝑥,𝑡) 𝜕𝑡 =𝐿𝑡[𝒖𝑚(𝑥,𝑡)]−𝜙𝑖𝑚[𝑔0𝒖𝑚(𝑥,𝑡)+𝜕𝑔(𝑡) 𝜕𝑡 ∗𝒖𝑚(𝑥,𝑡)+𝜅𝑔(𝑡)∗𝒖𝑚(𝑥,𝑡)] −𝜙𝑚𝜅𝒖𝑚(𝑥,𝑡) (B.6) in which, 𝑔(𝑡)=∫𝑓(𝛼) ∞ 0𝛼𝑒−𝜆𝑡𝑑𝛼 represents the memory function (Carrera et al., 1998). To characterize the total governing equation (B.6), we define characteristic length and characteristic transport time, written as 𝐿𝑐=𝜙𝑚𝐷 𝑞 (B.7) and APPENDIX B 50 𝑡𝑐=𝐿𝑐 𝑣=𝜙𝑚 2𝐷 𝑞2 (B.8) In one dimensional, the characteristic length is the longitudinal dispersivity, i.e. 𝐿𝑐=𝛼𝐿, due to the mechanical dispersion equals to longitudinal dispersivity multiply fluid velocity, i.e. 𝐷=𝛼𝐿𝑣. Introducing these characteristics into equation (B.6), the dimensionless form of the total governing equation (B.6) is obtained, it is given by 𝜕𝒖𝑚(𝑥𝐷,𝑡𝐷) 𝜕𝑡𝐷=𝜕2𝒖𝑚(𝑥𝐷,𝑡𝐷) 𝜕𝑥𝐷 2−𝜕𝒖𝑚(𝑥𝐷,𝑡𝐷) 𝜕𝑥𝐷 −𝜂[𝑔0𝒖𝑚(𝑥𝐷,𝑡𝐷)+𝜕𝑔(𝑡𝐷) 𝜕𝑡𝐷∗𝒖𝑚(𝑥𝐷,𝑡𝐷)+𝜅𝐷𝑔(𝑡𝐷)∗𝒖𝑚(𝑥𝐷,𝑡𝐷)] −𝜅𝐷𝒖𝑚(𝑥𝐷,𝑡𝐷) (B.9) with the definition of dimensionless variables as 𝜂=𝜙𝑖𝑚 𝜙𝑚, 𝑡𝐷=𝑡 𝑡𝑐, 𝑥𝐷=𝑥 𝐿𝑐, 𝜅𝐷=𝜅𝑡𝑐, 𝛼𝐷=𝛼𝑡𝑐, 𝜆𝐷=𝜆𝑡𝑐 (B.10) Given the initial and boundary conditions, 𝒖𝑚(𝑥𝐷,𝑡𝐷=0)=0, 𝑥𝐷≥0 𝒖𝑚(𝑥𝐷=0,𝑡𝐷)=𝒖 , 𝒖𝑚(𝑥𝐷=∞,𝑡𝐷)=0, 𝑡𝐷>0 (B.11) The solution of equation (B.9) is obtained in the Laplace domain, that is ℒ{𝒖𝑚}=𝒖 𝑠𝑒𝑥𝑝{[1−√1+4(𝑠+𝜅𝐷)(1+𝜂ℒ{𝑔})]𝑥𝐷 2} (B.12) in which, ℒ{𝒖𝑚} and ℒ{𝑔} indicate the Laplace transform of 𝒖𝑚(𝑥𝐷,𝑡𝐷) and 𝑔(𝑡𝐷), respectively. 51 C. Breakthrough curves of species A, B and C for sequential reactions Figure C.1 Breakthrough curves of species A, B and C for 𝑞𝑚=1.0 [LT-1], thus 𝑡𝑐=0.1 [T-1], 𝜏𝑟1𝐷=𝜏𝑟2𝐷=102, 103 and 104 from top to bottom. 59 Bibliography Adams, E. E., & Gelhar, L. W. (1992). Field study of dispersion in a heterogeneous aquifer: 2. Spatial moments analysis. Water Resources Research. https://doi.org/10.1029/92WR01757 Babaei, M., & Islam, A. (2018). Convective-Reactive CO2 Dissolution in Aquifers With Mass Transfer With Immobile Water. Water Resources Research. https://doi.org/10.1029/2018WR023150 Babey, T., de Dreuzy, J. R., & Casenave, C. (2015). Multi-Rate Mass Transfer (MRMT) models for general diffusive porosity structures. Advances in Water Resources. https://doi.org/10.1016/j.advwatres.2014.12.006 Bakke, R., Trulear, M. G., Robinson, J. A., & Characklis, W. G. (1984). Activity of Pseudomonas aeruginosa in biofilms: Steady state. Biotechnology and Bioengineering, 26(12). https://doi.org/10.1002/bit.260261204 Bea, S. A., Carrera, J., Ayora, C., Batlle, F., & Saaltink, M. W. (2009). CHEPROO: A Fortran 90 objectoriented module to solve chemical processes in Earth Science models. Computers and Geosciences. https://doi.org/10.1016/j.cageo.2008.08.010 Beisman, J. J., Maxwell, R. M., Navarre-Sitchler, A. K., Steefel, C. I., & Molins, S. (2015). ParCrunchFlow: an efficient, parallel reactive transport simulation tool for physically and chemically heterogeneous saturated subsurface environments. Computational Geosciences. https://doi.org/10.1007/s10596-015-9475-x Benson, D. A., & Meerschaert, M. M. (2009). A simple and efficient random walk solution of multirate mobile/immobile mass transport equations. Advances in Water Resources. https://doi.org/10.1016/j.advwatres.2009.01.002 Benson, D. A., Wheatcraft, S. W., & Meerschaert, M. M. (2000). The fractional-order governing equation of Levy motion. Water Resources Research. https://doi.org/10.1029/2000WR900032 Berkowitz, B., Cortis, A., Dentz, M., & Scher, H. (2006). Modeling Non-fickian transport in geological formations as a continuous time random walk. Reviews of Geophysics. https://doi.org/10.1029/2005RG000178 Berkowitz, B., Emmanuel, S., & Scher, H. (2008). Non-Fickian transport and multiple-rate mass transfer in porous media. Water Resources Research. https://doi.org/10.1029/2007WR005906 Berkowitz, B., & Scher, H. (1998). Theory of anomalous chemical transport in random fracture networks. Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics. https://doi.org/10.1103/PhysRevE.57.5858 Berkowitz, B., & Scher, H. (2009). Exploring the nature of non-Fickian transport in laboratory experiments. Advances in Water Resources. https://doi.org/10.1016/j.advwatres.2008.05.004 Bilke, L., Flemisch, B., Kalbacher, T., Kolditz, O., Helmig, R., & Nagel, T. (2019). Development of Open-Source Porous Media Simulators: Principles and Experiences. Transport in Porous Media. https://doi.org/10.1007/s11242-019-01310-1 Brangarí, A. C., Fernàndez-Garcia, D., Sanchez-Vila, X., & Manzoni, S. (2018). Ecological and soil hydraulic implications of microbial responses to stress – A modeling analysis. Advances in Water Resources, 116. https://doi.org/10.1016/j.advwatres.2017.11.005 Carles Brangarí, A., Sanchez-Vila, X., Freixa, A., M. Romaní, A., Rubol, S., & Fernàndez-Garcia, D. (2017). A mechanistic model (BCC-PSSICO) to predict changes in the hydraulic properties for bio-amended variably saturated soils. Water Resources Research, 53(1). https://doi.org/10.1002/2015WR018517 BIBLIOGRAPHY 60 Carrel, M., Morales, V. L., Dentz, M., Derlon, N., Morgenroth, E., & Holzner, M. (2018). Pore-Scale Hydrodynamics in a Progressively Bioclogged Three-Dimensional Porous Medium: 3-D Particle Tracking Experiments and Stochastic Transport Modeling. Water Resources Research, 54(3). https://doi.org/10.1002/2017WR021726 Carrel, Maxence, Morales, V. L., Beltran, M. A., Derlon, N., Kaufmann, R., Morgenroth, E., & Holzner, M. (2018). Biofilms in 3D porous media: Delineating the influence of the pore network geometry, flow and mass transfer on biofilm development. Water Research, 134. https://doi.org/10.1016/j.watres.2018.01.059 Carrera, J., Sánchez-Vila, X., Benet, I., Medina, A., Galarza, G., & Guinerà, J. (1998). On matrix diffusion: Formulations, solution methods and qualitative effects. Hydrogeology Journal. https://doi.org/10.1007/s100400050143 Chen-Charpentier, B. (1999). Numerical simulation of biofilm growth in porous media. Journal of Computational and Applied Mathematics. https://doi.org/10.1016/S0377-0427(98)00240-4 Christensen, T. H., Kjeldsen, P., Bjerg, P. L., Jensen, D. L., Christensen, J. B., Baun, A., Albrechtsen, H. J., & Heron, G. (2001). Biogeochemistry of landfill leachate plumes. In Applied Geochemistry (Vol. 16, Issues 7–8). https://doi.org/10.1016/S0883-2927(00)00082-2 Cirpka, O. A., Frind, E. O., & Helmig, R. (1999). Numerical simulation of biodegradation controlled by transverse mixing. Journal of Contaminant Hydrology, 40(2). https://doi.org/10.1016/S01697722(99)00044-3 Cortis, A., & Berkowitz, B. (2004). Anomalous Transport in “Classical” Soil and Sand Columns. Soil Science Society of America Journal. https://doi.org/10.2136/sssaj2004.1539 Costerton, J. W., Lewandowski, Z., Caldwell, D. E., Korber, D. R., & Lappin-Scott, H. M. (1995). Microbial biofilms. In Annual Review of Microbiology. https://doi.org/10.1146/annurev.mi.49.100195.003431 Cunningham, A. B., Characklls, W. G., Abedeen, F., & Crawford, D. (1991). Influence of Biofilm Accumulation on Porous Media Hydrodynamics. Environmental Science and Technology. https://doi.org/10.1021/es00019a013 Cunningham, A. B., Sharp, R. R., Hiebert, R., & James, G. (2003). Subsurface biofilm barriers for the containment and remediation of contaminated groundwater. Bioremediation Journal. https://doi.org/10.1080/713607982 De Dreuzy, J. R., Rapaport, A., Babey, T., & Harmand, J. (2013). Influence of porosity structures on mixing-induced reactivity at chemical equilibrium in mobile/immobile Multi-Rate Mass Transfer (MRMT) and Multiple INteracting Continua (MINC) models. Water Resources Research. https://doi.org/10.1002/2013WR013808 De Dreuzy, Jean Raynald, & Carrera, J. (2016). On the validity of effective formulations for transport through heterogeneous porous media. Hydrology and Earth System Sciences. https://doi.org/10.5194/hess-20-1319-2016 Deng, H., & Spycher, N. (2019). Modeling reactive transport processes in fractures. Reviews in Mineralogy and Geochemistry, 85(1). https://doi.org/10.2138/rmg.2019.85.3 Deng, W., Cardenas, M. B., Kirk, M. F., Altman, S. J., & Bennett, P. C. (2013). Effect of permeable biofilm on micro-and macro-scale flow and transport in bioclogged pores. Environmental Science and Technology. https://doi.org/10.1021/es402596v Dentz, M., & Berkowitz, B. (2003). Transport behavior of a passive solute in continuous time random walks and multirate mass transfer. Water Resources Research. https://doi.org/10.1029/2001WR001163 Dentz, M., Cortis, A., Scher, H., & Berkowitz, B. (2004). Time behavior of solute transport in heterogeneous media: Transition from anomalous to normal transport. Advances in Water Resources. https://doi.org/10.1016/j.advwatres.2003.11.002 61 Dentz, M., Gouze, P., & Carrera, J. (2011). Effective non-local reaction kinetics for transport in physically and chemically heterogeneous media. Journal of Contaminant Hydrology. https://doi.org/10.1016/j.jconhyd.2010.06.002 Dentz, M., Kang, P. K., & Le Borgne, T. (2015). Continuous time random walks for non-local radial solute transport. Advances in Water Resources. https://doi.org/10.1016/j.advwatres.2015.04.005 Donado, L. D., Sanchez-Vila, X., Dentz, M., Carrera, J., & Bolster, D. (2009). Multicomponent reactive transport in multicontinuum media. Water Resources Research. https://doi.org/10.1029/2008WR006823 Eberl, H. J., Picioreanu, C., Heijnen, J. J., & Van Loosdrecht, M. C. M. (2000). Three-dimensional numerical study on the correlation of spatial structure, hydrodynamic conditions, and mass transfer and conversion in biofilms. Chemical Engineering Science. https://doi.org/10.1016/S0009-2509(00)00169-X Esperanza Cortés, M., Consuegra Bonilla, J., & Dario Sinisterra, R. (2011). Biofilm formation, control and novel strategies for eradication. In Science against microbial pathogens: communicating current research and technological advances. Fernàndez-Garcia, D., & Sanchez-Vila, X. (2015). Mathematical equivalence between time-dependent single-rate and multirate mass transfer models. Water Resources Research. https://doi.org/10.1002/2014WR016348 Flemming, H. C., Neu, T. R., & Wozniak, D. J. (2007). The EPS matrix: The “House of Biofilm Cells.” In Journal of Bacteriology. https://doi.org/10.1128/JB.00858-07 Flemming, H. C., & Wingender, J. (2010). The biofilm matrix. In Nature Reviews Microbiology. https://doi.org/10.1038/nrmicro2415 Flemming, H. C., & Wuertz, S. (2019). Bacteria and archaea on Earth and their abundance in biofilms. Nature Reviews Microbiology, 17(4). https://doi.org/10.1038/s41579-019-0158-9 Gaebler, H. J., & Eberl, H. J. (2018). A simple model of biofilm growth in a porous medium that accounts for detachment and attachment of suspended biomass and their contribution to substrate degradation. European Journal of Applied Mathematics, 29(6). https://doi.org/10.1017/S0956792518000189 Gouze, P., Le Borgne, T., Leprovost, R., Lods, G., Poidras, T., & Pezard, P. (2008). Non-Fickian dispersion in porous media: 1. Multiscale measurements using single-well injection withdrawal tracer tests. Water Resources Research. https://doi.org/10.1029/2007WR006278 Gouze, Philippe, Melean, Y., Le Borgne, T., Dentz, M., & Carrera, J. (2008). Non-Fickian dispersion in porous media explained by heterogeneous microscale matrix diffusion. Water Resources Research. https://doi.org/10.1029/2007WR006690 Guimerà, J., & Carrera, J. (2000). A comparison of hydraulic and transport parameters measured in lowpermeability fractured media. Journal of Contaminant Hydrology, 41(3–4), 261–281. https://doi.org/10.1016/S0169-7722(99)00080-7 Haggerty, R., & Gorelick, S. M. (1995). Multiple‐Rate Mass Transfer for Modeling Diffusion and Surface Reactions in Media with Pore‐Scale Heterogeneity. Water Resources Research. https://doi.org/10.1029/95WR10583 Haggerty, R., Harvey, C. F., Von Schwerin, C. F., & Meigs, L. C. (2004). What controls the apparent timescale of solute mass transfer in aquifers and soils? A comparison of experimental results. Water Resources Research. https://doi.org/10.1029/2002WR001716 Haggerty, R., McKenna, S. A., & Meigs, L. C. (2000). On the late-time behavior of tracer test breakthrough curves. Water Resources Research. https://doi.org/10.1029/2000WR900214 Haggerty, R., Wondzell, S. M., & Johnson, M. A. (2002). Power-law residence time distribution in the hyporheic zone of a 2nd-order mountain stream. Geophysical Research Letters. BIBLIOGRAPHY 62 https://doi.org/10.1029/2002GL014743 Hammond, G. E., Lichtner, P. C., & Mills, R. T. (2014). Evaluating the performance of parallel subsurface simulators: An illustrative example with PFLOTRAN. Water Resources Research. https://doi.org/10.1002/2012WR013483 Iraola, A., Trinchero, P., Karra, S., & Molinero, J. (2019). Assessing dual continuum method for multicomponent reactive transport. Computers and Geosciences. https://doi.org/10.1016/j.cageo.2019.05.007 Kapellos, G. E., Alexiou, T. S., & Pavlou, S. (2015). Fluid-Biofilm Interactions in Porous Media. In Heat Transfer and Fluid Flow in Biological Processes. https://doi.org/10.1016/B978-0-12408077-5.00008-0 Kim, J. W., Choi, H., & Pachepsky, Y. A. (2010). Biofilm morphology as related to the porous media clogging. Water Research. https://doi.org/10.1016/j.watres.2009.05.049 Kitanidis, P. K. (1988). Prediction by the method of moments of transport in a heterogeneous formation. Journal of Hydrology. https://doi.org/10.1016/0022-1694(88)90111-4 Kolditz, O., Bauer, S., Bilke, L., Böttcher, N., Delfs, J. O., Fischer, T., Görke, U. J., Kalbacher, T., Kosakowski, G., McDermott, C. I., Park, C. H., Radu, F., Rink, K., Shao, H., Shao, H. B., Sun, F., Sun, Y. Y., Singh, A. K., Taron, J., … Zehner, B. (2012). OpenGeoSys: An open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (THM/C) processes in porous media. Environmental Earth Sciences. https://doi.org/10.1007/s12665-012-1546-x Kolditz, Olaf, Görke, U.-J., Shao, H., & Wang, W. (2012). Thermo-Hydro-Mechanical-Chemical Processes in Porous Media: Benchmarks and Examples. In Springer. Kone, T., Golfier, F., Orgogozo, L., Oltéan, C., Lefèvre, E., Block, J. C., & Buès, M. A. (2014). Impact of biofilm-induced heterogeneities on solute transport in porous media. Water Resources Research. https://doi.org/10.1002/2013WR015213 Kosakowski, G., Berkowitz, B., & Scher, H. (2001). Analysis of field observations of tracer transport in a fractured till. Journal of Contaminant Hydrology. https://doi.org/10.1016/S01697722(00)00140-6 Le Borgne, T., & Gouze, P. (2008). Non-Fickian dispersion in porous media: 2. Model validation from measurements at different scales. Water Resources Research. https://doi.org/10.1029/2007WR006279 Le Borgne, Tanguy, Dentz, M., & Carrera, J. (2008). Lagrangian statistical model for transport in highly heterogeneous velocity fields. Physical Review Letters, 101(9). https://doi.org/10.1103/PhysRevLett.101.090601 Levy, M., & Berkowitz, B. (2003). Measurement and analysis of non-Fickian dispersion in heterogeneous porous media. Journal of Contaminant Hydrology. https://doi.org/10.1016/S01697722(02)00204-8 Lewandowski, Z., & Boltz, J. P. (2011). Biofilms in Water and Wastewater Treatment. In Treatise on Water Science. https://doi.org/10.1016/B978-0-444-53199-5.00095-6 Liu, P., Zhang, T., & Sun, S. (2019). A tutorial review of reactive transport modeling and risk assessment for geologic CO 2 sequestration. In Computers and Geosciences. https://doi.org/10.1016/j.cageo.2019.02.007 Lopez-Peña, L. A., Meulenbroek, B., & Vermolen, F. (2019). A network model for the biofilm growth in porous media and its effects on permeability and porosity. Computing and Visualization in Science, 21(1–6). https://doi.org/10.1007/s00791-019-00316-y Marseguerra, M., & Zoia, A. (2008). Monte Carlo evaluation of FADE approach to anomalous kinetics. Mathematics and Computers in Simulation. https://doi.org/10.1016/j.matcom.2007.03.001 Meysman, F. J. R., Middelburg, J. J., Herman, P. M. J., & Heip, C. H. R. (2003a). Reactive transport in 63 surface sediments. I. Model complexity and software quality. Computers and Geosciences. https://doi.org/10.1016/S0098-3004(03)00006-2 Meysman, F. J. R., Middelburg, J. J., Herman, P. M. J., & Heip, C. H. R. (2003b). Reactive transport in surface sediments. II. Media: an object-oriented problem-solving environment for early diagenesis. Computers & Geosciences. https://doi.org/10.1016/s0098-3004(03)00007-4 Molins, S., Carrera, J., Ayora, C., & Saaltink, M. W. (2004). A formulation for decoupling components in reactive transport problems. Water Resources Research. https://doi.org/10.1029/2003WR002970 Morales, V. L., Parlange, J. Y., & Steenhuis, T. S. (2010). Are preferential flow paths perpetuated by microbial activity in the soil matrix? A review. Journal of Hydrology, 393(1–2). https://doi.org/10.1016/j.jhydrol.2009.12.048 Neuman, S. P., & Tartakovsky, D. M. (2009). Perspective on theories of non-Fickian transport in heterogeneous media. Advances in Water Resources. https://doi.org/10.1016/j.advwatres.2008.08.005 Nicolella, C., Van Loosdrecht, M. C. M., & Heijnen, J. J. (2000). Wastewater treatment with particulate biofilm reactors. In Journal of Biotechnology. https://doi.org/10.1016/S0168-1656(00)00229-7 Nordström, A., Hellman, M., Hallin, S., & Herbert, R. B. (2021). Microbial controls on net production of nitrous oxide in a denitrifying woodchip bioreactor. Journal of Environmental Quality, 50(1). https://doi.org/10.1002/jeq2.20181 Nordström, A., & Herbert, R. B. (2017). Denitrification in a low-temperature bioreactor system at two different hydraulic residence times: laboratory column studies. Environmental Technology (United Kingdom). https://doi.org/10.1080/09593330.2016.1228699 Nordström, A., & Herbert, R. B. (2018). Determination of major biogeochemical processes in a denitrifying woodchip bioreactor for treating mine drainage. Ecological Engineering, 110. https://doi.org/10.1016/j.ecoleng.2017.09.018 Nordström, A., & Herbert, R. B. (2019). Identification of the temporal control on nitrate removal rate variability in a denitrifying woodchip bioreactor. Ecological Engineering, 127. https://doi.org/10.1016/j.ecoleng.2018.11.015 Parkhurst, D. L., & Appelo, C. A. J. (2013). Description of Input and Examples for PHREEQC Version 3 — A Computer Program for Speciation , Batch-Reaction , One-Dimensional Transport , and Inverse Geochemical Calculations. In U.S. Geological Survey Techniques and Methods, book 6, chapter A43. Picioreanu, C., Kreft, J. U., & Van Loosdrecht, M. C. M. (2004). Particle-based multidimensional multispecies biofilm model. Applied and Environmental Microbiology. https://doi.org/10.1128/AEM.70.5.3024-3040.2004 Picioreanu, C., Van Loosdrecht, M. C. M., & Heijnen, J. J. (1998). Mathematical modeling of biofilm structure with a hybrid differentialdiscrete cellular automaton approach. Biotechnology and Bioengineering, 58(1), 101–116. https://doi.org/10.1002/(SICI)10970290(19980405)58:1<101::AID-BIT11>3.0.CO;2-M Prommer, H., Barry, D. A., & Zheng, C. (2003). MODFLOW/MT3DMS-based reactive multicomponent transport modeling. Ground Water. https://doi.org/10.1111/j.17456584.2003.tb02588.x Prommer, H., Barry, D., & Zheng, C. (2001). PHT3D-A MODFLOW/MT3DMS based reactive multicomponent transport model. MODFLOW 2001 and Other Modeling Odysseys. Pruess, K., & Narasimhan, T. N. (1985). PRACTICAL METHOD FOR MODELING FLUID AND HEAT FLOW IN FRACTURED POROUS MEDIA. Society of Petroleum Engineers Journal. Raiders, R. A., Knapp, R. M., & McInerney, M. J. (1989). Microbial selective plugging and enhanced BIBLIOGRAPHY 64 oil recovery. Journal of Industrial Microbiology. https://doi.org/10.1007/BF01574079 Rapaport, A., Rojas-Palma, A., De Dreuzy, J. R., & Ramirez, H. C. (2017). Equivalence of Finite Dimensional Input-Output Models of Solute Transport and Diffusion in Geosciences. IEEE Transactions on Automatic Control. https://doi.org/10.1109/TAC.2017.2701150 Rittmann, B. E., & McCarty, P. L. (2001). Environmental biotechnology : principles and applications. Current Opinion in Biotechnology, 7(3). Rodríguez-Escales, P., Folch, A., van Breukelen, B. M., Vidal-Gavilan, G., & Sanchez-Vila, X. (2016). Modeling long term Enhanced in situ Biodenitrification and induced heterogeneity in column experiments under different feeding strategies. Journal of Hydrology, 538. https://doi.org/10.1016/j.jhydrol.2016.04.012 Saaltink, M. W., Batlle, F., Ayora, C., Carrera, J., & Olivella, S. (2004). RETRASO, a code for modeling reactive transport in saturated and unsaturated porous media. Geologica Acta. https://doi.org/10.1344/105.000001430 Saaltink, Maarten W., Ayora, C., & Carrera, J. (1998). A mathematical formulation for reactive transport that eliminates mineral concentrations. Water Resources Research. https://doi.org/10.1029/98WR00552 Saaltink, Maarten W., Carrera, J., & Ayora, C. (2001). On the behavior of approaches to simulate reactive transport. Journal of Contaminant Hydrology. https://doi.org/10.1016/S01697722(00)00172-8 Salamon, P., Fernàndez-Garcia, D., & Gómez-Hernández, J. J. (2006). Modeling mass transfer processes using random walk particle tracking. Water Resources Research. https://doi.org/10.1029/2006WR004927 Satpathy, S., Sen, S. K., Pattanaik, S., & Raut, S. (2016). Review on bacterial biofilm: An universal cause of contamination. In Biocatalysis and Agricultural Biotechnology. https://doi.org/10.1016/j.bcab.2016.05.002 Schumer, R., Benson, D. A., Meerschaert, M. M., & Baeumer, B. (2003). Fractal mobile/immobile solute transport. Water Resources Research. https://doi.org/10.1029/2003WR002141 Seifert, D., & Engesgaard, P. (2007). Use of tracer tests to investigate changes in flow and transport properties due to bioclogging of porous media. Journal of Contaminant Hydrology, 93(1–4). https://doi.org/10.1016/j.jconhyd.2007.01.014 Sethi, R., & Di Molfetta, A. (2019). Remediation of Contaminated Groundwater. In Springer Tracts in Civil Engineering. https://doi.org/10.1007/978-3-030-20516-4_17 Shafahi, M., & Vafai, K. (2009). Biofilm affected characteristics of porous structures. International Journal of Heat and Mass Transfer. https://doi.org/10.1016/j.ijheatmasstransfer.2008.07.013 Silva, O., Carrera, J., Dentz, M., Kumar, S., Alcolea, A., & Willmann, M. (2009). A general real-time formulation for multi-rate mass transfer problems. Hydrology and Earth System Sciences. https://doi.org/10.5194/hess-13-1399-2009 Soler-Sagarra, J., Luquot, L., Martínez-Pérez, L., Saaltink, M. W., De Gaspari, F., & Carrera, J. (2016). Simulation of chemical reaction localization using a multi-porosity reactive transport approach. International Journal of Greenhouse Gas Control. https://doi.org/10.1016/j.ijggc.2016.01.026 Steefel, C. I., Appelo, C. A. J., Arora, B., Jacques, D., Kalbacher, T., Kolditz, O., Lagneau, V., Lichtner, P. C., Mayer, K. U., Meeussen, J. C. L., Molins, S., Moulton, D., Shao, H., Šimůnek, J., Spycher, N., Yabusaki, S. B., & Yeh, G. T. (2015). Reactive transport codes for subsurface environmental simulation. Computational Geosciences. https://doi.org/10.1007/s10596-014-9443-x Steefel, C. I., & Lasaga, A. C. (1994). A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive flow in single phase hydrothermal systems. American Journal of Science. https://doi.org/10.2475/ajs.294.5.529 65 Steefel, Carl I., & MacQuarrie, K. T. B. (1996). Approaches to modeling of reactive transport in porous media. Reviews in Mineralogy, 34. Taylor, S. W., & Jaffé, P. R. (1990a). Biofilm growth and the related changes in the physical properties of a porous medium: 1. Experimental investigation. Water Resources Research. https://doi.org/10.1029/WR026i009p02153 Taylor, S. W., & Jaffé, P. R. (1990b). Biofilm growth and the related changes in the physical properties of a porous medium: 3. Dispersivity and model verification. Water Resources Research. https://doi.org/10.1029/WR026i009p02171 Taylor, S. W., Milly, P. C. D., & Jaffé, P. R. (1990). Biofilm growth and the related changes in the physical properties of a porous medium: 2. Permeability. Water Resources Research. https://doi.org/10.1029/WR026i009p02161 Thullner, M. (2010). Comparison of bioclogging effects in saturated porous media within oneand twodimensional flow systems. Ecological Engineering. https://doi.org/10.1016/j.ecoleng.2008.12.037 Tiwari, S. K., & Bowers, K. L. (2001). Modeling biofilm growth for porous media applications. Mathematical and Computer Modelling, 33(1–3). https://doi.org/10.1016/S0895-7177(00)002466 Valhondo, C., Carrera, J., Martínez-Landa, L., Wang, J., Amalfitano, S., Levantesi, C., & Diaz-Cruz, M. S. (2020). Reactive barriers for renaturalization of reclaimed water during soil aquifer treatment. Water (Switzerland). https://doi.org/10.3390/W12041012 Valhondo, C., Martinez-Landa, L., Carrera, J., Ayora, C., Nödler, K., & Licha, T. (2018). Evaluation of EOC removal processes during artificial recharge through a reactive barrier. Science of the Total Environment, 612. https://doi.org/10.1016/j.scitotenv.2017.08.054 Valhondo, C., Martínez-Landa, L., Carrera, J., Díaz-Cruz, S. M., Amalfitano, S., & Levantesi, C. (2020). Six artificial recharge pilot replicates to gain insight into water quality enhancement processes. Chemosphere, 240. https://doi.org/10.1016/j.chemosphere.2019.124826 Valocchi, A. J. (1985). Validity of the Local Equilibrium Assumption for Modeling Sorbing Solute Transport Through Homogeneous Soils. Water Resources Research. https://doi.org/10.1029/WR021i006p00808 Van Hamme, J. D., Singh, A., & Ward, O. P. (2003). Recent Advances in Petroleum Microbiology. Microbiology and Molecular Biology Reviews. https://doi.org/10.1128/mmbr.67.4.503-549.2003 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. https://doi.org/10.1002/aic.11674 Wang, P. P., Zheng, C., & Gorelick, S. M. (2005). A general approach to advective-dispersive transport with multirate mass transfer. Advances in Water Resources. https://doi.org/10.1016/j.advwatres.2004.10.003 Willmann, M., Carrera, J., & Sánchez-Vila, X. (2008). Transport upscaling in heterogeneous aquifers: What physical parameters control memory functions? Water Resources Research. https://doi.org/10.1029/2007WR006531 Willmann, M., Carrera, J., Sanchez-Vila, X., Silva, O., & Dentz, M. (2010). Coupling of mass transfer and reactive transport for nonlinear reactions in heterogeneous media. Water Resources Research. https://doi.org/10.1029/2009WR007739 Wu, Y., Cai, P., Jing, X., Niu, X., Ji, D., Ashry, N. M., Gao, C., & Huang, Q. (2019). Soil biofilm formation enhances microbial community diversity and metabolic activity. Environment International, 132. https://doi.org/10.1016/j.envint.2019.105116 Zhang, Y., Benson, D. A., & Baeumer, B. (2007). Predicting the tails of breakthrough curves in BIBLIOGRAPHY 66 regional-scale alluvial systems. Ground Water. https://doi.org/10.1111/j.1745-6584.2007.00320.x Zinn, B., Meigs, L. C., Harvey, C. F., Haggerty, R., Peplinski, W. J., & Von Schwerin, C. F. (2004). Experimental visualization of solute transport and mass transfer processes in two-dimensional conductivity fields with connected regions of high conductivity. Environmental Science and Technology. https://doi.org/10.1021/es034958g