Reactive Transport in Evolving Porous Media

Reactive transport modeling is a process-based approach that accounts for advection, diffusion, dispersion and a multitude of biogeochemical reactions. The occurrence of these reactions, by nature, tends to affect the properties of porous media in many ways (Tenthorey and Gerald 2006). If these alteration reactions are significant, then feedback mechanisms could occur that influence the flow of groundwater as well as the migration of solutes and gases through porous media (Le Gallo et al. 1998; Kaszuba et al. 2005; Jin et al. 2013). In addition, changes induced by the reactions on the solid grains can also affect the rates of the reactions themselves (Hao et al. 2012; Harrison et al. 2017). A prime example for reactive transport in evolving porous media is the dissolution of mineral phases. If dissolution reactions are substantial, the porosity, i.e., the void space between grains or apertures of fractures in jointed rocks, will increase. Such an increase in porosity commonly has secondary effects, by altering the connectivity or larger scale pores in the porous medium under consideration (Navarre-Sitchler et al. 2009). Together, these changes in porosity and connectivity can substantially affect flow and transport processes by modifying the key transport parameters such as the medium's permeability and tortuosity, leading to alteration of the groundwater flow regime and modification of transport pathways. The impact of these changes can affect transport in the water phase as well as in the gas phase. In addition, because mineral dissolution reshapes the surface of the dissolving phases or leads to the complete dissolution of smaller particles, the system's reactivity can be affected as well, leading to a direct feedback on reaction progress and rates (Noiriel et al. 2009).

The dissolution of minerals is only one example for evolving porous media. Flow and transport properties can be affected by a multitude of other processes such as mineral precipitation, possibly leading to clogging (Gaucher and Blanc 2006; Brovelli et al. 2009; Chagneau et al. 2015) or modification of the pore size distribution (Emmanuel et al. 2010, 2015), microbial activity leading to the growth of biofilms (Kim and Fogler 2000; Ezeuko et al. 2011) or bioclogging (Thullner et al. 2002), swelling of clays (Herbert et al. 2008; Sedighi and Thomas 2014; Wang et al. 2014), and variations in surface loading and temperatures leading to the expansion or contraction of porous media and possibly the formation or closure of fractures (Pfingsten 2002; MacQuarrie and Mayer 2005; Tian et al. 2014). In this chapter, we will focus on evolving porous media as affected by biogeochemical reactions with an emphasis on mineral dissolution and precipitation.

Specific examples where biogeochemical reactions result in the evolution of porous media include wellbore integrity studies. It has been shown that the productivity of wells can decrease due to clogging. Houben (2003) showed that chemical processes including corrosion and mineral precipitation account for more than 90% of ageing of wells, leading to a decrease in their productivity. In addition, geochemical interactions at interfaces between cementitious materials and clay can lead to the formation of an impermeable layer of calcite and other secondary minerals, effectively blocking solute transfer across the interfaces (Dauzères et al. 2010, 2019). These processes are important in the context of the safety assessment for the long-term storage of spent nuclear fuel (Atkinson et al. 1987; Spycher et al. 2003; Soler and Mäder 2005; Yang et al. 2008). An additional example where evolution of porous media is of importance is the long-term evolution of mine waste deposits. In these systems, the formation of alteration rims and secondary mineral formation has a strong impact on reactivity, leading to a decline of acid generation and metal release over time due to surface passivation (Blowes and Jambor 1990; Wunderly et al. 1996; Johnson et al. 2000). In extreme situations, extensive secondary mineral formation can lead to the formation of hard-pans (Blowes et al. 1991), with a pronounced effect on transport parameters. Furthermore, systems designed for the treatment of contaminated groundwater are also often affected by evolving transport and reactivity parameters. The precipitation of secondary mineral phases leads to passivation of carbonate minerals and limited pH-buffering capacity in limestone drains for the treatment of acid rock drainage (Hammarstrom et al. 2003; Rötting et al. 2008). Groundwater treatment by permeable reactive barriers (Blowes et al. 2000) or by bioremediation (Hazen and Fliermans 1995; Ritter and Scarborough 1995) is frequently affected by the growth of biofilms, which can have detrimental effects on the efficacy of the remediation effort and modifies permeability and transport parameters (Scherer et al. 2000; Li et al. 2006; Wantanaphong et al. 2006). The precipitation of secondary oxide phases in case of groundwater remediation via in-situ chemical oxidation with permanganate also may result in less effective contaminant treatment or contaminant encapsulation (Nelson et al. 2001; Henderson et al. 2009). Furthermore, evolving porous media plays a role in the clogging and passivation of permeable reactive barriers used for the treatment of chlorinated solvents and metals, which occurs due to pronounced changes in pH and redox conditions causing secondary mineral formation (Ouellet et al. 2006; Jeen et al. 2007).

Accounting for the evolution of porous media is not only of importance in man-made and engineered systems, but is equally relevant in natural environments, where transport and chemical properties can exhibit significant variation in space and time (Gouze and Coudrain-Ribstein 2002; White et al. 2005; Jin et al. 2013; Opolot and Finke 2015). A key example is the progress of weathering and soil formation. Also, karst formation caused by dissolution reactions of soluble minerals (predominantly carbonate minerals and gypsum) induced by the ingress of slightly acidic and dilute recharge. Birk et al. (2003, 2005) reported that the evolution of heterogeneity patterns strongly depends on flow and weathering rates. Karst development is obtained by positive feedback between reaction and flow (Hartmann et al. 2014): higher flow rates, lead to higher dissolution rates, resulting in higher porosity and permeability, in turn further increasing flow rates (Xiao et al. 2007). Kaufmann et al. (2010) showed that in karst systems permeability can vary by up to 9 orders of magnitude between fractures and the solid matrix. The morphology of the karst systems strongly depends on dissolution kinetics (Raines and Dewers 1997).

In addition, precipitation reactions in natural systems can have significant consequences for porous media evolution. Cementation (Giles et al. 2000; Putnis and Mauthe 2001; Pape et al. 2005) can reduce porosity and permeability, while strengthening the soil mechanical properties. It has also been shown that soil hydraulic properties can be decreased due to biogeochemical reactions causing bioclogging or biocementation (Ivanov and Chu 2008). Sedimentary rocks are often cemented due to the precipitation of clay minerals on the solid grain surfaces (Peters 2009). Such precipitation reactions can also have an impact on the reactivity of the primary minerals. Peters (2009) showed that while secondary minerals only occupied a volume fraction of 5-30%, they represented 65-86% of the reactive surface. Precipitation of secondary minerals can also lead to the sealing of fractures (Noiriel et al. 2010; Zhang 2011, 2013; Phillips et al. 2016). A decrease of permeability is often observed in certain hydrothermally altered rocks induced by mineral precipitation (Fontaine et al. 2001; Dobson et al. 2003a; Polak et al. 2003; Giger et al. 2007). Mineral precipitation can also lead to the formation of skarn (Meinert et al. 2003) or calcrete (Wang et al. 1994). In addition, several authors reported a correlation between chemical weathering and landslide occurrence (Jaboyedoff et al. 2004; Watanabe et al. 2005; Opfergelt et al. 2006; Regmi et al. 2013), implying that mechanical properties are evolving as chemical weathering progresses.

Despite the widespread relevance of reactive transport processes leading to changes in porous media properties, parameters such as porosity, permeability, effective diffusion coefficients and reactivity have often been considered constant in time in previous modeling efforts. The most common exception to this assumption has been in the consideration of the impact of changes in water content (Millington and Quirk 1961). Modeling studies of reaction induced porosity and permeability change date from the early 1990s, with Steefel and Lasaga (1992) presenting what appears to be the first numerical study of the reactive infiltration instability discussed by Ortoleva et al. (1987). In a later study of reactive transport in hydrothermal systems, Steefel and Lasaga (1994) analyzed the impact of quartz precipitation in fractures occurring as fluids cooled. They found that, unlike the reactive infiltration instability or “wormholing” effect due to dissolution of the rock matrix wherein reactive flow becomes increasingly focused in a smaller number of flow pathways, quartz precipitation had the effect of diffusing the flow paths over an increasingly wider volume, to the point where a convection cell might even reorganize completely. In another modeling study, Steefel and Lichtner (1994) showed how reactions in the rock matrix bordering a fracture could lead to “armoring” of the fracture, or isolation of the fracture volume from the surrounding rock matrix. Modeling was also used to show that the delicate balance between cementation of the rock matrix and the fracture itself can determine the subsequent evolution of the entire fractured rock system (Steefel and Lichtner 1998). The effect of evolving pore connectivity and tortuosity and their effect on diffusivity was considered in the modeling of chemical weathering of basaltic fragments (Navarre-Sitchler et al. 2011). Evolving reactivity has been less commonly treated in reactive transport models, except for the relatively straightforward case of dissolving (disappearing) phases where the reactive surface area necessarily must go to zero. However, one modeling study based on reactive flow experiments with calcite conducted in capillary tubes showed that it was necessary to consider the formation of high surface area precipitates in order to capture the evolving reactivity of the system (Noiriel et al. 2012). In a study of CO₂ injection into reactive volcanogenic sands, it was shown the use of a grain size distribution could capture at least part of the evolution of the reactivity of the system (Beckingham et al. 2017).

The simplification of the assumption of constant transport and reactivity parameters in modeling is due in part to the inherent difficulty of describing the temporal evolution of these parameters due to changes in porous media structure and composition. Traditionally, reactive transport has been simulated using upscaled parameters, i.e., on the continuum or Representative Elementary Volume (REV) scale. On this scale, relationships have been developed that describe interdependencies between transport parameters and porosity (e.g., Akanni et al. 1987; Hommel et al. 2018) as well as evolving surface area and reactivity. However, these relationships are commonly empirical in nature and have been developed based on fitting to experimental data (e.g., Low 1981; Tomadakis and Sotirchos 1993), or are based on theoretical considerations (e.g., Petersen 1958).

In reality, geochemical reactions and the processes controlling the evolution of porous media occur at the pore scale. Given the complexity and heterogeneity of the pore structure, the wide range of chemical perturbations and the coupling between these different processes, obtaining a representative evolution of the effective properties of the porous medium in a continuum approach can be a challenge. With the recent advent of improved imaging techniques such as X-ray microtomography (e.g., Werth et al. 2010), the development of improved pore-scale model formulations (e.g., Navarre-Sitchler et al. 2009; Dewanckele et al. 2012; Beckingham et al. 2013; Vilcáez et al. 2017), and substantial advances in computational power, the exploration of processes on the pore scale has become possible from both an experimental and modeling perspective. As a result, recent studies have focused more on describing these phenomena at the pore-scale (Molins et al. 2012, 2014; Steefel et al. 2013; Trebotich et al. 2014; Molins 2015; Deng et al. 2018). Although many advances have been made in pore scale investigations in recent years, it is difficult to simulate reactive transport on larger scales (on the orders of m's to km's) using pore scale approaches, although recently developed hybrid multi-scale approaches might handle this problem. For such studies, it is necessary to continue relying on continuum approaches. However, pore scale studies are providing deeper understanding and insights, which in turn can be incorporated to advance continuum approaches.

In this chapter, we will briefly review the governing equations of multiphase flow and reactive transport on the continuum scale, emphasizing the parameters that are affected by biogeochemical reactions in evolving porous media, and defining them from a pore scale perspective. We will then review recent studies carried out on the pore scale to provide insights and fundamental understanding on processes that affect and control porous media parameters. Subsequently, we will return to the continuum scale to provide a summary of the common approaches used to describe evolving flow, transport and reaction parameters such as porosity, permeability, tortuosity and reactive surface area. In this context, we review different approaches taken by the community to take these phenomena into considerations. We will then follow up with selected reactive transport modeling case studies, accounting for evolving properties in porous media. Finally, we will provide an outlook on future work and opportunities to work towards an integrated model for describing reactive transport in evolving porous media and challenges associated with filling the gap between pore scale and continuum approaches.

In porous media, variably saturated flow is modeled by the coupled two-phase flow equations (Xu et al. 2011), described by Sin et al. (2017) as:

where t [s] is time, ϕ is porosity [–], α denotes the fluid phase, either liquid (l) or gas (g), S_α is the saturation of the fluid phase α [–], k_rα is the phase relative permeability [–]. p_α is the fluid pressure [kg m⁻¹ s⁻²], and Q_α is a source-sink term for the fluid phase α [s⁻¹]. ρ_α is the fluid density [kg m⁻³], g [m s⁻²] is the gravitational acceleration, μ_α [kg m⁻¹ s⁻¹] is the fluid dynamic viscosity and k is the intrinsic permeability tensor [m²]. However, some of the reactive transport simulators do not model the coupled two-phase flow (Steefel et al. 2015a), but instead use a simplified approach through Richards’ equation (Neuman 1973; Panday et al. 1993):

where h is the hydraulic head [m] and S_s defines the specific storage coefficient [m⁻¹]. The latter approach does not consider gas advection.

In variably saturated media, the reactive transport equations in global implicit form can be written as (Lichtner 1996; Mayer et al. 2002)

where q_l and q_g are the Darcy flux vectors, representing the specific discharges in both phases [m s⁻¹]. For single phase flow formulation using Richards' equation (Eqn. 2), the gas specific discharge q_g is not considered. $Tjl$ [mol L⁻¹ H₂O] and $Tjg$ [mol L⁻¹ gas] define the total aqueous/gaseous component concentrations for the j-th component. $Qjint$ [mol dm⁻³ porous medium] are internal source and sink terms due to all kinetically controlled reactions, including mineral dissolution–precipitation reactions. $Qjext$ [mol dm⁻³ porous medium] are external source and sink terms and define mass fluxes across the domain boundaries for the aqueous and gas phases and N_c defines the number of components. D_l and D_g are the hydrodynamic dispersion tensors [m² s⁻¹] for the aqueous and gas phases, respectively. The dispersion tensors include the contributions of molecular diffusion and mechanical dispersion, a process by which dissolved species migrate at different rates than the average linear groundwater velocity due the small-scale heterogeneities. In Equation (3), molecular diffusion is integrated into the hydrodynamic dispersion tensor in the form of the pore diffusion coefficient, defined as:

where $Dαp$ is the pore diffusion coefficient of phase α (aqueous or gas phase), τ_α is the tortuosity of phase α [–], and $Dα0$ is the free phase diffusion coefficient in phase α. The tortuosity accounts for diffusion along tortuous pathways through the porous medium due to the presence of grains and partial phase saturation, leading to a reduction of the pore diffusion coefficient in comparison to the free phase diffusion coefficient. Following previous definitions, tortuosity is here defined as a value smaller than unity (as in Steefel and Maher 2009), capturing that the presence of tortuous pathways leads to a decline in the pore diffusion coefficient. However, it should be noted that in other studies, an inverse definition has been used, in which case tortuosity is proportional to the length of the tortuous diffusion pathway in relation to the direct pathway, therefore giving values greater than unity (Walsh and Brace 1984; Hommel et al. 2018). Using this alternative formulation, Equation (4) would have to be modified, since for any formulation, the pore diffusion coefficient is smaller than the free phase diffusion coefficient. For completeness, it must also be noted that the effective diffusion coefficient $Dαe$⁠, representative for the rate of solute or gas diffusion through a porous medium, is obtained by multiplying the pore diffusion coefficient (Eqn. 4) with the phase saturation and porosity, which is accounted for in Equation (3).

In addition to the mass balance equations for the aqueous and gas phase (Eqn. 3), it is also necessary to define mass balance equations for reactive minerals (Steefel and Lasaga 1994):

where φ_i is the mineral volume fraction [dm³ mineral dm⁻³ porous medium] of a mineral phase, $Vim$ is the molar volume of the mineral [dm³ mineral mol⁻¹], $Rim$ is the mineral i overall dissolution–precipitation rate [mol dm⁻³ porous medium s⁻¹], and N_m is the number of mineral phases. The reaction rate $Rim$ is commonly a function of solution composition including pH, but also the exposure of the mineral to the pore water, defining its intrinsic reactivity. Reactivity is commonly controlled by the reactive surface area. A mineral dissolution precipitation rate can therefore be expressed as:

where S_i [m² mineral dm⁻³ bulk] is the reactive surface area for the mineral i, k_i is the surface area normalized rate constant [mol m⁻² mineral s⁻¹], C^l_j is the concentration of a chemical species in solution [mol L⁻¹ H₂O] and N_s is the number of dissolved species.

Equations (1–6) include several parameters that can evolve in response to biogeochemical reactions, namely:

• porosity ϕ permeability tensor k;

• tortuosity τ_α;

• mineral reactive surface area S_i as a surrogate for reactivity.

Figure 1 depicts the dependencies of groundwater flow, transport processes and geochemical reactions on these parameters and illustrates how changes of these parameters lead to interdependencies and feedback mechanisms between the flow, transport and reaction domains. Since the dissolution of mineral phases will lead to a decline in the volume filled with the solid phase, an increase in porosity will result. The increase in porosity, in turn, will tend to cause changes in permeability and tortuosity. In addition, the reactive surface area may also be altered as a mineral progressively dissolves. Based on this analysis and considering Equations (1–6) together, it becomes clear that mineral dissolution has the potential to directly affect groundwater flow (both porosity ϕ and k are featured in Eqn. 1 and 2), advective solute/gas transport, since the specific discharges in Equation (3) depends on k, and diffusive solute and gas transport, since the diffusion coefficients are a function of τ_α (Eqn. 4). In addition, the fact that porosity is included in some terms of the multicomponent reactive transport Equation (4), but not in others, reveals the effect of changes in the water-rock ratios on the mass balance. Similar considerations are valid for precipitation reactions or the formation of biofilms. It becomes evident that biogeochemical processes have the potential to cause several feedback mechanisms, increasing the inter-dependencies of the equations governing flow and reactive transport (Fig. 1). In order to understand the evolution of these parameters at the continuum scale, it is instructive to review how biogeochemical processes control porous media evolution on the pore scale.

As most of the relevant physicochemical processes occur at the pore-scale, an understanding of processes that occur at this scale is essential for developing an appreciation on how flow, transport and reaction parameters evolve, how these processes and parameters are linked to each other, and how these parameters affect reactive transport both at the pore and continuum scales.

With the advent of modern, non-invasive imaging techniques, it has become possible to visualize the structure of porous media on the pore scale. Werth et al. (2010) provide a review of these methods including X-Ray microtomography, Nuclear Magnetic Resonance Imaging (NMR) and optical imaging methods, as well as their applications in hydrogeology and reactive transport. Specifically, the use of 3D X-Ray microtomography has found increasing applications due to its high resolution. These imaging techniques allow visualizing pore structure at sub-micron resolution, providing a window into the complex nature of porous media.

Figure 2a depicts a 3D micro-CT image of the pore space of a limestone (Blunt et al. 2013). This limestone sample is characterized by large well-connected pores (Blunt et al. 2013). Based on such images, it is possible to construct pore network models (Fig. 2b), in which the pores are represented as spheres, while the pore throats are shown as cylinders, connecting the pores with each other (Blunt et al. 2013). In many situations, a fraction of the pores can be isolated, i.e., disconnected (trapped porosity), while other pores are difficult to reach (Navarre-Sitchler et al. 2009). Although connected to the pore-network, they are “dead-end” pores, and therefore do not contribute substantially to the “effective porosity” (Navarre-Sitchler et al. 2009). The connectivity, in principle, is defined by the ratio between the connected (i.e., effective) porosity and the overall porosity.

Examining the pore network shown in Figure 2 clarifies that flow and transport through such a material will heavily depend on the size of the pores, the diameter of the pore throats, and the connectivity between the pores. Large pores and direct connectivity in three dimensions implies a low degree of tortuosity, facilitating flow, and both diffusive as well as advective transport with relative ease.

In the case of an evolving porous medium affected by mineral dissolution–precipitation reactions, it is evident that both porosity and connectivity can be affected. The dissolution of a solid phase is directly linked to the increase of the overall porosity. However, dissolution does not necessarily increase effective porosity to the same degree. The effect can be substantial, if dissolution widens and opens up pore throats or opens up new diffusion pathways (Knackstedt and Zhang 1994; Moreira and Coury 2004), or it can be more limited, if dissolution is more equally distributed (Algive et al. 2010). Similarly, mineral precipitation might have a limited effect on connectivity, if minerals form on existing surfaces without clogging pore throats. On the other hand, precipitation may lead to the disconnection of existing solute migration pathways (Algive et al. 2010), with a pronounced impact on flow and transport parameters. These theoretical considerations imply that for the same change in total porosity, the dissolution and precipitation of minerals can lead to very different impacts on effective porosity and connectivity. As a consequence, flow and transport will also be affected in different ways due to the evolution of tortuosity and permeability as a function of effective porosity and connectivity, but not simply total porosity (Navarre-Sitchler et al. 2009).

In addition to changes in transport properties, the dissolution and precipitation of mineral phases, can result in changes in reactivity, leading to enhancement or reduction of reaction rates Normally, the dissolution of a mineral phase is controlled by mineral-specific surface-area-normalized reaction rates and by the reactive surface area of the mineral (Lasaga 1981, 2014).

The reactive surface area is controlled by a variety of factors, such as grain size (Emmanuel and Berkowitz 2007; Liu et al. 2008), exposure of the mineral and contact with the pore water (Landrot et al. 2012), the degree of occlusion (Peters 2009; Lai and Krevor 2014), surface roughness and presence of imperfections or etch-pits (Deng et al. 2018). Moreover, when multiple minerals are present, the surface fraction occupied by a certain mineral might significantly differ from its volume fraction (Lai and Krevor 2014). Also, the flow conditions can lead to more reactive zones, contributing to commonly observed deviations of reactive from geometric surface area. Together, these factors control the magnitude of the effective reactive surface area (Jung and Navarre-Sitchler 2018).

Any of these factors can change over time, as a porous medium progressively evolves, simply through shrinking of grain size, the enhancement of surface roughness and etch pit formation during dissolution (Teng 2004; Buss et al. 2007), incongruent dissolution as shown in Figure 3a for pyrite oxidation, or the formation of precipitates on the surface of a dissolving mineral. Secondary Fe(III)-oxy-hydroxides tend to precipitate on the surface of pyrite grains (Fig. 3a). Similarly, Figure 3b shows the precipitation of secondary oxide and clay minerals (goethite and vermiculite) on the surface of biotite, which is undergoing weathering. Incongruent dissolution and the formation of surface precipitates limit the interaction of the bulk solution with the surface of the minerals and in this way, tends to limit further mineral dissolution (Daval et al. 2009).

Considering that mineral-dissolution precipitation reactions have a pronounced impact on porous media evolution, we review these processes in more detail in the following sections, making a distinction between dissolution and precipitation reactions. In this context, we emphasize the evolution of porosity and transport parameters, as well as the interactions between the rate of flow and transport processes and reactions and the resulting effect on porous media evolution. Processes are discussed predominantly at the pore scale, keeping implications for the continuum scale in mind.

In many situations, sedimentary deposits or fractured rocks are affected by weathering reactions involving the dissolution and alteration of their solid structure (Bjørlykke et al. 1989; Aharonov et al. 2004; Navarre-Sitchler et al. 2011; Jin et al. 2013). The oil and gas industry performed a substantial amount of research dedicated to understand and optimize the process of permeability enhancement via mineral dissolution with the goal of increasing production rates (Wang et al. 1993; Fredd and Fogler 1999; McDuff et al. 2010). Another example with relevance to the mining industry is the enhancement of in-situ uranium recovery (Simon et al. 2014; Regnault et al. 2015; Lagneau 2019, this volume) by using acidic solutions to dissolve uranium-containing minerals. In the case of construction materials, the dissolution of primary minerals is commonly a process that is unwanted. Indeed, the increase of porosity usually induces a decrease of the material's intended performance. Negative effects of leaching and acid attack on the properties of cement and concr have been demonstrated (Atkinson et al. 1987; Bentz and Garboczi 1992; Bertron et al. 2005).

As dissolution reactions can impact groundwater flow and solute transport through porous media, reactive transport modelling of the aforementioned problems has to account for the impact of dissolution on the relevant material properties. To predict how a porous medium's properties will evolve through dissolution reactions, it is necessary to understand the processes that control the development of dissolution patterns. Previous investigations of mineral dissolution have revealed that, under certain conditions, stable dissolution fronts develop (Barlet-Gouedard et al. 2006; Hesse et al. 2013; Walsh et al. 2013), while in other cases, such as karst development, dissolution fronts become unstable, leading to the development of fingering or wormholing (Hoefner and Fogler 1988; Daccord et al. 1989; Steefel and Lasaga 1990; Golfier et al. 2002).

Aharonov et al. 2004; Algive et al. 2012 and Luhmann et al. 2014 have shown that the dissolution regime depends on the characteristic times for advection, diffusion and reaction. Typically, relationships between characteristic times of transport and reaction processes are defined by the Peclet (Pe) and Damköhler (Da) numbers. The Peclet number addresses the relative importance of the advective and diffusive transport (Algive et al. 2010), while the Damköhler numbers reflect the ratio between the characteristic times of reaction compared to that of transport (Tartakovsky et al. 2007; Min et al. 2016). The Peclet number is defined as:

where v represents the groundwater flow velocity, D the diffusion coefficient and L is a characteristic length of interest in the porous media. A high Peclet number reflects an advection-dominated transport regime. Two versions of the Damköhler number exist, which here will be referred to as the advective and diffusive Damköhler numbers, respectively. The Da numbers provide a relationship between the reaction rate and the transport (advection–diffusion) rates

Experimental and numerical research has focused on the prediction of wormhole formation in carbonate rocks (Fredd et al. 2000; Buijse et al. 2005; Quintard et al. 2007; Budek and Szymczak 2012).

Fredd et al. (2000) reviewed the response of carbonate rocks to ingress of acidic waters. Figure 4 depicts the different dissolution patterns for different acid injection rates (increasing from left to right). At very low injection rates (low Pe, high Da_adv), uniform face dissolution occurs near the injection point where all the acid is consumed. Increasing Pe facilitates deeper ingress of acidic waters along the main flow pathway where reactions take place, leading to a conical or dominant wormhole. The distinction between the conical/dominant wormholes and ramified or uniform dissolution was further evaluated by Detwiler et al. (2003). When the injection rate is further increased (to the point where Da_adv becomes ≪ 1), acidic waters are forced through the entire porous medium without having time to react, leading to the formation of a ramified wormhole and ultimately uniform dissolution. The effects of Peclet and Damköhler number were analyzed with numerical modeling as early as 1990 (Steefel and Lasaga 1994), who showed the transition from dominant wormholes to ramified wormholes and uniform dissolution with decreasing Damköhler number.

These diverse dissolution patterns lead to different evolution of transport properties. As uniform facial dissolution leads to a focused increased in permeability near the injection point, the macroscopic transport properties of the entire material do not change significantly. On the other hand, uniform dissolution tends to increase permeability throughout the entire sample. For situations between the two extremes cases, the development of a broad and relatively linear transport pathway ensues, through which most solutes can be transported. The latter case results in the formation of wormholes.

Min et al. (2016) performed a pore-scale study represented in Figure 5, which emphasizes the role of Da_diff. Figure 5a depicts a reaction-limited regime (Da_diff ≪ 1) subject to diffusion-dominated transport (Pe ≪ 1). In this case, dissolution is relatively slow, leading to solute release throughout the domain. However, since diffusion is the dominant transport process, the embedded solid inclusions dissolve rather uniformly. Figure 5b considers a diffusion-limited regime (Da_diff > 1) with low flow (Pe ≪1). Reaction rates are faster than diffusion rates and one can see that dissolution reactions occur locally, leading to a steep concentration front. In Figure 5c, (Pe and Da_diff > 1), front instabilities are developed due to the advection dominated transport regime. Min et al. (2016) also simulated the evolution of permeability. The case depicted in Figure 5a corresponds to the largest increase in permeability, while the scenario shown in Figure 5b leads to the lowest increase in permeability.

The complexity of the dissolution regimes reported by Fredd et al. (2000) and Min et al. (2016) demonstrates that a quantitative assessment of dissolution patterns and rates, and feedback on porosity and permeability is difficult solely based on theoretical considerations. To properly take pore scale processes into account, pore-scale simulations of processes for relevant pore geometries provides a sensible path forward for deriving upscaled relationships for continuum scale modeling of evolving porous media. Recent developments in imaging and pore-scale modeling methods have made it possible to investigate these phenomena in more depth (Pereira Nunes et al. 2016; Soulaine and Tchelepi 2016; Tian and Wang 2018). Gao et al. (2017) and Gray et al. (2018) were able to approximately reproduce the experimentally observed dissolution patterns in microstructures obtained from CT scan images. Pre-and post-experimental X-ray synchrotron imaging of reactions in a fracture developed in Duperow Dolomite quantified the wormholing taking place due to the influx of acidic fluid, and this effect was captured in both reduced dimension models (Deng et al. 2016) and in pore scale models (Steefel et al. 2013).

The precipitation of secondary minerals often impacts groundwater flow and solute transport. In some cases, precipitation reactions lead to the self-sealing of fractures (Steefel and Lasaga 1994; Dobson et al. 2003b; MacQuarrie and Mayer, 2005; Gherardi et al. 2007) or the reinforcement of the solid structure (Bickmore et al. 2001; Pfingsten 2002; Jacquemet et al. 2012). Similarly, carbonation of cementitious materials under water can lead to a calcite crust forming on the surface of the material (Dauzères et al. 2010; Seigneur et al. 2017), making the material almost impermeable. Crust formation has also been observed as a result of gypsum precipitation on historical monuments (Dewanckele et al. 2012). Furthermore, permeability decrease in sedimentary basins (Aharonov et al. 2004) or soils (White et al. 2005) due to mineral precipitation have been observed. The formation of secondary minerals may lead to fractures through generated internal mechanical stresses. This can lead to the deterioration of mechanical properties (Chagneau et al. 2015; Li et al. 2017). For cement materials, the delayed formation of ettringite (Taylor et al. 2001) or atmospheric carbonation (Auroy et al. 2013) are typical examples of precipitation processes that can induce fracturing, which can be detrimental in the context of CO₂ sequestration or radioactive waste management. In all the aforementioned cases, precipitation can lead to clogging, which significantly affects flow and solute transport. It is therefore important to understand the role of precipitation reactions on transport properties and to consider these evolving parameters in modeling analyses.

Kang et al. (2005) performed numerical studies of precipitation reactions involved in CO₂ sequestration. Working with 2D simplified images, they studied the influence of Pe and the Da_diff number on mineral precipitation reactions and subsequent permeability reduction. For a fixed high Pe, these simulations indicate that Da_diff controls the time required to plug the medium. Larger Da_diff lead to more localized precipitation near the inlet and more efficient clogging, reducing permeability to zero. Conversely, increasing Pe (for a fixed Da_diff) has the opposite effect. The efficient clogging of the medium relies on moderate Da_diff numbers and small Pe numbers. A small amount of precipitates can render the pore structure impermeable. For higher Pe, it takes more time (and more reactants) to clog the porous medium, as the precipitation front is wider.

Tartakovsky et al. (2007) performed a similar study and obtained comparable results, which are depicted in Figure 6. A high Pe number combined with a low Da_diff number resulted in precipitation throughout the porous medium (Fig. 6a). For high Pe and increased Da_diff, precipitation mostly occurred in the main flow pathways (“dendrite-like precipitation” depicted in Fig. 6b). In the case of low Pe and high Da_diff, most of the precipitation took place near the inlet, rapidly clogging the medium, leading to a rapid and sharp decrease of flowrates, as depicted in Figure 6c. Tartakovsky et al. (2008) obtained similar results for mixing-induced precipitation involving different solutes. For high Da_diff, the precipitation zone is thin, restricted to the region where both reactants meet. For lower Da_diff, as diffusion is more dominant, the precipitation zone is widened.

These examples emphasize the complex interactions between transport and reaction processes and relations between their time scales. Even if the total amount of precipitates is correctly captured, it is not straightforward to predict the evolution of transport parameters, which depend strongly on the distribution of the precipitates and not only their amount, especially in clogging situations.

While wormholing dissolution patterns tend to exert a positive feedback loop in terms of permeability enhancement, precipitation reactions can have the opposite effect by preferentially occurring in the high-flow zones (Steefel and Lasaga 1994; Aharonov et al. 1998; Mehmani et al. 2012), causing flow to diverge towards zones of lower permeability. This behavior can be described as a “dewormholing” process.

However, other factors can also affect the location and distribution of precipitate formation. Rajyaguru et al. (2019) studied through-diffusion in chalk samples. Precipitation of barite and gypsum were investigated under a purely diffusive regime. Gypsum is characterized by a rather high solubility and kinetic rate compared to barite. It was observed through micro-CT analysis that the precipitation patterns were significantly different for barite and gypsum. While barite precipitated as a thin continuous disk at the center of the chalk sample (even though its intrinsic rate is lower), gypsum mostly precipitated as large isolated spheres. Even though the overall porosity decrease was similar (~2%), the diffusive properties were more substantially impacted in the barite experiment. The authors suggest that the barite precipitation was the result of both homogeneous and heterogeneous nucleation. Conversely, the heterogeneous precipitation of gypsum seemed to have been highly influenced by local pore structure heterogeneities. These experiments reveal the complexity of mineral precipitation within pore structures, even in simple diffusion-controlled transport regimes. Chagneau et al. (2015) showed that precipitation of celestite in a diffusion-reaction cell occurred almost exclusively in macropores based on the differential transport of anions (excluded from the smallest pores due to electrostatic effects) and cations.

Other studies also showed the importance of spatial heterogeneities for determining the distribution of precipitation, which in turn will control the evolution of tortuosity, permeability and reactivity of porous media (Andreani et al. 2009; Fox et al. 2016; Garcia-Rios et al. 2017).

In many situations, mineral precipitation reactions occur in tandem or in response to dissolution reactions, leading to feedback between these reactions. Singurindy and Berkowitz (2003) provided images of the coupled dissolution/precipitation of carbonate and gypsum and showed how the dissolution patterns influenced the precipitation patterns. The evolution of a porous medium at the grain scale subject to dissolution and precipitation reactions is conceptually represented in two dimensions in Figure 7.

The initial pore structure depicted in Figure 7 (top left quadrant) exhibits several features. The region surrounding point (a) is critical, as it is well connected, providing a main pathway for water flow and solute migration. The region around point (b) represents a rather tortuous pathway for diffusion and advection. The region around (c) on the other hand, consists of large pores, which are poorly connected to the main transport pathways. Figure 7 also emphasizes that modifications to a pore structure can have varied impacts on macroscopic properties.

In a dissolution-dominated regime (top right quadrant of Fig. 7), dissolution of primary minerals can lead to the creation of new flow and transport pathways (1) or increase the size of critical pore throats (2). This evolution can have significant impacts on flow and transport processes. In addition, dissolution will impact the reactive surface of the initial solid structure (3) or can lead to the formation of leached layers in the case of incongruent dissolution (6).

In a precipitation-dominated regime (bottom left quadrant of Fig. 7), precipitation of secondary minerals can have opposite effects: decline of connectivity, leading to a reduction of flow and transport pathways (1), decrease of critical pore size (2), trapping of initially connected porosity (4). In addition, precipitation reactions can limit access to (4) or passivate (5) a fraction of the initial reactive surface.

In a regime that is both affected by dissolution and precipitation (depicted in the lower right quadrant of Fig. 7), coupled dissolution–precipitation reaction can lead to more complex effects. Figure 7 clarifies that the modification of the pore structure depends on the respective dissolution and precipitation kinetics and interactions between dissolution and precipitation reactions.

Despite obvious challenges with the parameterization of evolving porous media, simulation tools are needed at the continuum scale to assess a variety of relevant environment and engineering problems. In the following, we provide a review of the most common approaches found in the literature to parameterize evolving parameters for reactive transport simulations at the continuum scale.

As discussed above, the dissolution and precipitation of minerals lead to changes in porosity. From a modeling perspective, it is fairly straightforward to account for total porosity changes at the continuum scale. At any point in time, the porosity can be calculated based on the mineral volume fractions that constitute the grain matrix. In this context, porosity is defined as:

It must be kept in mind that porosity defines total porosity and does not take into consideration the fraction of disconnected porosity, dead-end pores and effective porosity. In many situations, the mineralogical composition of porous media is not well known and/or only a few mineral phases contribute significantly to porosity changes. Under such conditions, it is more practical to define porosity at a specific location in space and point in time as:

where φ_NR is the volume fraction occupied by non-reactive minerals and N_{m, r} is the number of reactive minerals. Based on these considerations, the change in porosity can then be estimated as follows:

Similar approaches can be taken to update effective porosity due to bioclogging, i.e., the formation of biofilms. While the previous equation accurately describes the physical processes at play, it does not provide any information about the structural evolution of the pore structure. Looking back at the few referenced studies of the previous section, this obviously constitutes a shortcoming. However, as we will discuss in the following sections, it is possible to attribute different continuum relationships for the porous media properties, in order to represent indirectly the pore structure evolution. This approach has been commonly used in reactive transport modeling and the evolving porosity has been used to update transport parameters such as tortuosity and permeability.

In the absence of significant advective fluxes (e.g.: in low permeability materials, low hydraulic head gradients), diffusion is the main driving process for solute transport. Diffusion-dominated transport is relevant for deep geologic repositories for long-term storage of radioactive waste. In this context, reactive transport modelling has been used to develop an improved understanding of the degradation of cementitious materials at cement/clay interfaces (Gaucher and Blanc 2006; Savage 2013; Jenni et al. 2017; Dauzères et al. 2019). Since these interfaces are affected by dissolution–precipitation reactions that induce porosity changes, diffusion properties of the materials are also affected, directly through changes of porosity, as well as through the evolution of tortuosity. Despite known limitations, it is therefore necessary to provide relationships able to describe the evolution of tortuosity as a function of changes in porosity.

A review of the literature reveals that many tortuosity–porosity relationships exist to describe the dependence of diffusion on porosity. These relationships have been developed to estimate effective diffusion coefficients based on porosity, a parameter that is easier to measure than the effective diffusion coefficient. These relationships are also commonly used in evolving porous media. Akanni et al. 1987; Boudreau 1996 and Ray et al. 2018 present a summary of the commonly used relationships. However, most reactive transport simulations in saturated porous media rely on Archie's law to account for the porosity dependence of the pore diffusion coefficient (Xie et al. 2015), which, in its simplest form can be written as:

with n being the cementation factor. Let us recall that, in this study, the effective diffusion coefficient for a fully saturated aqueous phase is defined as:

Some confusion usually exists about these different relations. In some cases, studies relate a relationship between porosity and effective diffusion coefficient, instead of tortuosity. Other forms of Archie's law can be found, some of them describing a percolation threshold under which diffusion cannot occur (e.g., Bentz and Garboczi 1991 for cement materials). Another example is the relation developed for the effective diffusion coefficient in the liquid phase by Navarre-Sitchler et al. (2009), which will be further studied in a case study :

where $Dlmin⁡$ represents the effective diffusion coefficient of the matrix, and where the effective porosity, ϕ_e, is defined as:

Another commonly used tortuosity–porosity relationship assumes an exponential relationship to porosity:

where a is an experimental fitting factor. Let us emphasize that these laws were usually developed to link the effective diffusion coefficient to porosity.

Under unsaturated conditions, it has been shown that tortuosity also depends on phase saturation. A formulation that is commonly used to estimate porosity was derived by Millington and Quirk (1961):

This relationship, which shows a strong and non-linear dependency of tortuosity on phase saturation S_α, was derived for the gas diffusion in porous media (Millington 1959). However, most reactive transport codes adopt the Millington–Quirk relation (Millington and Quirk 1961) to model diffusion through the liquid phase as well. Considering the opposite wetting behavior of the gas and liquid phase, the use of similar relations to describe tortuosity in these different phases remains questionable. Chou et al. 2012 provides a review of different tortuosity models for unsaturated conditions and concludes that the use of the Millington–Quirk relation overestimates experimental data.

Figure 8 depicts tortuosity as a function of porosity based on common relationships for saturated media, some of which are documented in the literature (Maxwell 1881; Petersen 1958; Tomadakis and Sotirchos 1993), while others are variations of Archie's law (Bentz and Garboczi 1991) or represent an exponential dependence on porosity (Mainguy et al. 2000; Seigneur et al. 2017). It can be observed that tortuosities obtained from the various relationships differ substantially, in particular for low porosities. The range between minimum and maximum tortuosities are approximately 2 orders of magnitude for a porosity of 0.5, while calculated tortuosities differ by more than 5 orders of magnitude for porosities less than 0.1. These pronounced differences suggest that tortuosity–porosity relationships are specific to porous media. The fact that some relationships do not yield a tortuosity of “zero” for a porosity of “zero” (even though the effective diffusion coefficient does go to 0) implies that these relationships are not applicable for conditions of complete pore clogging. However, there is experimental evidence that in some cases mineral precipitation will not completely clog the porous media. Li and co-workers observed that despite the formation of a distinct calcite cement band in cement subject to CO₂ attack, the porosity did not go completely to zero as evidenced by the continued progress of the primary portlandite dissolution front (Li et al. 2017). They attributed this to be due to either the lack of precipitation in nanopores that were still able to allow for diffusive transport of the uncharged H₂CO₃ molecule into the unaltered cement, or alternatively due to continued micro-fracturing. The lack of clogging of the nanopores is compatible with the interpretation of anion (and cation) transport in clay-rich material in which anion exclusion occurs (Chagneau et al. 2015).

Leaching of cementitious materials has been studied extensively, since dissolution reactions lead to increased diffusion into these materials, enhancing their deterioration. The prediction of degradation depths of these materials strongly depends on which cementation (or exponential) factor (Shao et al. 2013) or which feedback law is considered (Seigneur et al. 2017; Georget et al. 2018). This is due to the fact that dissolution rates of primary minerals are affected by the diffusion of calcium through the degraded layers. Several reactive transport studies managed to adequately reproduce the extent of measured degradation depths using either Archie's law (Galíndez et al. 2006; De Windt and Badreddine, 2007; Chagneau et al. 2015; Li et al. 2017; Georget et al. 2018) or an exponential relation (Mainguy et al. 2000; Seigneur et al. 2017) to describe tortuosity evolution as a function of changing porosity. However, the choice of the cementation factor n or the exponential factor a differs between the studies in a way that is not straightforward to predict, depending on the conditions and the materials considered.

Most of these modelling studies are based on an initial porosity and tortuosity representative for the investigated materials, and only use these empirical laws to compute the evolution of diffusion parameters. In these studies, the initial diffusive properties of the material are either known independently or were fitted to early time data. Since dissolution under a diffusion-controlled regime occurs in a fairly uniform manner and leads to a relatively mild increase in tortuosity, this approach was able to produce reasonable matches with observational data, albeit with different fitting factors.

However, when precipitation reactions occur simultaneously, it becomes more complicated to assess the relative dynamics between dissolution and precipitation reactions and their effect on tortuosity and diffusion properties. Several studies had to artificially increase the empirical factors when describing precipitation reactions (Huet et al. 2010; Brunet et al. 2013; Chagneau et al. 2015; Lalan 2016) or the reactivity (Jacquemet et al. 2012), or use other empirical assumptions (Walsh et al. 2013) to be able to reproduce breakthrough curves and/or degradation depths when simulating the carbonation of these materials.

These results are consistent with pore scale analyses (see above), implying that limited and localized mineral precipitation can significantly alter diffusive migration paths. Hence, the impact of precipitation reactions on diffusion is more difficult to assess than that of dissolution. These observations also imply that different materials are affected in distinctive ways, which is inherently related to the initial organization of their pore structure. Materials with a high degree of initial tortuosity, containing small pore throats, will be more sensitive to the effects of precipitation reactions, which was observed both experimentally and in numerical modelling studies (Kutchko et al. 2007; Brunet et al. 2013; Huber et al. 2014; Seigneur et al. 2017).

These results suggest that it might be necessary to use different tortuosity–porosity relationships for materials affected by dissolution or precipitation reactions. However, since these reactions are most commonly coupled and occur simultaneously, a proper assessment of the effects of these reactions on diffusivity remains a challenging issue.

Moreover, previous studies have shown that obtaining a better representation for the evolution of diffusion (sometimes through the use of unrealistic high cementation factors) can make it difficult to reproduce observed mineral precipitation and dissolution patterns. Jacquemet et al. 2012 indicated that the evolution of reactivity also plays a role in the degradation depth, although most of the research to date has focused on the evolution of the transport properties.

Permeability and its evolution plays a significant role in advection-driven transport regimes, including the injection of fluids in disequilibrium with the surrounding rock formation (acidizing carbonate rocks or uranium recovery via sulfuric acid leaching). In these regimes, the modification of porosity can also impact permeability, and therefore groundwater flow patterns, and in turn the degree of chemical perturbations which alter the porous medium. It is therefore necessary to utilize relationships capable of capturing the impact of evolving porosity on permeability.

Like relationships between tortuosity and porosity, permeability–porosity relationships have been developed independently from the field of reactive transport modeling. The Kozeny–Carman relationship is the most commonly used permeability–porosity relationship in reactive transport simulations (Poonoosamy et al. 2018) for simulating the evolution of permeability affected by mineral dissolution–precipitation reactions. It is commonly used in the following form:

where k₀ and ϕ₀ respectively represent the initial or reference permeability and porosity. Xu and Yu (2008) and Hommel et al. (2018) provide an extensive review of the most common permeability–porosity relationships. In general, power-laws as a function of porosity have been widely used to describe the evolution of permeability:

As an alternative, Verma and Pruess (1988) have introduced a critical porosity ϕ_crit in the power law:

which corresponds to the porosity for which permeability becomes zero. The parameter n represents an empirical fitting parameter, and could be considered as similar in concept to thresholds for diffusion found in modifications of Archie's Law (Bentz and Garboczi 1992; Navarre-Sitchler et al. 2009). Values between 3 and 75 have been reported for the fitting parameter n (Brunet et al. 2013). Gouze and Luquot (2011) suggested that this wide range of power values can be explained by the fact that permeability is controlled by two parameters: the first parameter corresponds to the tortuosity of the pore space. The second parameter corresponds to the pore throat size, or the “hydraulic radius” (Knackstedt and Zhang 1994; Moreira and Coury 2004) representing the characteristic radius of the pore channels. Mercury intrusion porosimetry and the “ink-bottle effect” (Diamond 2000) provide useful experimental illustrations of the pore-throat size control on permeability. During the first stages of wormhole development, an increase of hydraulic radii occurs, while when channeling develops, tortuosity (as defined in this chapter) increases significantly. In more general terms, Hao et al. (2013) and Menke et al. (2016) suggest that the power value is representative of the heterogeneity of the pore structure.

Figure 9a depicts the common permeability–porosity relationships for an arbitrary reference porosity of 0.3. Similar to the dependence of tortuosity on porosity, it can be observed that the various relationships lead to widely differing results. A porosity increase, representing the case of mineral dissolution, leads to differences on the order of one order of magnitude between the relationships for a porosity enhancement of 0.3. On the other hand, much more significant differences are seen for the case of precipitation and porosity reduction. It is worth pointing out that in certain cases, inverse correlations were found between permeability and porosity (Luquot et al. 2013, 2016; Garing et al. 2015). Such correlations usually arise when coupled dissolution and precipitation processes are occurring at the same time, but do not involve the same pore sizes (Luquot et al. 2014; Garing et al. 2015). None of the permeability–porosity relationships introduced above could represent these inverse relationships.

The evolution of permeability is more difficult to capture than the evolution of diffusivity. Indeed, in diffusion-dominated regimes, reaction fronts tend to be rather stable. Conversely, in advection-dominated regimes, pore-scale heterogeneities affect the distribution of flowrates and can induce unstable fronts and channeling.

Power laws are often used to describe permeability–porosity relationships. These relationships strongly depend on the initial material's pore structure and its heterogeneity (Rötting et al. 2015; Menke et al. 2016; Wolterbeek and Raoof 2018). Civan (2001) uses three different power laws to describe Fontainebleau sandstones that are similar in nature but have different initial porosities.

It has also been shown that the fitting parameter n depends on the Pe and Da numbers, since they influence the dissolution and precipitation patterns (Fig. 9b). Egermann et al. (2010) showed that the permeability evolution in five different dissolution experiments was well represented by power laws; however, each required a different fitting parameter. The values of the power increased from 0.5 (relatively uniform dissolution of grains) to 31 (dissolution focused on pore throats). Colón et al. (2004); Panga et al. (2005); Menke et al. (2016) provide additional examples with similar outcomes.

In addition, Smith et al. (2013) observed that the value of the fitting parameter n evolved with time as dissolution progressed. This was likely due to variations in heterogeneity caused by dissolution reactions in advective regimes. In other words, a constant fitting parameter n might not be suitable to describe the evolution for substantial modifications of porosity (Pape et al. 1998; Bernabé et al. 2003; Nogues et al. 2013; Ovaysi and Piri 2014; Vialle et al. 2014).

Considering that mineral dissolution reactions can lead to wormhole formations and precipitation reactions can lead to the closure of wormholes and fractures, the form of permeability–porosity relationships are typically different for precipitation and dissolution reactions (Bernabé et al. 2003; Li et al. 2010), as shown in Figure 9b. Furthermore, in cases of coupled dissolution and precipitation or dissolution–precipitation cycles, hysteresis was observed in the permeability–porosity curves (Algive et al. 2012; De Boever et al. 2012; Van der Land et al. 2013), whose importance depends on the Pe and Da numbers (Algive et al. 2012). In general, hysteresis is limited in the case of diffusion-dominated regimes with low kinetic rates (low Pe and high Da_diff), consistent with expected behavior, as discussed above in the section “Effect of mineral dissolution and precipitation on porous media” (Tartakovsky et al. 2008; Min et al. 2016).

These considerations explain why the permeability–porosity relationships have limited predictive capabilities in the context of evolving porous media. Indeed, for identical pore structure and otherwise identical conditions, considering diverse minerals with different kinetic rates (i.e., different Da numbers), experimental data can only be reproduced with different fitting parameters n. It seems unlikely that a single model, however complex, could predict changes in permeability in evolving porous media under different flow and transport conditions. These observations emphasize the need to choose an appropriate power-value, based on a priori knowledge or constrained by observational data.

In fractures, it is customary to use the cubic law for permeability, which is related to the local fracture aperture (and thus, local fracture porosity). This relation has a long history in the literature, with the most cited in the geosciences being perhaps the study by Witherspoon et al. (1980). Generalization to a local cubic law for rough fractures was discussed by Oron and Berkowitz (1998). The interested reader is referred to the chapter by Deng and Spycher (2019, this volume) for further discussion.

In addition to the evolution of transport parameters, it is equally important to describe the evolution of reaction rates, i.e., the reactivity in evolving porous media. An adequate description of rates is necessary to capture long-term solute concentrations and loadings, pH-buffer reactions, and to adequately account for feedback with flow and transport processes.

Reactivity is often described via a reactive surface area of the minerals. Several models are commonly used to describe the evolution of the surface area S. The simplest formulation considers the uniform dissolution of spherical grains, in which case the surface area is inversely proportional to its radius:

Based on this dependence, a two-third power law relationship can be derived to describe the evolution of surface area with respect to a mineral's volume fraction (Steefel and Lichtner 1998):

where S₀ is a reference surface area and φ₀ defines the corresponding reference mineral volume fraction. This relationship is commonly used to describe surface controlled mineral dissolution reactions, accounting for a decrease in rate as dissolution proceeds (Fig. 10). Although this relationship is derived based on the dissolution of a single grain, it has been used frequently on the continuum scale in the form of Equation (23) (Steefel and Lichtner 1998). A decline in reactivity and surface area is intuitive; however, it does not hold under all circumstances. Dissolution reactions can also lead to an increase of reactive surface area, as illustrated by the sugar-lump model developed by Noiriel et al. (2009). Similarly, the formation of etch pits on the mineral surface may also lead to an increase in reactive surface area, although the reactive surface area (like the mineral volume fraction) must ultimately go to zero with continued dissolution.

Luquot and Gouze (2009) proposed a power-law relationship to describe the evolution of surface area as a function of porosity:

This approach allowed Luhmann et al. (2014) to fit the decrease of reactivity as evidenced by decreasing concentrations in the breakthrough curves of their experiments.

However, surface reactions are also affected by pore scale flow dynamics. Li et al. (2008) considered the role of a diffusion boundary layer in the dissolution of calcite within a single cylindrical pore. Noiriel et al. (2012) also considered the role of a diffusion-controlled boundary layer on mineral dissolution rates, an effect that is expected in all but the most vigorously stirred system. Molins et al. (2014) performed pore-scale simulation of the injection of CO₂ into crushed calcite. While only a small amount of calcite dissolved in these simulations (porosity variation negligible), the results showed that mass transport limitation to the reactive surface also result in an important control on mineral dissolution. These conclusions were reinforced by a recent pore scale study on flow and transport in discrete reactive fractures (Molins et al. 2019). The shrinking sphere model depicted in Figure 10 assumes that concentrations in the bulk solution and in water directly in contact with the mineral surface are fully mixed and identical, which constitutes a simplification of the processes at work.

Minerals often do not undergo complete dissolution, but dissolve incongruently. Such dissolution reactions have been described by the shrinking core model (Levenspiel 1998), as depicted in Figure 11, in particular with applications to sulfide mineral oxidation (Ritchie 1994; Wunderly et al. 1996; Lefebvre et al. 2001; Mayer et al. 2002). In the case of sulfide oxidation, the rate-determining step is the diffusion of dissolved oxygen towards the mineral surface through a leached (oxidized) layer. As sulfide oxidation progresses, the thickness of the leached layer grows, and continued oxidative dissolution slows substantially, leading to a negative feedback on the grain scale. This model is not restricted to sulfide mineral oxidation, it can be applied to different minerals, also subject to incongruent dissolution and leached layer formation, or composite mineral grains, representing conditions when reactants diffuse through the porous leached layer generated by the dissolution of one of the minerals present in the grain (Steefel and Lichtner 1994, 1998; Deng et al. 2016). The shrinking core model can also be used to represent the effect of the formation of surface coatings associated with the precipitation of secondary minerals. In general, these considerations emphasize the mass transport limitations can affect reactivity and can significantly decrease the reaction rates (Dentz et al. 2011).

Assuming spherical particles, an effective first order rate constant capturing the effect of the formation of a leached layer or surface coating can be described as:

where $Die,m⁡$ represents the effective diffusion coefficient of the species driving mineral dissolution through the protective surface layer, $rip⁡$ is the initial radius of the particle and $rir⁡$ is the radius of the unreacted portion of the particle. Although the shrinking core model is derived for a single particle, this formulation makes it possible to extend its application to the continuum scale.

The application of the shrinking core model was shown to capture experimental observations when applied to the weathering of mine wastes (Wunderly et al. 1996; Brookfield et al. 2006; Demers et al. 2013; Soleimani et al. 2013; Doulati Ardejani et al. 2014; Smith et al. 2018), leaching of ore minerals (Ferrier et al. 2016), and for thermochemical sulfate reduction (Bildstein et al. 2001).

Although practical, the shrinking core model has several drawbacks: its formulation commonly does not account for a local mass balance during incongruent dissolution and is normally not accounting for the elements that are left behind in the leached layer and the elements that are released into the bulk solution. In addition, the formulation requires assumption of an initial leached layer thickness, which further adds to the empiricism of the model. An alternative approach would be to develop a formulation that is able to account for transition from surface-controlled dissolution to diffusion- (transport-) controlled dissolution (Fig. 12). However, such a formulation adds additional complications. For surface-controlled reactions at the continuum scale, it is commonly assumed that dissolved species concentrations on the reactive mineral surface are identical to species concentrations in the bulk solution (Fig. 12a), while for the shrinking core model, it is assumed that the species driving the dissolution of the mineral is consumed instantaneously when arriving at the mineral surface (Fig. 12c). However, for a transition regime, the concentration of the species driving the dissolution of the mineral at the reactive surface is unknown (Fig. 12b). Due to this additional complexity, this formulation has not been further developed to date.

The formation of secondary minerals may also affect the reactivity of primary mineral phases and their dissolution rates. The precipitation of these secondary phases on the surface of primary minerals can lead to surface passivation (e.g., Harrison et al. 2016) induced by heterogeneous nucleation (e.g., precipitation of carbonate minerals in permeable reactive barriers, precipitation of hydrated Mg-carbonates during carbon sequestration, precipitation of Fe and Mn oxides). Surface passivation is not straightforward to take into account, due to complex surface morphology of secondary phases and limited availability of data. Conceptually, the reactivity of primary mineral phases can be linked to the volume fraction of secondary precipitates. Jeen et al. (2007) used an exponential formulation to describe evolving reactivity of iron within a permeable reactive barrier as a decline in reactive surface area:

where φ_p represents the total volume fraction of precipitated secondary carbonate minerals. Jeen et al. (2007) showed a significantly better agreement with observed data using the surface passivation model than for using a two-third power law formulation (Eqn. 22). Similarly, Harrison et al. (2016) described the surface passivation using a volume-fraction threshold model, implying that dissolution rates approached effectively zero above a threshold volume fraction for secondary precipitates. Also, Daval et al. (2009) described the passivation of wollastonite dissolution by the observed precipitation of calcite, by modeling a linear decrease in surface area with respect to calcite volume fraction.

Noiriel et al. (2012) studied the evolution of surface area during carbonate precipitation, based on imaging techniques, BET analysis and a micro-continuum model. The main conclusions were that the local increase in surface area arising from mineral precipitation had a significant impact on the integrated mineral concentrations and the breakthrough curves. Modelling the experiments with a constant surface area corresponding to the measured initial value failed to represent the observed data, indicating the necessity of taking into account the evolution of surface area during precipitation processes (Fig. 13). As discussed in more detail below (see case study by Beckingham et al. 2016, 2017), mineralogical heterogeneity and the evolution of reactive surface area can also substantially affect overall rock reactivity and solute release.

In the following, three case studies are presented covering the three main aspects of evolving porous media covered in this chapter, namely evolving permeability, diffusivity and reactivity. The first case study by Poonoosamy et al. (2015, 2016) focuses on permeability evolution in response to mineral-dissolution precipitation reactions. The second case study by Navarre-Sitchler et al. (2009) evaluates the effect of weathering on diffusivity on the pore scale and the continuum scale, while the third case study by Beckingham et al. (2016, 2017) focuses on the progress and evolution of mineral dissolution reactions in a heterogeneous rock matrix. These three case studies have provided successful comparison between experimental data and continuum approaches, incorporating details from pore-scale observations.

The experimental and modeling study by Poonoosamy et al. (2015, 2016) provides an example on the effect of coupled mineral dissolution–precipitation reactions, leading to a localized decrease of porosity with observable effects on flow and solute transport.

The 2D experimental setup, depicted in Figure 14, considers a layered porous medium composed on non-reactive quartz sand in which a reactive layer of celestite (SrSO₄), composed of two different particle sizes, was embedded. A concentrated solution of barium chloride was injected at a constant rate (100 µL/min) at the bottom left corner of the experimental reactor to induce celestite dissolution and barite precipitation (BaSO₄), associated with a porosity decrease. The effect of the porosity decrease on permeability was assessed by measuring fluid pressure upgradient and downgradient of the celestite layer over time.

Tracer tests were performed before and after completion of the experiment and showed significant discrepancies due to the porosity evolution of the reactive material (not shown). Barite precipitation was recorded throughout the reactive zones, with an enrichment near the left interface between the quartz sand and the region containing fine-grained celestite (Fig. 15a). High resolution μXRF imaging presented in Poonoosamy et al. (2016) provided insights into the dynamics of the dissolution/precipitation reactions. Dissolution of the small celestite grains was rapid and barite precipitation locally filled the open pore space. Larger celestite grains; however, were passivated by barite precipitation on the surface (Fig. 15b). Image analysis showed that barite precipitation led to partial surface passivation of the celestite grains as well as cementation of the solid matrix, disconnecting initially connected pores (Fig. 15b).

This zone of higher barite precipitation led to significant modifications of permeability and the flow patterns, as recorded by the measured pressure increase (Fig. 16a), despite maintaining constant flow rates. Figure 16b shows the dissolved quantities of celestite replaced by barite precipitation for the duration of the experiment.

These experimental results were analyzed with a reactive transport model to reproduce the observed mineral dissolution and precipitation patterns and the decrease in permeability. Reactive transport analysis required the ability to model density-dependent flow and transport, the evolution of porosity and the feedback of pore clogging on permeability and diffusion. Poonoosamy et al. (2015) presented simulations that reproduced the experimental results: evolution of the chloride concentration at the outlet (not shown), the increase in pressure due to pore clogging (Fig. 16a), and the amount of dissolved celestite and precipitated barite (Fig. 16b). The impact of evolving porosity on tortuosity and permeability were modeled using Archie's law and the Kozeny–Carman relationship. To obtain agreement between the experimental results and the simulations, it was necessary to include a fine-grained celestite layer. Since a larger amount of small celestite grains are located in this zone, locally higher reaction rates were simulated which led to a more significant porosity and permeability reduction in a narrow region, which had a substantial impact on flow, matching the experimental values.

The benchmarking study of Poonoosamy et al. (2018) was based on these experiments. Several codes (CORE2D, MIN3P-THCm, OpenGeoSys-GEM, PFLOTRAN, TOUGHREACT), showed similar results regarding the evolution of porosity, mineral volume fraction, permeability and concentrations of the main elements. However, discrepancies in the spatial porosity and permeability patterns were observed. Therefore, the evolution of the flow patterns differ slightly between the codes, and discrepancies increase with time, and become more important when strong porosity changes occur (cases 3a and 3b from Poonoosamy et al. 2018).

To understand and predict the rates of chemical weathering under diffusion-controlled conditions, Navarre-Sitchler et al. (2009) carried out a combined X-ray synchrotron microtomography and tracer diffusion experiment on a sample of weathered basalt from a Costa Rican chronosequence. They used the resulting tortuosity-diffusivity model in the reactive transport code CrunchFlow (Steefel et al. 2015a) to simulate the evolution of the weathering rind over 340 thousand years (Navarre-Sitchler et al. 2011). Using a sample of basalt from a 125 thousand years Costa Rican alluvial terrace, they conducted a tracer experiment that was imaged with mXRF. Figure 17 shows the distribution of the bromide tracer (blue) after 7 days resulting from diffusion from left to right. Navarre-Sitchler et al. (2009) took the additional step of comparing these results to estimates of the diffusivity from pore-network modeling (see discussion below).

Following the tracer experiment, Navarre-Sitchler et al. (2009) made use of X-ray synchrotron microtomography to map the pore structure of samples of the weathered basalt. Using a simple implementation of thresholding to map basalt versus pores, they were able to delineate chemical weathering related macroporous zones (> 4.4 μm voxel resolution) that were connected in 3-D. Navarre-Sitchler et al. (2009) then carried out numerical tracer diffusion experiments in 3-D cubes of weathered basalt by assuming a low diffusivity of 1.75 × 10⁻¹⁴ m² s⁻¹ for the largely unconnected pore structure of unaltered basalt based on separate tracer diffusion experiments (not shown) and a free ion diffusivity of 10⁻⁹ m² s⁻¹ (corresponding to a tortuosity of 1.0) for connected pores that can be fully resolved with the 4.4 mm discretization. Implemented in this way, the results of a 3-D numerical tracer diffusion experiment using CrunchFlow is similar to what could be obtained from a pore-network model. A 2-D slice through the skeletonized pore structure of one of the weathered zones is shown in Figure 18a (basalt in blue, pores in red). Results of the numerical tracer experiment are shown in Figure 18b, with diffusion of the tracer from the bottom of the figure towards the top. Note that in these 3-D tracer diffusion simulations, only two distinct tortuosities (or diffusion coefficients) are used and they are based on the segmented 3-D porosity map of the weathered sample.

The segmented microtomography data on the weathered basalts was then used to quantify connected versus total porosity in the samples using the burning algorithm Percolate (Bentz et al. 2002). The data on all Costa Rican basalt samples analyzed indicate that connected porosity at > 4.4 µm increases markedly at about 9% total porosity (Fig. 19a).

The effective diffusion coefficient can be described with a modified Archie's Law (Eqn. 15) model that incorporates a critical porosity threshold ϕ_crit value of 9% (Eqn. 16), below which the porosity is considered to be largely unconnected. The parameter a in Equation (15) is taken as 1.3 while the parameter n is assumed to be 1.0 for 3-D volumes measuring 220 mm on a side (Navarre-Sitchler et al. 2009). Figure 19b compares experimental points to this relationship with the usual Archie's law described by Equation (13), which used the total porosity as mapped with the X-ray synchrotron microtomography with a 4.4 mm voxel resolution.

Based on this estimate of the effective diffusion coefficient shown in Figure 19b and captured in Equations (15) and (16), it was possible to compare the results of reactive transport simulations to observations of weathering rind thickness of basalt clasts in the Costa Rican alluvial terraces. As shown in the study, it is only possible to capture the four characteristic features of the weathered basalts, including (1) the mineralogy of weathering products, (2) weathering rind thickness, (3) the coincidence of plagioclase and pyroxene reaction fronts, and (4) the thickness of the weathering rinds, by including the modified Archie's Law formulation of effective diffusion coefficient. In particular, the coincidence of the pyroxene and plagioclase fronts, which have differing intrinsic solubilities and reactivities, is possible only with a threshold model of the kind implemented here (Navarre-Sitchler et al. 2011). In an unmodified Archie's Law formulation, the pyroxene and plagioclase fronts would move at different rates, so it is because of the generation of connected porosity as a result of mineral dissolution (weathering) that the mineral fronts advance in lock step. Figure 20 shows the position of the weathering fronts (data and solid lines model fit) as quantified using the mass transfer coefficients for calcium relative to immobile titanium at 240, 125 and 7 thousand years for quaternary terraces 1, 2 and 3, respectively.

A case study for evolving reactivity is provided by Beckingham et al. (2016, 2017) on highly reactive volcaniclastic sediments from the Nagaoka pilot CO₂ injection site in western Japan. The investigators carried out batch and coreflood experiments on a sample of sediment from the site to determine the time evolution of the sediment and fluid as a result of injection of CO₂-saturated fluid. They combined the experiments with multi-scale image analysis using synchrotron X-ray microCT, SEM QEMSCAN, XRD, SANS, and FIB-SEM to map the initial mineralogy and pore structure of the sediment (Fig. 21). The mineral mapping provided the initial conditions for continuum reactive transport modeling with the code CrunchFlow (Steefel et al. 2015a), and the simulation results were compared against effluent chemistry over the course of the batch and coreflood experiments.

The mapping with a multi-scale image analysis described provides both the modal mineralogy of the complex sediment and the co-location of connected porosity and reactive phases (Beckingham et al. 2016), thus defining more rigorously its bulk effective reactive surface area. In addition, the imaging was used to create a grain size distribution of the most reactive phases that could be used directly in the modeling. This was considered necessary because the fluid chemistry in both the batch and coreflood experiments evolved significantly over time, with early times showing very high cation (Ca²⁺ and Mg²⁺) and silica concentrations that were as much as 10× higher than what they were at the end of the ~650 hours experiment.

The hypothesis was that the fine fraction of phases like pyroxene and plagioclase could account for the high reactivity early in the experiment because of their high specific surface area (Fig. 22). Using the mapping to constrain the volume fraction of the pyroxene and plagioclase (Table 1), it was possible to use the reactive transport simulations in both the batch and coreflood experiments to determine whether the pronounced evolved reactivity could be attributed to the early reaction and then disappearance of the finest mineral fractions. While the trends were correct for pyroxene and plagioclase, it was found that it was also necessary to include a reactive glass phase to account for the high concentrations and reactivity (Fig. 23).

For the aforementioned reasons and given the large variety of pore structures and heterogeneities, it is evident that obtaining continuum scale relationships for effective transport properties and reactivities in porous media, in particular under dynamically evolving conditions, constitutes a major challenge.

Accurate simulation of breakthrough curves, mineral abundances and concentration distributions relies heavily on the adequate evaluation of the evolution of the porous medium's effective porosity, transport and reaction properties. As discussed above, a wide range of feedback relationships exists, relating transport parameters back to porosity; however, these relationships provide diverse predictions of the transport properties, in some cases differing by several orders of magnitude. In addition, it has been shown that identical systems might evolve differently as a function of transport regime, as well as the ratio between the time scales of reaction and transport. Furthermore, it is also common that relationships between transport parameters and porosity themselves evolve, as the pore structure is modified. Keeping in mind that the evolution of porous media is predominantly affected by processes occurring at the pore scale, and that mineralogical compositions and reactivities can vary widely, it becomes clear that the predictive capabilities of relationships describing the evolution of porosity, tortuosity, permeability and reactivity are limited. It is necessary to constrain these relationships by observational data on a case-by-case basis, although we are confident that some generalizations can eventually be drawn in terms of processes and materials. To this regard, the presented case studies show adequate approaches on how to incorporate pore-scale observations into the continuum approach. Although some theoretical underpinnings exist, relationships to describe the evolution of transport and reactivity parameters remain empirical in nature and are case specific. However, despite these shortcomings, empirical relationships, often in the form of power laws, have proven useful to provide a quantitative interpretation of observational data under dynamic and evolving conditions. Examples include problems involving diffusive mass transport (Seigneur et al. 2017; Georget et al. 2018), changing permeability (Mok et al. 2002; Luhmann et al. 2014) and evolving reactivity (Jeen et al. 2007; Noiriel et al. 2009). Although empirical in nature, these relationships provide insight and an opportunity for quantitative interpretation of feedback mechanisms.

It is relatively trivial to implement these empirical relationships in continuum scale reactive transport codes. However, the simulation of extreme cases of evolving porous media poses a challenge. On the one hand, complete clogging causes the disappearance of the water phase, which most geochemical solvers used in simulators are not able to deal with. Typically, reactive transport codes deal with the clogging problem by setting an arbitrary porosity threshold value, below which groundwater flow and solute migration does not occur. Although benchmarking exercises have been successful at verifying the correct implementation of permeability–porosity and tortuosity–porosity relationships, obtaining agreement between different reactive transport codes under clogging conditions has been more difficult (Xie et al. 2015). On the other hand, simulations can yield a substantial increase of porosity that, in fact, might not be physically representative, since the solid matrix would eventually collapse. This constitutes a mechanical issue which at this point in time is not commonly considered by reactive transport simulators.

Modelling of unsaturated cases, whether through Richards equation (Eqn. 2) or two-phase approaches (Eqn. 1), requires the use of relationships to take into account the material's hydraulic and transport properties. The models of Brooks Corey, van Genuchten provide the water retention characteristics (capillary-pressure), as well as the relative permeabilities of the fluid phases. These models involve several different empirical parameters. The latter have been assessed for different porous media and general understanding of their dependence on physical properties has been acquired (e.g., Zhang et al. 2017 for the dependence on grain sizes). Della Vecchia et al. (2015) also described a model linking pore size distribution to water retention models, specifically for clays. Similarly, Suazo et al. (2016) studied the evolution of the water retention curves during cement hydration. However, only a limited amount of research has been dedicated to the prediction of how these parameters evolve as the porous media undergoes chemical perturbations. Furthermore, no clear relationship between the unsaturated parameters and porosity has been described. For these reasons, most reactive transport simulations consider these parameters to be constant.

The reactive transport problem is commonly solved by either the global implicit approach (e.g., Crunchflow, MIN3P), the sequential-iterative approach (Hytec), or the sequential non-iterative approach (e.g., PHT3D, HP1) (Steefel et al. 2015a). Independent of the numerical methods used, permeabilities and reactivities are typically updated in a time-lagged fashion (e.g., Mayer et al. 2002 for MIN3P). While this is not necessarily the case for porosity and diffusion coefficient (Lagneau and Van Der Lee 2010), this implies that at each time step, small mass balance errors are introduced.

An additional complication occurs due to the feedback of evolving porous media parameters between reactive transport and flow. In most models, groundwater flow and multicomponent reactive transport are solved sequentially. This approach is adequate if there is no feedback of the reactive transport solution on the flow problem. However, as depicted in Figure 1, this is not the case for evolving porous media, since a change in porosity affects the storage term in the flow equation (Eqn. 1 and Eqn. 2) and the flux term, by modifying the permeability. The decoupled treatment of groundwater flow and reactive transport makes it necessary to take relatively small time steps to avoid the introduction of mass balance errors. However, recent efforts to recouple multiphase flow and reactive transport in case of variable porosity have been implemented in Hytec (Seigneur et al. 2018).

As outlined above, the concepts of surface area, tortuosity and (effective) porosity play a significant role in evolving porous media. Usually, tortuosity and surface area are linked to diffusion and reactivity. Empirical relationships are generally used to update their evolution. Relationships for different evolving parameters are commonly independent from one another, despite the fact that they are inherently correlated, as suggested by numerical studies by Koponen et al. (1997). While general and universal correlations between these relationships are difficult to obtain, common trends can be estimated. Typically, as porosity increases, surface area tends to decrease (Koponen et al. 1997), while diffusivity (and hence, as defined in this chapter, tortuosity) increases (Knackstedt and Zhang 1994; Koponen et al. 1997). Saripalli et al. (2002) have also shown that tortuosity can be related to surface area (a higher surface area leads to a lower diffusivity/tortuosity).

Most of the efforts to estimate effective transport properties such as tortuosity and permeability have focused on developing relationships invoking porosity, tortuosity and/or surface area. While permeability and diffusivity represent fundamentally different transport processes, their dependence on tortuosity (hence surface area) and on the pore throat size are similar.

As was observed for the simulation of clogging reactions during cement carbonation, systematic errors are encountered when attempting to reproduce experimentally observed alteration depths. A failure to reproduce observed degradation depths, while managing to reproduce breakthrough curves, might be an indication of the need to work towards a more integrated framework of evolving parameters.

While the mathematical description of such an integrated approach is not yet determined, it is evident than an approach linked to porosity evolution alone is not suited to properly describe the processes which occur at the pore-scale and the continuum scale. Instead, geometrical considerations about the structural organization of the pore scale might be more appropriate. As an example, the surface area in itself carries more information about the potential interactions between the fluid and the solid matrix. On the one hand, high surface area indicates tortuous pathways which are linked with lower transport properties. In addition, reactivity is also intrinsically linked to the surface area.

Let us consider the observed inverse correlations between permeability and porosity. An integrated generalized approach based on surface area might be able to reproduce these results in a better way than one based on porosity alone.

In recent years, imaging techniques have helped to improve our understanding of pore structure and its evolution, as well as the effective transport properties. Such techniques, combined with high-resolution of reactive transport at the pore scale (Molins et al. 2012; Steefel et al. 2013; Molins et al. 2014) are destined to be used increasingly, as they allow for a process-based assessment of the evolution of the pore structure. Imaging and pore scale modeling techniques have already provided outstanding results and insights. Molins et al. (2017) reproduced dissolution processes in a multi-mineral fracture zone by incorporating the heterogeneity of the spatial distribution of minerals. In another study, it was possible to reproduce the wormholing processes by direct simulation at the pore scale (Molins 2015). In addition, pore scale modeling has provided a means to document and quantitatively assess the effect of dissolution and precipitation reactions on the pore structure, as they occur within a porous material (Molins 2015; Noiriel 2015; Bultreys et al. 2016; Pereira Nunes et al. 2016). Future improvement in the resolution of imaging techniques might even make it possible to characterize porous media at even smaller scales and improve pore scale models built from these methods. However, it is unlikely that pore scale simulations will be able soon to deal with the large time and spatial scales required for engineering purposes.

Pore scale modelling is a promising tool and may eventually allow the derivation and upscaling of transport parameters for use in a continuum approach. This modeling technique may be able to help overcome the biggest gaps which still exist for the continuum approach, in particular when applied to evolving porous media. In the future, it is likely that combined approaches (micro-continuum) will provide a tool to bring the knowledge acquired from the pore space to the macroscale (Steefel et al. 2015b).

Biogeochemical reactions in porous media, involving mineral dissolution and precipitation, can lead to significant modifications of the pore structure with a potential dramatic evolution of the material's properties in a wide variety of natural and engineered systems. Reactive transport modelling is a useful deterministic tool to investigate these strongly coupled nonlinear effects. However, an accurate description of these phenomena requires a thorough understanding of their dynamics at the pore scale.

The last decade has provided major improvements with pore scale imaging and simulation. First, imaging techniques have allowed to reveal more detailed information about pore structures and their evolutions at a sub-micron level. Second, the parallel increase in computing capabilities opened the door to pore-scale simulations which can now be performed on 3D image samples. These tools, used together, provide precious data for developing an improved understanding of how porous media evolve under different conditions.

However, there are still major shortcomings to overcome. First, the resolution of imaging techniques is currently insufficient to provide adequate insights into the 3D-structure of fine grained materials. Indeed, porous media in which nanometric pores constitute the building blocks of the pore structure (clays, cement), cannot be accurately resolved by these imaging techniques at this point in time. The second problem, potentially more challenging, is that reactive transport simulations relevant to engineering problems usually involve scales (m–km) which are orders of magnitude larger than the scales at which processes relevant to porous media evolution occur (nm–µm). Continuum approaches for these engineering purposes likely still have a promising future, despite their inherent limitations.

Some significant limitations remain. The clogging phenomenon still constitutes a challenge due to the inherent complexity of capturing the effect on transport processes when near zero porosities are approached and the associated numerical issues. Furthermore, there is still a lack of experimental data to calibrate reactive transport models. Additionally, the development of fractures due to internal stresses (either due to heat generation, mineral dissolution or precipitation) have been observed in different experimental setups. Unfortunately, a full coupled problem involving heat and mechanics is not available.

Despite obvious challenges inherent to the continuum approach regarding the evolution of relevant parameters describing a porous medium's properties (tortuosity, permeability, reactivity and surface area), quantitative reactive transport analysis of evolving systems has been successful. Additionally, significant knowledge has been acquired through recently developed micro-continuum approaches. Empirical relations for transport (such as Archie's law, the Kozeny–Carman relationships or various power laws) or reactivity (threshold passivation or shrinking-core model) are not meant to be physically representative of pore-scale process. It must be acknowledged that it is not realistic to expect that using these relationships would lead to a fully process-based assessment of experimental data. Nevertheless, we now possess sufficient knowledge about how certain processes can be modelled through the use of these empirical relations.

The contributions by Nicolas Seigneur were supported by a postdoctoral fellowship through the TERRE-NET program, a strategic network funded by NSERC (Natural Science and Engineering Research Council of Canada).

K. Ulrich Mayer would like to acknowledge funding by NSERC through a Discovery Grant.

The contributions by Carl I. Steefel were supported by the Director, Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231 to Lawrence Berkeley National Laboratory.

We would like to acknowledge Jenna Poonoosamy for providing additional material for Figure 15b. We would like to thank Tom Al (UOttawa) for providing the images for Figure 3.