The field of reactive transport lies at the intersection of several disciplines in the Earth and Environmental sciences, including hydrology, geochemistry, biology and geology. The processes in natural and engineered media that are the focus of study of these disciplines take place over a wide range of spatial and temporal scales. Specifically, geological media are characterized by their physical and mineralogical heterogeneity at spatial scales from nanometers to hundreds of meters and beyond. Flow and advection of solutes take place at the scale of individual pores but are commonly represented at the Darcy scale where the porous medium is treated as a continuum. A large contrast is often observed between fluid residence times in regions of enhanced permeability such as fractures or macropores and less permeable media where diffusion may be the dominating solute transport process. Understanding of reactive processes, including those mediated by microorganisms, is often developed at the molecular scale in the laboratory but their impact in the environment is observed at larger spatial scales.

In addition to considering the scales of the individual processes, reactive transport must also consider how these different scales interact with one another to give rise to the overall coupled behavior. In fact, in many instances considering the processes at the observation (or native) scales has limited applicability in subsurface environments. For example, reaction rates derived from laboratory studies show large discrepancies from those observed in natural environments, where transport processes and accessibility to reactive areas control effective rates. Reactive transport models, thus, even in a simple form, must make assumptions regarding the scales associated with each process and how they interact with each other. Implicit in any model is also the assumption that the models for each process are applicable at the same spatial scale as the other processes represented. For example, local geochemical equilibrium may only be assumed where reaction rates are faster than transport rates such that the solution reaches equilibrium over a characteristic spatial scale.

Reactive transport modeling, as a tool to integrate knowledge and develop mechanistic understanding, seeks to incorporate improved process model representations that reflect our advances in fundamental understanding (Druhan and Tournassat 2019 and references therein). The multiscale nature of reactive transport is one of the most prominent aspects in models and, hence, modeling approaches that address this multiscale nature are an increasingly important component of the toolset needed by researchers (Scheibe et al. 2015a). In particular, multiscale approaches make it possible to incorporate process representations at the appropriate native scale in models intended to simulate the coupled problem at a different spatial or temporal scale. A variety of approaches have been brought to bear that range from conceptual to mathematical to numerical. Specific goals of multiscale models are also diverse and include using the appropriate coupling between processes, capturing the processes at the relevant scale in different regions, capturing the physical and mineralogical heterogeneity at multiple scales or incorporating fine-scale information in larger-scale applications. Ultimately, there is the need to identify what processes and at what scale are controlling overall system behavior, and hence the appropriate spatial and temporal scales to represent each process.

Multiscale modeling is a very broad topic with applications across many disciplines (Tomin and Lunati 2013, 2016; Scheibe et al. 2015a; Amanbek et al. 2019). In this chapter, we specifically review multiscale approaches for reactive transport modeling from the conceptual and mathematical perspectives. Many of the approaches have also been used in the individual disciplines reactive transport draws from. They are here discussed in a general manner here but also specifically in relation to reactive transport applications. We include approaches where the multiscale nature is reflected in a continuum-mechanics-based model and we discuss numerical aspects of these approaches where needed. Approaches that incorporate upscaling procedures in the numerical solution process such as multiscale finite element or finite volume methods, or based on numerical upscaling however are not included, e.g., Efendiev and Hou (2009).

We start by establishing the equations that describe the processes of interest at a single scale and discussing the multiscale aspects associated with process coupling at a single scale. We relate these equations to the two scales commonly identified in porous media—the pore scale and the Darcy scale—but are generally applicable at a range of spatial scales, from fluid in pores to streams and rivers. Derivation of effective models by upscaling pore-scale equations to the Darcy continuum scale is used specifically to motivate the need for multiscale approaches. First, we describe approaches that use different scale representations in different regions in the domain. We continue with approaches based on the existence of two or more porous continua in the same region of the domain. We give examples of the use of the multiscale approaches described in selected literature applications before making some concluding remarks.


Reactive transport models simulate flow, solute transport and geochemical reactions. In this chapter, for simplicity, we will consider only single-phase (aqueous) flow and transport. These processes are typically described by two sets of equations, one for the conservation of mass and momentum of the fluid and the other for the conservation of mass of the reactive components. The form of these equations depends on continuum of reference for which they are written: a fluid continuum, a solid continuum or a continuum that includes both fluid and solid phases. For example, flow in streams and rivers may be represented by considering the fluid phase as the continuum of reference. In porous and fractured media in the subsurface, when the fluid and the solid are both treated as separate continua, we refer to the scale of observation as the pore scale, while when the porous medium is the continuum, we refer to it as the Darcy scale. For convenience, in the derivation that follows we will focus on porous and fractured media, but these single-scale equations can be read more generally as applicable at a range of spatial scales. Further, we will not consider here the scales where the medium is not characterized as a continuum. For example, we do not discuss characterizations at the molecular or atomistic scale, or the organism level in the case of microbial processes.

Separate fluid and solid continua: Pore-scale equations

When the individual pores are represented explicitly, and the solid–fluid interfaces are the boundaries of the domain considered (Fig. 1), flow may be described with the Stokes equations


where u is the fluid velocity (with |u| =0 on the solid–fluid boundary), and m, and p are the fluid viscosity, and pressure, respectively. The equations that describe the mass balance of chemical species subject to advective–diffusive transport and heterogeneous reactions at the fluid–solid surface may be written as


where c is the solute concentration, D is the diffusion coefficient of the solute in the solution, and r is the surface reaction rate. Equation (4) expresses the mass balance at the fluid–solid interface for the aqueous species involved in the heterogeneous reaction, where n denotes the unit normal pointing from solid to liquid.

Single porous continuum: Darcy-scale equations

When the porous medium is treated as a continuum, i.e., when within an elementary representative volume (or REV) properties that described the medium such as porosity (q) or permeability (k) may be assumed constant (Fig. 1), flow in porous media can be described by:


where q is the Darcy velocity vector, which is calculated with Darcy's law (Eqn. 6), q is the porosity (water content in fully-saturated conditions), r is the fluid density, g is the gravitational constant, and z is the vertical coordinate. The equations that describe the mass balance of chemical species subject to advective-dispersive transport and heterogeneous reactions in porous media may be written as


where C is the concentration at the Darcy-scale, D* is the effective diffusion/dispersion tensor, and R is the bulk reaction rate. In the view presented in this section, the properties that characterize the porous medium, i.e., θ, k, or D*, or the bulk rates are assumed known (for example, empirically) and applicable at this scale. In general, however, they encapsulate information of the processes that take place at the scale of individual pores. In the section devoted to upscaling, we discuss how one may formally derive these parameters from the pore-scale description.

Multiscale aspects of process coupling

The flow equations and the reactive transport equations are in general coupled via the composition of the fluid, as well as hydraulic and geometric properties of the porous medium. For example, changes in the fluid composition caused by geochemical reactions affect the fluid density, or dissolution–precipitation reactions change the pore space geometry which in turn affects the fluid velocity. The evolution of these properties is in many applications relatively slow compared to flow and transport and the coupling is assumed weak. In some other applications, the feedback processes are significant and are addressed specifically in the chapter devoted to porous media evolution in this volume (Seigneur et al. 2019, this volume).

Here we will briefly focus on the coupling of time scales associated with transport and reactions and how they affect the process representation at a specific spatial scale. For convenience, we re-write Equation (7) generically as


where ℒ() is the transport operator and ℛ() the reaction operator. Transport is important as it provides the driving force for reaction but also because it provides a characteristic time scale to which the time scale of the reaction is compared. In an open system, as implied by Equation (8), if the characteristic time of transport (for a given characteristic length scale), τ, is larger than that of the reaction, τ, the solution reaches equilibrium with itself or with a mineral phase. In Darcy-scale models, if this characteristic length scale is that over which the REV is defined, this makes it possible to assume local equilibrium. In the local equilibrium assumption (LEA) (Lichtner 1996), the rate of reaction is thus determined by the rate of transport of matter across the boundaries of the domain (Fig. 2a).

In pore-scale models (Eqns. 34 may also be written in a form similar to Eqn. 8), the concept of local equilibrium is different (Lichtner 1996). At the Darcy scale, geochemical equilibrium is attained at the REV scale which includes many pores. At the pore scale, the detailed pore space geometry plays a role. Although equilibrium may be attained at mineral surfaces, geochemical gradients may still be present in individual pores (Fig. 2b).


The macro-(Darcy) scale model in Equations (56) and (7) may be formally derived by scaling-up the micro-(pore) scale Equations (12) and (34). Upscaling has been the subject of intensive research for at least 50 years, with volume averaging being a commonly used approach for this purpose. In classical averaging theory, concentrations and fluxes are averaged over an REV composed of fluid and solid phases. The averaging of the pore-scale equations makes it possible to write conservation equations for these averaged quantities (see, e.g., Gray and Miller 2014). In these equations, however, new terms appear that still depend on pore-scale quantities and hence for which closure relations must be postulated. The mathematical theory of (periodic) homogenization is used to solve the closure problem (Hornung 1997).

In the most simple form, periodic homogenization relies on the spatial periodicity of the domain Ω, which is conceived to be composed of shifted and e-scaled copies Yɛ of an REV Y (Fig. 3). The REV Y is made up of solid (Ys) and liquid (Yl), separated by an interface G. Accordingly, the porous media domain Ωɛ is composed of solid phase, with Ys,ɛ being the union of the shifted and e-scaled Ys and the liquid phase. The subdomain Yl,ɛ is defined similarly and assumed to be connected such that flow can take place. Let's assume a coefficient (φ) that oscillates with a periodicity of length e such as the diffusion coefficient D, with D > 0 in the fluid and D = 0 in the solid. The aim of homogenization is to consider ɛ → 0 (i.e., to zoom out) in order to see what description holds for the emerging homogeneous medium. Formally, this can be done assuming a two-scale asymptotic expansion for all unknown quantities φ in the form


where y : = x / ɛ is a fast spatial variable in the sense that y covers Y, if x covers Yɛ. Depending on the problem, the interplay between the processes in e and separation of scales, it is possible to derive equations in which only the Y-averaged 0-th order terms appear. This is done by inserting the expansion in the governing (pore-scale) equation such that each e-power gives rise to a new equation for the coefficient that is being investigated.

For example, one can derive the upscaled, Darcy-scale form of the diffusion equation (i.e., Eqn. 7 with q = 0 and no reaction term) from the pore-scale diffusion counterpart (i.e., Eqn. 3) with |u| = 0). Solving the closure equations on Y numerically makes it possible then to obtain the effective diffusion tensor D*, which encodes the information about the pore geometry, see, e.g., Ray et al. (2018). The closure equations are typically solved using idealized or randomly generated porous geometries (Fig. 4). From these calculations, it is possible to describe the dependence on pore geometry solely by a macroscopic parameter such as porosity (Fig. 5).

The above procedure can be accomplished when there is scale separation (i.e., the micro-and macro-scales can be described independently from each other) but in general it also depends on the proper scaling of the parameters in the e-problem. For example, to deduce Darcy's law from the Stokes equation at the pore scale (i.e., with a computable permeability), a scaling of the viscosity m to em is required. Further, when several different processes are in play (e.g., advection and diffusion), characteristic numbers such as the Péclet number (Pe) may be needed in the efficient coefficients, e.g., in the derivation of mechanical dispersion (Mikelić et al. 2006). For reactive transport problems, efforts in upscaling have focused on heterogeneous reactions. Here we present a short derivation of upscaled relationships for (first-order) sorption reactions and dissolution-precipitation reactions first with non-evolving geometries and then with evolving geometries.

Upscaling interfacial reactions

We consider a surface reaction such that Equations (34) may be written on Γɛ as

cmt=(ucmDcm)in Yl,ε
Dcmn=εr on Γl,ε

where cm is the concentration in the fluid, and the rate is expressed as an exchange between the fluid and the surface concentration, r = α(cmcim), with the immobile concentration (cim) being a surface concentration such that


where α is the rate of mass transfer between the bulk fluid and the surface, with α > 0. A noflow boundary condition (u = 0) is assumed on Γl, ɛ. The upscaled model then takes the following macroscale form


where σ is the specific surface (σ = Γ/Y, which along with θ and D*are computed from the pore geometry, and the macro-scale rate is described by


where micro- and macro-scale problems are separate. However, if surface diffusion is considered at Γɛ such that the micro-scale Equation (12) is written as


where Ds, is the surface diffusion coefficient, it is not possible any longer to obtain the average of cim and a coupled macro-scale/micro-scale model is obtained

cim(y,t;x)tDsΔycim=α(kcmcim) for yΓ,xΩ

where the micro problem depends at every point x on the macroscopic concentration Cm(x,t). (Note that the emerging equations for Cm and Cim = Cim (y,t;x) read as Eqn. 35, with G in lieu of Wx).

Upscaling interfacial reactions with evolving geometries

Precipitation–dissolution reactions are of interest as they change the pore geometry and can have a positive feedback on flow and transport processes, leading in some case to clogging or wormholing. In the classical approach this is reflected in the change of porosity as described by Ray et al. (2015)


where ρs is a surface density used as conversion factor between mass and volume. To close the model, the specific surface σ must be related to θ, which requires assumptions on the evolving micro geometry. In the case of known geometries (e.g., spheres, cubes, etc), this leads to the commonly used relation (Lichtner 1996)


where Cmin is the mineral concentration and β is a factor to convert from units of mass to volume

The evolution of the interface as a result of the reaction can be described with a sharp interface approach (e.g., a level set) or approximated by a phase field model (Bringedal et al. 2019). The normal component of the velocity of the interface is denoted by vn


This velocity is a function of the reaction rate R and the density of the precipitate ρs. Considering the reaction to be a single-component reaction for simplicity, Equation (11) can be replaced with


where the precipitate can form the outer surface of the solid grain. By upscaling we can again recover Equation (13), but here the closure comes from a detailed description of the interface via the evolution of a level-set function L0 in the REV Y, described by (Ray et al. 2015)

tL01ρR(Cm(x,t))|yL0|=0 in Y

yielding θ, σ, D needed in Eqn. 13 but dependent now on the macroscopic concentration Cm(x, t). Thus, we obtain again a coupled micro-macro model. The solution of Equation (21) is not straight-forward and can be derived by a two-scale asymptotic procedure (Ray et al. 2015).

Extension of this approach to multicomponent problems that include homogeneous and heterogeneous reactions, in equilibrium or kinetic, is still in the early stages of development. While homogeneous or sorption reactions are tractable, development of passivation layers in systems with multiple mineral reactions, e.g., Daval et al. (2009), is considerably more difficult to handle in this framework.


The ability of pore-scale models to explicitly resolve individual pores make them suitable to simulate flow and reactive transport without using bulk parameters to characterize the medium. The computational cost of pore-scale simulations, however, is very high if one wants to cover volumes of porous media large enough for relevant applications. As noted in the previous section, the macro (Darcy)-scale problem can only be formulated separately from the micro (pore)-scale problem under a number of simplifying assumptions, see also Battiato and Tartakovsky (2011). Rather than formulating and solving closure equations for periodic unit cells of idealized geometries in these cases, one may want to solve the pore-scale problem directly. In applications where the pore-scale characterization is needed only in certain regions rather the entire domain, an attractive approach is then to combine a pore-scale description in these regions and revert to a Darcy-scale description elsewhere. Broadly, two approaches have been used for this purpose: hybrid models and the Darcy–Brinkman–Stokes approach.

Hybrid multiscale models

Hybrid multiscale models combine different scale representations in a single simulation (Battiato et al. 2011; Roubinet and Tartakovsky 2013; Yousefzadeh and Battiato 2017). In this approach, the domain is divided in two or more regions where different scale representations are used (Fig. 6a). Because often this implies that different process models, and thus potentially also the numerical solution method, are used in each of the regions, hybrid models are sometimes known as multi-algorithm or algorithm refinement models. An advantage of hybrid models from the multi-algorithm perspective is that it is possible to consider different spatial and temporal discretization in different portions of the domain, or the use of different dimensionality for each of the sub-domains (e.g., Fig. 6b).

A hybrid model of reactive transport on a domain Ω composed of a pore-scale domain Ωp and a Darcy-scale domain ΩD, such that Ω = Ωp∪ ΩD, entails the solution of Equations (14) in Ωp and Equations (57) in ΩD. The pore-scale and Darcy-scale simulations are coupled by enforcing the continuity of mass (concentration) and mass flux (its normal component) along the interface G between Ωp and ΩD:


where Qp(x, t) and QD(x, t) are the normal components of the pore-scale and Darcy-scale mass fluxes, Qp(x, t) = ucDc and QD = qCD*∇C, respectively.

From a numerical perspective, hybrid models need to consider the interface between the subdomains explicitly and coupling the sub-problems at the interface. Typically, concentrations and fluxes are used as coupling unknowns for Equations (22) and (23). This adds some complexity in the implementation of hybrid models, especially when the dimensionality or discretization on both sides of the interfaces are different. An example of a numerical method designed to handle these exchanges in a general and flexible way is the mortar method (Balhoff et al. 2008; Mehmani et al. 2012). In this method, coupling between sub-domains is accomplished by using finite-element (FEM) spaces to determine interface conditions (Fig. 7). Using the flow problem as example (with pressure being the unknown of the problem), the pressure field in the mortar space (noted by p) is a linear combination of finite element (FE) basis functions (ϕi):


The basis functions may be constant, linear, quadratic, or higher order functions. The solution is obtained by determining the coefficients (ζi) that result in matching of fluxes at interface. The mortar solution then describes the pressure field only at the interface which is used as a boundary condition and projected onto the individual sub-domain. The sub-models are then solved using the appropriate algorithm.

Darcy–Brinkman–Stokes approach

An approach for combining pore- and Darcy-scale representations that has received increasing attention is the one conceptualized by the Darcy–Brinkman–Stokes equation (Golfier et al. 2002; Popov et al. 2009; Gulbransen et al. 2010; Yang et al. 2014; Soulaine and Tchelepi 2016; Soulaine et al. 2017). Darcy–Brinkman–Stokes describes flow in open pore space and in a porous continuum with a single equation (written here to recover the transient incompressible Navier–Stokes equations in the fluid domains in lieu of Eqn. 2):


where ɛ is the porosity of the medium, with ɛ =1 in the pore spaces, 0 < ɛ < 1 in a porous continuum and ɛ = 0 in the solid phase. The permeability k of the medium requires a constitutive relation linking it to ɛ, for example the Kozeny–Carman equation. As a result, the terms associated with porous-media flow become negligible in the pore scale, while the terms associated with pore-scale flow become negligible in the porous continuum or the solid phase (Fig. 8). The transport of the aqueous species is described by a locally averaged equation


where dissolution is here described as a source-sink term (R) as in Equation (7).

The use of a single equation simplifies the numerical implementation of this method. Further, it does not require explicit representation of an interface between pore-scale and Darcy-scale, which is especially convenient in problems with evolving media. In fact, the distinction between pore-scale and Darcy-scale is only conceptual. In the pore-scale limit, i.e., when e and k are very small, Equations (25) and (26) recover the pore-scale description in Equations (13) and it may be used as a pore-scale method (Golfier et al. 2002; Soulaine et al. 2017).


Natural porous media are characterized by the existence of porosity at multiple spatial scales. The conceptual model presented in Figure 1 only considers inter-granular porosity, that is, porosity available between solid grains where in general fluid velocities are appreciable. We update here this model to include intra-granular porosity, that is, porosity that exists at small spatial scales in regions, which we will for convenience refer to as aggregates, where fluid flow is in general very slow and transport processes are dominated by diffusion (Fig. 9). In this media, the mass of each constituent is distributed between inter-porosity and intra-porosity, or mobile and immobile regions, between which mass is exchanged.

Multi-rate models

The concept of mass exchange between mobile and immobile regions is analogous to considering heterogeneous reactions between a mobile and immobile species, such as in sorption reactions, see Equations (1012). We will use this idea here to present single- and multi-rate models, see, e.g., Haggerty and Gorelick (1995), Hollenbeck et al. (1999). We start by writing the mass balance of a species as (omitting the species index for simplicity)


where Cm and Cim are the mobile and immobile concentrations of the species, where we assume the immobile concentration is not subject to transport ℒ(). The mobile (θm) and immobile (θim) porosities are the corresponding conversion factors to normalize the concentrations to the volume of an REV (if Cim is a surface concentration, then θim = ρb, with ρb being the bulk density). Equation (27) needs a closure relation that links the two concentrations, which may be a quasi-stationary approximation


describing a linear equilibrium sorption reaction (κ ≡ Kd), or a first-order kinetic exchange such as


Integration of Equation (29) makes the memory effect of the linear mass transfer process evident


For more flexibility in the description of this memory (or tailing) effect, Cim can be further subdivided in different fractions (βk) with different kinetic behavior in a multi-rate formulation. Equations (27) and (29) can be generalized to


with βk0,Σk=1Nβk=1.

Applying (30) to each of the fraction and expanding the sum in (31), one can see that the exponential kernel in the memory term becomes multimodal. The model can be also used to describe physical non-equilibrium, where now the mass transfer is between the mobile region and each of the fractions of the immobile region (or classes of micro-porosity, each giving a different memory effect). All these fractions (or sorption sites in the sorption interpretation), exchange mass with the mobile region directly based on the same Cm.

A generalization of Equations (31) and (32) is given by extending them from a finite number of fractions (i.e., sorption sites/micro-porosity classes) to an infinite number (or a continuum) of them by re-writing them as


with 0β(k)dk=1

Extension of multi-rate mass transfer models to multicomponent reactive transport is possible; see, e.g., Donado et al. (2009) where the equilibrium is assumed in each the mobile and immobile domains such that Equations (27) and (30) may be written in terms of total concentrations and speciation calculations performed once the equations are solved for these total concentrations.

Multi-continua models

Multi-rate models are derived for two continua: a mobile and immobile region. This can be generalized for any number of continua. One can transition from multi-rate models to multi-continua models (MC). We can start by assuming that in one aggregate class the diffusive transport within the aggregate is taken into account. For the ease of notation, we consider only one rate (N = 1). The averaged concentration Cim = Cim(x,t) is now replaced by cim = Cim(y,t;x), y ∈ Ωx where Ωx is a representative aggregate and the total mass conservation from (31) reads


Equation (31) is substituted by a partial differential equation on Ωx

dimycimny=α(Cm(x¯,t)cim(y,t;x¯)) for yΩx¯

where dim > 0 denotes the molecular diffusion coefficient in the aggregate and ny is the unit normal at ∂Ωx pointing out of Ωx.

If we have several points x = xk, where the detailed dynamics must be considered, we go from a two-continuum/region model to a multiple continuum/regions models with a considerable increase in numerical complexity. The extreme is to do this for every (discretization) point x = x ∈ Ω. In this case, we can again arrive at a micro-macro model as in Equations (13) and (16). To avoid enormous numerical complexity in simulating such a model, on can try solving Equations (36) and (37) analytically with appropriate initial conditions. This is possible in special cases, e.g., when Ωx is a sphere. This solution representation, which depends on cm(x, t), may be viewed as analogous to the multi-model extension of Equation (30), where now the integral in (30) has to be substituted by an integral of the type


to close the relation in Equation (35), with cm(y, 0, x) = 0 for simplicity.

Multiple interacting continua

The multiple interacting continua (MINC) is a multi-continuum approach developed specifically as a discretization method for fractured media. In this approach, the matrix is represented by multiple continua, each one further away from the fracture continuum. The idea is that changes in fluid conditions propagate more slowly in less permeable matrix blocks compared to the smaller fracture volumes. This approach allows for accurate resolution of the gradients in pressures and concentrations into the matrix. It is distinct from the MRMT in that matrix continua are connected in series to account mass transfer between, while in the MRMT mass transfer is always between the mobile region and the immobile region according to a number of parallel rates. While the MINC method was developed for fractured media, it may be used as general multi-continuum discretization approach, and in some MINC implementations multi-rate mass transfer models can be obtained as a specific case.


In the previous sections, we have reviewed approaches to simulate reactive transport processes in media characterized by the multiplicity of scales. These approaches have been brought to bear on several applications in the field. Here we presented a selection of applications where one or more of these approaches have been used. In some of these applications, to evaluate the ability of multiscale methods to capture the processes of interest, they are compared to micro-scale simulations, macro-scale simulations, or simulations with other multiscale approaches. We organize this section by first distinguishing between applications in granular porous media and fractured media. We move on to applications in integrated hydrology that connect surface and subsurface compartments, which share similarities to the multiscale concepts discussed in this chapter.

Granular porous media

Dissolution of rocks often involves the development of altered rinds in the grains that make up the rock, e.g., Navarre-Sitchler et al. (2009). These layers may be characterized by changes in porosity or diffusivity that determine the rate at which reactant accesses the reactive mineral. Rates of dissolution observed at the large scale depend on diffusion-reaction processes that take place at the scale of micrometers. Although the method of multiple interacting continua was developed for fractured media it can be viewed as similar to the “shrinking core” model (Wunderly et al. 1996). In this view, the grains can be conceptualized as consisting of several continua connected in series. Aradóttir et al. (2013) used this approach to capture the microscopic dynamics of basaltic glass dissolution in a Darcy-scale model (Fig. 10). The MINC method involves dividing the system up to ambient fluid and grains, using a specific surface area to describe the interface between the two. The various grains and regions within grains can then be described by dividing them into continua separated by dividing surfaces. Millions of grains can thus be considered within the method without the need to explicitly discretizing them. Four continua were used for describing a dissolving basaltic glass grain; the first one describes the ambient fluid around the grain, while the second, third and fourth continuum refer to a diffusive leached layer, the dissolving part of the grain and the inert part of the grain, respectively.

The physical heterogeneity of natural porous media often leads to poorly mixed conditions such that different conditions exists within relatively small volumes of the media. From a reactive transport perspective, the existence of micro-environments implies that geochemical gradients can develop. For example, diffusion-dominated micro-environments conditions may be at equilibrium or close to equilibrium with certain minerals, while in advection-dominated pore spaces reactions are far from equilibrium. Mineralogical heterogeneity at different scales compounds to this effect where a correlation may exist between a certain mineralogy and enhanced reactivity due to availability of micro-porosity (Landrot et al. 2012). Poorly mixed conditions may also exist within individual pores, which under certain conditions lead to scale dependent rates (Li et al. 2008).

Although pore-scale simulations make it possible to reproduce geochemical gradients at the micrometer scale (Molins et al. 2012), multiscale models offer an alternative that is less computationally demanding. Liu et al. (2015) designed a micromodel experiment coated with hematite where macro- and micro-porosity domains where present (Fig. 11a). Three separate models were used to describe transport processes and reductive dissolution of hematite: a pore-scale model, a 1D single-continuum model, and a 1D triple-continua model. The predictions from the pore-scale reactive transport model predicted reasonably well the measured pore-scale rates of hematite reduction. Geochemical gradients within the domain (Fig. 11a) made it necessary to divide the domain in three continua: one that captured advection-dominated domain, one to capture the diffusive gradients within the macro-pore and a third on to capture diffusive limitations in the micro-pores. While the rate of hematite reduction in the advection-dominated and macro-pore domains was affected by the flow rate, the rate in the micropore domain was not, however, as reactant diffusion was rate-limiting. Results from the single domain model deviated significantly from the pore-scale results.

Pore-scale is not always available due to the large dimensions of the domain and the associated large computation costs. In these circumstances, multiscale approaches that retain a pore-scale description for part of the domain while using Darcy-scale for the rest help bridge the trade-off between process resolution and domain size. Yan et al. (2017) used a Darcy–Brinkman–Stokes-based approach to simulate biogeochemical reaction rates in heterogeneous sediments. An X-ray computed tomography image was used to construct a 3D multiscale domain where e was assigned a value from grayscale image (Fig. 12a), and in turn, a permeability value calculated from e. A critical aspect of the model was, however, the assumption that the distribution of soil organic carbon (SOC) and biomass was correlated inversely with e. That is, these variables were high near or on solid surfaces while low in large pore spaces. Additional simulations with single- and dual-domain models were performed to test this assumption. These simulations show that only when a large fraction of the soil organic carbon and biomass was placed in the immobile domain, dual-domain models were able to capture effluent concentrations of nitrate (Fig. 12b). In fact, single-domain models captured well effluent concentrations for non-reactive tracers. The multiscale aspect of the problem appeared only in the reactive transport component.

These applications highlight the importance of mixing processes in reactive transport in heterogeneous porous media. In some instances, however, simulation of mixing processes in relatively homogeneous media may need of a multiscale approach when they are coupled to reactive processes. An example of this are mixing-controlled reactions. When two solutions mix such that a precipitate may form that has the potential to clog the pore space, it may be necessary to perform pore scale simulations. Scheibe et al. (2015b) presented a hybrid model that performed Darcy-scale simulations everywhere in the domain, and based on an incomplete mixing conditions, dynamically performed additional pore-scale simulations in a narrow region of the domain where precipitation occurred as a result of the mixing. This overlapping or hierarchical approach eliminated the need for matching boundary conditions between pore-scale and Darcy-scale domains. Hybrid simulations showed a sharper reaction front than equivalent Darcy-scale simulations, although some instabilities were observed in the hybrid approach (Fig. 13).

Fractured media

Flow and transport in fractured media occur primarily through a network of fractures, while flow in the matrix may be significantly slower with transport dominated by diffusive processes. While the fractures account for most of the flow and transport they typically make up a small portion of the overall volume of the medium. One could argue that to simulate fracture systems and incorporate this disparate scale, most fracture models has in one way or other multiscale aspects. Specific approaches to simulate fractures are reviewed in detail in a chapter of this volume (Deng and Spycher 2019, this volume). Here we describe the work of Molins et al. (2019) to develop a hybrid multiscale of fractured media as an example of the two separate scales. In this hybrid model, a pore-scale component captures Navier–Stokes flow, multi-component transport and aqueous equilibrium in the fracture, while a Darcy-scale component captures multi-component diffusive transport, aqueous equilibrium and mineral reactions in the porous matrix (Fig. 14). The interface between the sub-models, the fracture surface, is represented by an embedded-boundary. To simplify exchange of concentrations and fluxes at this interface, adaptive mesh refinement is used such that resolutions of the sub-models match at the interface while still using coarser resolution away from the interface when not needed in the Darcy-scale domain. The multiscale model is capable to capture flow channelization observed in an experimental fractured core and, at the same time, limitations in the dissolution of calcite by diffusive transport through an altered porous layer.

Surface–subsurface hydrologic coupling

Reactive transport of geochemical species in streams is result from an interplay between biogeochemical processes and mass exchange between the stream and the subsurface. The saturated sediment adjacent to the stream is therefore an important region for understanding the composition and evolution of water in the stream. For its role, the hyporheic zone has been the focus of study to understand flow and solute transport. Increasingly, there is interest to simulate reactive transport in the context of integrated surface–subsurface processes where both compartments are considered.

Although we have motivated the need for multiscale approaches in porous media from the pore- to Darcy-scale models in porous media, coupling of surface and subsurface processes requires the solution of similar equations in a coupled manner. The understanding of multiscale approaches in this sense is related to that of multi-physics, where the processes of interest are described by different equations. These processes may be characterized by different time scales, e.g., fast overland flow compared to long residence times for subsurface flow. Conceptually, these systems are similar to some of the subsurface systems considered in this chapter such as fractured media with fast flow in fractures compared to long residence times in the rock matrix. As a result, the approaches to coupling processes between the different compartments fall within those described in this review, including multi-rate approaches and hybrid approaches that require enforcement of continuity of mass and fluxes across interfaces. Examples of each of them are presented in what follows.

Painter (2018) use the residence time concept to develop a multiscale model for hyporheic exchange considering biogeochemical reactions. In this approach, the channel is a one-dimensional domain in which each cell of the discretization is connected to one dimensional sub-grid model for reactive transport (Fig. 15a), which are convolution representations of the exchange of solute with the hyporheic. This approach is mathematically equivalent to multi-rate mass transfer formulations such as Equations (3334). In Painter (2018), the sub-grid model is generalized to include multicomponent reactive transport with general nonlinear reactions. Hyporheic zone denitrification is simulated with these non-linear models to demonstrate the approach (Fig. 15b).

Hybrid approaches are also being brought to bear on the surface–subsurface hydrologic exchange. Bao et al. (2018) developed a one-way coupled surface and subsurface water flow model to simulate a 7-km long reach along the Columbia River (Fig. 16). While the subsurface model is a Darcy-scale model, the surface component is solved with computational fluid dynamics software that solves a form of the Navier–Stokes equations with a free surface boundary. The model was employed to investigate surface water fluid dynamics and the impact of subsurface structures on the hydrologic exchange.


The issue of scale is central in reactive transport modeling and the different disciplines it draws from. It is often implicitly behind common assumptions used in most models, for example to couple different processes at the appropriate scale. The term multiscale is often used to specifically refer to applications that explicitly consider the different scales present in a given application. Here we have motivated the use of multiscale approaches by the need to scale-up processes that take place at the pore scale, where the medium is treated as composed of separate fluid and solid phases, to the Darcy-scale continuum, where the porous medium is characterized as a continuum. As shown, formally upscaling reactive transport only results in separate macro- and micro-problems under a limiting set of conditions. In general, however the macro and micro problem are coupled.

Although insightful, mathematical methods for upscaling reactive transport rely on idealized representation of porous media and rarely used to investigate experimental or field systems. Simulation of these systems usually requires bringing to bear methods that use different scale descriptions within the same model such that micro-scale descriptions are used where it is needed; for example, to represent processes without making assumption of constitutive models for bulk parameters. Two approaches are prominent to simulate reactive transport in domains consisting of pore-scale and Darcy-scale sub-domains, one that uses a single equation to describe processes at both scales, and one that solves the problem in each domain separately and applies appropriate coupling conditions at the interface. These methods have also been used to find a compromise between the process and domain resolution and the size of domain that can be addressed with limited computational resources.

Multi-rate and multi-continua constitute another class of multiscale models that conceptualize the Darcy-scale porous medium as composed of two or more sub-regions that exists at each point in space and exchange mass according to a single-rate, multi-rate or continuum of rates. Mass exchanges between sub-regions or continua can represent micro-environments in which geochemical conditions may be different. Often these differences imply that a certain correlation exists between their hydrological accessibility and their mineralogical composition and reactivity.

Applications of these approaches in porous and fractured media show that often more than one approach can be used for a specific application. As application of reactive transport expands to consider surface processes, in addition to subsurface processes, available approaches are adapted and brought to bear on the surface–subsurface coupling. Generally multiscale approaches seek to represent processes at the appropriate scale, but different reasons are identified for this need. In some cases, the heterogeneous structure and composition of the porous medium necessitates micro-scale representation of the processes, while in other cases, the processes lead to the need for this micro-scale representation. For example, when solutions of different compositions mix, zones of active biogeochemical processes may require a micro-scale representation, especially when the media evolves as a result of it. While there is interest in capturing accurately these hot spots of biogeochemical activity, there is also interest in multiscale methods that can dynamically adjust process representation in time to capture periods of increased biogeochemical activity or hot moments.


The contribution of Sergi Molins to this chapter was supported by the Interoperable Design of Extreme-scale Application Software Watersheds (IDEAS Watersheds) project funded by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research under Award Number DE-AC02-05CH11231.

Open access: Article available to all readers online.