Mesoscale and Hybrid Models of Fluid Flow and Solute Transport

Fluid flow and reactive transport is relevant to many subsurface applications including CO₂ sequestration, miscible/immiscible displacements in enhanced oil recovery, wellbore acidization, pollutant transport, and leakage/remediation of nuclear waste repositories. In all these scenarios, one or more fluid phases flow through the complicated geometry of the pore space, while advecting one or more chemical species along their flow streamlines. Simultaneously, the chemical species undergo molecular diffusion, due to their Brownian motion, allowing them to randomly jump from one streamline to the next. In the case of fluid–fluid or fluid–mineral reactions, chemical species may be transformed, potentially leading to precipitation and/or dissolution of solid minerals that alter the geometry/topology of the pore space. This in turn affects the velocity field of flow, and thus transport via advection/diffusion. Such complicated feedback between these pore-scale processes could give rise to “emergent” manifestations at larger scales. These manifestations are referred to as “emergent” because they cannot be foreseen from the behavior of the individual pore-scale mechanisms involved. In order to make reliable predictions of flow and transport at any scale of interest, accurate models need to be developed.

Two spatial scales are commonly identified with a porous medium: the “micro/pore scale” (1–100 μm) and the “macro/continuum scale” (>1 m). The former is the fundamental scale in which physical processes (flow, transport, and geochemistry) take place, and the porous medium is regarded as discrete in nature (void space vs. grain space). The latter is a more practical scale, where we would ultimately like to have a reliable description of flow and reactive transport, and the porous medium is regarded as a continuum. The macroscopic parameters appearing in the description of continuum models, such as permeability or dispersion coefficient, are typically extracted from experiments or stand-alone pore-scale simulations. While such a “hierarchical” upscaling approach is often appropriate, there are many practical problems for which this is not the case. In these problems, perturbations of state variables (e.g., concentration) at the pore scale are tightly coupled to their averaged quantities at the continuum scale. In other words, a clear separation between the temporal scales and/or the spatial scales does not exist. An often encountered scenario is solute transport in the presence of strong fluid–mineral reactions (Kechagia et al. 2002; Battiato and Tartakovsky 2011). Any improvement to the continuum description of these processes requires a fundamental understanding of the relevant physics at the micro scale. In this effort, pore-scale models are an essential asset.

This chapter is divided into two parts. In the first, we review various pore-scale modeling approaches for simulating flow and reactive transport, and discuss their advantages and limitations. We focus primarily on single-phase flow and transport, and give a brief discussion on multiphase flow for completeness. Several review articles on the latter are provided as reference for the interested reader. Among the various pore-scale models, we place the greatest emphasis on pore-network models. Pore-network models are “intermediate-scale”, or “mesoscale”, models bridging the gap between the pore scale (1–100 μm) and the core scale (10 cm⁻¹ m). This is an important scale because it incorporates several thousand pores, where most macroscopic behaviors begin to emerge. Network models are also attractive because of their simplicity, which is why they can handle domain sizes much larger than other pore-scale models. However, their simplicity can also be the source of inaccuracy and ambiguity. In some cases, one can modify the simplifying assumptions of traditional network models (at the expense of increasing their complexity) and substantially increase their accuracy, while preserving computational efficiency. In other cases, such a modification may be inherently limited, and a different pore-scale approach has to be sought. We give examples of both scenarios in this chapter. In the second part, we turn our attention to the problem of transferring information obtained from pore-scale models to larger continuum models. This is an essential step, since the latter is the most practical scale of interest to most geochemists, hydrologists, and subsurface engineers. This chapter focuses only on a review of modeling strategies that are designed for scenarios in which the pore and macro scales are tightly coupled, and are not amenable to hierarchical upscaling. Instead these approaches seek to establish a dynamic “two-way communication” between the pore and the continuum scales, typically over a small region of the porous medium. In doing so, both pore-scale and continuum models are used within the same computational domain, which is why they are commonly referred to as “hybrid” models. The hope is that the following material familiarizes beginning researchers with current advances in pore-scale as well as hybrid modeling, and aids them in adopting and improving upon those methods that best fit their research interests.

Over the past few decades modeling fluid flow and solute transport at the pore scale has seen the adoption and development of various computational methods. The first and perhaps most critical step prior to any modeling effort is the accurate characterization of the pore-space geometry/topology. Imaging techniques such as X-ray microtomography (XMT) (Wildenschild and Sheppard 2013; Noiriel 2015, this volume) have made it possible to obtain accurate 3-D characterizations of the complex pore-space geometry of rock samples. These techniques have also been used to visualize and quantify pore structural changes as a result of biofilm growth (e.g., Iltis et al. 2011) and reactive transport (e.g., Noiriel et al. 2005) as shown in Figure 1. For granular media, Monte Carlo (Maier et al. 2003), cooperative rearrangement (Thane 2006), and sequential sedimentation (Coelho et al. 1997; Øren and Bakke 2002) algorithms are commonly used to digitally reconstruct representations of various grain packs, in which the pore space geometry is well-defined. Such reconstructions, although approximations of actual porous media, provide valuable insights into the link between depositional processes (e.g., cementation and compaction) and hydraulic/transport properties (e.g., permeability) of granular media (Bryant et al. 1993a,b; Bakke and Øren 1997).

Modeling can then proceed either by simulating directly on the complex void geometry, or on a simplification thereof. The former approach is typically referred to as direct modeling, while the latter is closely associated with pore-network modeling (see the Pore-network models section). In this section, we review a handful of direct modeling approaches recently used at the pore scale. These include computational fluid dynamics (CFD) (Wendt 2009; Molins et al. 2014; Trebotich et al. 2014; Molins 2015, this volume), Lattice-Boltzmann (LB) (Chen and Doolen 1998; Yoon et al. 2015, this volume), and smoothed particle hydrodynamics (SPH) (Monaghan 1992; Tartakovsky et al. 2008b). These methods can be used to solve the governing Equations (1–2) for single-phase, incompressible, Newtonian flow on the pore space.

Equations (1a) and (1b) represent the Navier–Stokes and continuity equations, respectively. v, p, ρ, μ, c denote fluid velocity, pressure, density, viscosity, and solute concentration. f_b denotes body forces exerted on the fluid. Equation (1c) is the no-slip boundary condition at the rock–fluid interface, Γ_fs. Since typical Reynolds numbers (Re = ρνd/μ; d is a characteristic length of the medium typically the grain size) in the subsurface are low (<1), the inertial term (i.e., v·∇v) is commonly ignored. Equation (2a) represents the species-balance equation with a homogeneous reaction rate of r_hm (D_m is the molecular diffusion coefficient), and Equation (2b) represents the heterogeneous reactive boundary condition (of rate r_ht) at the rock–fluid interface. While both fluid flow (i.e., Eqn. 1) and solute transport (i.e., Eqn. 2) can be solved via either of the foregoing direct modeling approaches, Eulerian (or mesh-based) methods (e.g., CFD) are typically preferred for computing flow and Lagrangian (or particle-based) methods (e.g., SPH) for transport, since they are devoid of numerical dispersion. Another popular Lagrangian method used for simulating transport at the pore scale is particle tracking (PT). In this method, particles (or random walkers) are propagated through the void space via a deterministic advection step that follows the streamlines, followed by a stochastic diffusion step that obeys Einstein’s equation; i.e., δr² = 2dD_mδt, where d is the dimensionality of the problem, and δr and δt are the stochastic spatial and temporal step sizes, respectively.

Maier et al. (2000, 2003) studied dispersion on digitally generated sphere packs using LB and PT to solve flow and transport, respectively, and reported good agreement with measurements from nuclear magnetic resonance (NMR) spectroscopy experiments. Molins et al. (2012) used an adaptive meshing CFD method to study the “lab-field discrepancy” of geochemical reaction rates (Fig. 2i). Their study provides valuable insight into the vital role of pore-scale modeling as a crucial guide towards predictive macro-scale modeling. Yang et al. (2013) used CFD to solve the Navier–Stokes flow equations on a bead pack and successfully compared the velocity field to that obtained from magnetic resonance velocimetry measurements. Zaretskiy et al. (2010) used a finite-element–finite-volume method to study longitudinal dispersion on a digitized sample of Fontainebleau sandstone. They highlighted the sensitivity of the pore-scale modeling results to the computational mesh employed. Mostaghimi et al. (2012) used a combination of finite difference and PT for flow and transport, respectively, on micro-CT images of Berea sandstone, and produced longitudinal dispersion coefficients in quantitative agreement with experimental data. Ovaysi and Piri (2011) used a Lagrangian approach referred to as the modified moving particle semi-implicit (MMPS) method to solve flow and transport on sandstone samples. They obtained longitudinal dispersion coefficients in favorable agreement with experimental data (Fig. 2iii), and included non-inertial effects in their simulations.

Kang et al. (2010) and Yoon et al. (2012) used LB to perform single-phase multi-component reactive transport simulations involving calcite precipitation and dissolution. Kang et al. (2010) additionally pointed out discrepancies between continuum and pore-scale predictions of the reactive front velocity. Multiphase flow and multicomponent reactive transport simulations using direct modeling approaches have also been performed in the literature (Parmigiani et al. 2011; Chen et al. 2013, 2014). For example, Parmigiani et al. (2011) used LB to simulate drainage accompanied by thermally induced melting of the solid phase (Fig. 2ii). Since the focus of this paper is primarily on single-phase transport phenomena, the reader is referred to Meakin and Tartakovsky (2009) for a complete formulation of the governing equations as well as modeling methods for multiphase flow and reactive transport at the pore scale.

Despite the high fidelity of direct modeling predictions and the fundamental insights they provide about the underlying physical mechanisms of a given phenomenon, they often demand high performance computational resources and massive parallelism (Oostrom et al. 2014). This poses a limit on the size and scale of the problems to which they can be applied. Such limitations have made pore-network modeling very popular over the past few decades, which overcome this limitation by simplifying the complex void space geometry while preserving essential features thereof. In the following section, we discuss pore-network models and provide examples of what we mean by “essential features”.

The popularity of pore-network models arose out of the pioneering works of Fatt (1956a,b,c), who studied two-phase drainage on a 2-D lattice network of randomly assigned throat radii. This was a radical shift from the hitherto bundle-of-tubes representation of porous media for computation of various macroscopic properties. Pore networks are simplified representations of the complex pore-space geometry and consist of an interconnected network of pores (or nodes) and throats (or bonds) (Fig. 3i). This geometric simplification is the reason why network models are more computationally efficient than direct pore-scale models, and applicable to much larger domain sizes. The elements of a pore network are typically assigned relatively simple shapes amenable to analytical treatment, e.g., spheres for pores and cylinders for throats. Scenarios in which either pores or throats are assigned zero volumes have also been considered. A summary of various pore/throat shapes used in the literature can be found in Joekar-Niasar (2010). The manner in which pores are connected to their nearest neighbors constitutes the network topology (coordination number is a parameter closely related to the network topology and is defined as the number of connected neighbors to a given pore) while the specific geometric idealizations used to represent pores/throats comprise the network geometry. Proper characterization of both the topological and geometric aspects of porous samples is the first necessary step, if pore-network models are to be predictive.

Early works on pore-network characterization typically involved statistical mappings of pore/throat properties (e.g., radii, coordination number) onto a lattice structure (Payatakes et al. 1980; Mohanty et al. 1987), or adjusting pore/throat properties to match one set of measurements (e.g., capillary pressure) followed by predictions of more difficult-to-measure properties (e.g., relative permeability) (Fischer and Celia 1999; Vogel 2000). However, statistically mapped networks often ignore spatial correlations present in real rocks due to the random assignment of pore/throat properties, and adjusted networks (to match one set of measurements) suffer from non-uniqueness in their representation. For these reasons, statistically generated networks are typically not considered to be predictive. Bryant and coworkers (Bryant and Blunt 1992; Bryant et al. 1993a,b) introduced the concept of physically representative networks, which marked the beginning of truly predictive network modeling. In their work, permeability, relative permeability, capillary pressure, and permeability–porosity relationships, were all accurately and directly predicted without adjustable parameters using only grain positions in a disordered sphere pack (measured by Finney 1968). The network was extracted via Delaunay Tessellation of the sphere centers (Fig. 3i), and its main attraction was that spatial correlations of the pore space were naturally imbedded into the extracted network. Recent advances in computer and imaging facilities has given rise to various image analysis techniques whereby physically representative networks can be obtained from digitized images of real samples. These include medial-axis (Thovert et al. 1993; Lindquist et al. 1996; Lindquist and Venkatarangan 1999; Prodanovič et al. 2006), watershed (Sheppard et al. 2006; Thompson et al. 2008), and maximalball (Silin and Patzek 2006; Dong and Blunt 2009) algorithms. These techniques have resulted in further quantitative predictions of single-phase and multi-phase flow properties in water-wet and mixed-wet media (e.g., Bakke and Øren 1997; Øren et al. 1998; Patzek 2001).

Since Fatt’s (1950a,b,c) works, pore-network models have been used in many other applications including: non-Newtonian flow (Lopez et al. 2003; Balhoff 2005), non-Darcy flow (Thauvin and Mohanty 1998; Balhoff and Wheeler 2009), solute dispersion (Sahimi 1986; Sorbie and Clifford 1991; Bijeljic et al. 2004, Bijeljic and Blunt 2007; Acharya et al. 2007b), reactive transport (Hoefner and Fogler 1988; Li et al. 2006; Algive et al. 2010; Kim et al. 2011) (Fig. 3iii), multi-phase flow (Koplik and Lasseter 1985; Al-Gharbi and Blunt 2005; Piri and Blunt 2005; Joekar-Niasar et al. 2010; Ellis et al. 2011) (Fig. 3iv), biofilm growth (Suchomel et al. 1998a,b), and modeling dual-scale media containing micron- and submicron-scale porosity (Mehmani et al. 2013; Mehmani and Prodanovič 2014a,b; Prodanovič et al. 2014) (Fig. 3ii). A larger portion of this literature is devoted to two-phase flow (drainage and imbibition), due to its significance in petroleum engineering (oil recovery) and soil sciences. The reader is referred to comprehensive reviews of pore-scale modeling of multiphase flow by Celia et al. (1995) and Blunt (2001) for more information. An excellent review by Berkowitz and Ewing (1998) further provides insights into the close connection between (quasi-static) two-phase flow and (invasion) percolation theory. It is noteworthy, that while the development of multiphase flow network models as well as single-phase reactive/passive transport network models have separately seen considerable progress over the past few decades, a merger between the two is rarely observed in the literature. Since the geologic storage of CO₂, for example, involves the flow of supercritical CO₂ and brine in the presence of multicomponent geochemical reactions with the potential for mineral precipitation/dissolution, such a merger is very much needed.

Pore-network models are often regarded as “mesoscale” models acting as a critical bridge between our fundamental understanding of fluid flow/transport physics at the pore scale (1–100 μm) and our practical desire for modeling at the continuum scale (>1 m). The reliability of such intermediate-scale models, however, hinges upon their underlying simplifying assumptions. There are two levels of approximations made in every pore-network model. The first is geometric, and occurs in replacing the complex void space geometry with an “equivalent” pore network. The other is physical, and occurs when we make simplifying assumptions to describe the pore-scale physics on this idealized geometry. The key to constructing predictive network models is to identify and capture the most essential geometric and physical features associated with a given phenomenon (as alluded to in Direct pore-scale modeling). These essential features are not invariant from one problem to the next. In multiphase flow, for example, it is very important to include angular geometries in the description of the pore network. Excluding this feature ignores the possibility of corner flow of the wetting phase, and snap-off of the non-wetting phase because of it. On the other hand, including inertial effects in the description of two-phase flow in a pore-network model seems to be of little importance if residual saturation trapped behind an advancing drainage front is to be predicted (Moebius and Or 2014). In single-phase flow, angularities are unimportant and only accurate values for hydraulic conductivities of the network elements are required. When geochemical reactions are present, neglecting to accurately quantify fluid–mineral interfacial areas will lead to large errors in macroscopic predictions. In passive solute transport, while assuming perfect mixing conditions within individual pores is inappropriate for ordered media, it appears to be a good assumption for disordered media (see Network modeling of solute transport).

In identifying these essential geometric and physical features, direct pore-scale modeling as well as controlled (e.g., micromodel) experiments play an indispensable role. Once an essential feature is identified, traditional network models may be modified to attain improved predictive accuracy. Although such a modification may not always be possible, in which case a different mesoscale modeling strategy has to be sought. In the following section, we review network modeling of single-phase solute transport and provide two examples where direct pore-scale modeling was used to establish the need for modifying traditional mesoscale assumptions. In the first, modifications lead to substantial improvements in the predictive capacity of the traditional network model, while in the second improvements appear to be inherently limited.

Pore-network modeling of solute transport has received special interest from many authors in the past few decades and several methodologies have been proposed. A prerequisite to simulating transport is the computation of the velocity field within throats. The procedure is quite standard and involves: a) imposing pressure boundary conditions on the pore network, b) describing flow rates within throats via a constitutive equation e.g., Hagen–Poiseuille for Newtonian fluid in a cylindrical throat, c) enforcing mass balance at each pore e.g., Equation (3), where g_ij is the throat hydraulic conductivity, Δp_ij the pressure drop across the throat, and N_i^th the number of throats connected to pore i, d) solving the resultant system of (linear or nonlinear, depending on fluid rheology and/or flow regime) equations for pore pressures, and e) computing throat flow rates/velocities from the aforementioned constitutive throat equation.

The simulation of transport then resumes using computed throat flow rates/velocities. Below in this section, we review several Eulerian and Lagrangian approaches developed in the literature for modeling transport on pore networks. Advantages and limitations of each method are outlined and discussed. We then focus on the main sources of error and difficulty within each modeling class, and the possibility of minimizing them. Finally, we provide a brief review of modeling homogeneous/heterogeneous reactions on pore networks.

Bryntesson (2002), Acharya et al. (2005, 2007b), Li et al. (2006), Kim et al. (2011), Mehmani et al. (2012), and Nogues et al. (2013) are among those who have adopted the popular mixed-cell method (MCM), in which solute-balance equations are written for each pore, i.e., Equation (4).

In Equation (4)c_i and V_pi are the concentration and volume of pore i; q_ij, a_ij, and l_ij are the throat flow rate, cross-sectional area, and length, respectively; R(c_i) is the reaction term. In this method, throats are assumed to be volumeless and the solute within the pores perfectly mixed (i.e., single concentration value assigned to each pore). Solute flow rates within throats are formulated as the algebraic sum of an upwinded advection term (first two terms on the RHS of Eqn. 4) and a linearly varying diffusion term (third term on the RHS of Eqn. 4). Essentially, MCM can be regarded as a low-order, finite-volume method on an unstructured grid, i.e., the pore network. The advantage of MCM is that it is very computationally efficient and highly adaptable to various transport scenarios. For example, Acharya et al. (2005) used MCM to study non-linearly adsorbing solute transport, and determined that more than a million pores were required for their results to be statistically representative. Li et al. (2006) and Kim et al. (2011) studied complex geochemical reaction kinetics of anorthite and kaolinite precipitation/ dissolution relevant to CO₂ sequestration. Nogues et al. (2013) studied porosity/permeability evolutions in carbonates due to carbonic-acid-driven precipitation/dissolution reactions. They considered 18 aqueous species and 5 mineral species undergoing 14 independent reactions. The flexibility and computational efficiency of MCM is why such complex systems acting on sufficiently large pore-scale domains can even be considered.

A number of variants and/or modifications of MCM have also been developed in the literature. For example, Raoof et al. (2013) assign volume to both pores and throats (solute still perfectly mixed in both) and sub-discretize the wetting filaments in the corners of partially drained pores to account for the partial mixing of solute within them (Fig. 4iv). Van Milligen and Bons (2014) proposed a modification to the throat rate expressions used in MCM (i.e., first three terms on the RHS of Eqn. 4) by deriving analytical expressions based on a steady-state plug-flow assumption within throats. They occasionally further sub-discretized throats into smaller “pores” for increased modeling resolution. A generalization of this method for non-uniform velocity profiles was given by Mehmani (2014) (referred to as the rate-modified MCM), but it was shown that these modifications offered little improvement in model predictions compared to MCM. Algive et al. (2010, 2012) and Varloteaux et al. (2013) similarly modified throat solute rate expressions in the context of reactive transport. These works employed moment theory to derive corrected macroscopic parameters (i.e., solute mean velocity, dispersion coefficient, reaction source term) for each pore/throat element in the long-time asymptotic regime. The model was used to study the effects of dissolution/ precipitation reactions on macroscopic single- and two-phase flow properties in the contexts of CO₂ sequestration and diagenetic alterations in carbonate rocks. Suchomel et al. (1998b) developed a model in which pores were assumed volumeless and throats were sub-discretized into finite-difference grids. Interpore diffusion (equivalent of the third term in the RHS of Eqn. 4) was implicitly incorporated by adjusting numerical diffusion via grid and time-step sizes (although this does not account for diffusion countercurrent to the flow direction), and perfect mixing was assumed within pores. The model was used to study permeability/porosity alterations during bacterial biofilm growth in porous media. A simple but interesting model was developed by Martins et al. (2009), in which the solute balance equations at the pores were formulated as a system of delay-differential equations, i.e., Equation (5).

In Equation (5), τ_ij is the throat residence time, which accounts for the delay in transport from one pore to the next. In essence, the model uses upstream concentrations from earlier time steps to compute solute flow rates flowing into pores at later times. This is implemented by sub-discretizing the throats and marching pore concentrations forward within the sub-discretized segments (akin to a traveling wave). However, several limiting assumptions were made including the neglect of diffusion, plug-flow within throats, and perfect mixing within pores.

A different set of models formulate transport equations in the Laplace transform domain with respect to the time variable. There is a certain extent of elegance and convenience associated with working in the Laplace domain, mainly due to the fact that convolutions of transit-time probabilities are converted to simple algebraic multiplications. Another advantage is that computation of temporal moments becomes rather straightforward (the Laplace transform of transit-time distributions act as moment generating functions). De Arcangelis et al. (1986) developed first-passage-time probabilities for tracer particles that move through throats connecting two neighboring pores. They assumed perfect mixing at the (volumeless) pores and plug flow with molecular diffusion at the throats. Under these conditions, they derived exact transit probabilities in the Laplace domain for particle motions in a network. They then used a “probability propagation” algorithm to determine the first-passage-time distribution of a 10 × 10 diamond lattice network and computed longitudinal dispersion coefficients for various Péclet numbers (= advection/diffusion). At no point was the time-domain concentration field computed. Koplik et al. (1988) used the same set of equations as de Arcangelis et al. (1986) but diverged in their analysis by writing species balance equations for each pore in the Laplace domain. The resulting linear system of equations was then solved and numerically inverted into the time domain using the Stehfest (1970) algorithm. A strategy for computing higher-order moments of first-passage-time distributions of the network was further outlined. This method was later extended by Alvarado et al. (1997) to reversible adsorption/reaction scenarios, where they arrived at the interesting conclusion that dispersion coefficients depend on the degree of spatial heterogeneity of reactive sites in a porous sample and scale non-linearly with Péclet number. However, they noted that the numerical Laplace inversion step was prohibitive for networks larger than 20 × 20 pores and inaccurate for large Péclet numbers (>10). Indeed numerical inversion of the Laplace transform is known to be notoriously difficult (often unstable) and ill-posed in computational and applied mathematics. For this reason, although valuable for performing moment analyses, time-domain predictions via these methods on representative sample sizes seem impractical and unlikely.

A more natural description of transport is provided by Lagrangian models, among which particle tracking (PT) is almost exclusively employed. In this method the steady-state flow equation is first solved on the network, as described in the beginning of this section, to obtain mean fluid velocities within each throat. Then, depending on the specific throat geometry, a (rectilinear) velocity profile is assumed and used to track particles from pore to pore, subject to simultaneous convection and molecular diffusion (Fig. 4i). Lagrangian network models are generally more computationally expensive than Eulerian network models, although more accurate. The limiting factor is the large number of solute particles often required in Lagrangian network models to obtain statistically converged results. PT methods on pore networks can be divided into two categories: a) those that trace particle motions in detail within throats following a discrete-time random walk process (Bruderer and Bernabé 2001; Bijeljic et al. 2004; Acharya et al. 2007a; Jha et al. 2011), and b) those that perform continuous-time random hops from one pore to the next (without explicit throat-level simulations) using throat transit-time distributions (Sahimi et al. 1986; Sorbie and Clifford 1991; Rhodes and Blunt 2006; Bijeljic and Blunt 2006; Picard and Frey 2007). We refer to the first class as DPT (discrete-time particle tracking) and to the second as CPT (continuous-time particle tracking).

Compared to CPT, DPT simulations are more time consuming since computational performance is limited by the time step size (controlled by the minimum throat residence time throughout the network). However, DPT can be quite accurate and has been used to successfully predict dispersion coefficients in unconsolidated granular media (e.g., Bijeljic et al. 2004; Jha et al. 2011) (Fig. 4ii). It is also very flexible in the sense that various velocity profiles (e.g., parabolic or plug-flow) within the throats can be considered, and one has substantial control on how to reassign particles to new outlet throats upon their arrival at the pores. Specifically, particles can be reassigned to any desired throat connected to the arrival pore, and even any desired location on the cross-section of that throat. In contrast, CPT is computationally more efficient but comparatively less flexible and less accurate especially in ordered media, as discussed in Difficulties and Sources of Error. The efficiency is due to the fact that particle motions within throats are not explicitly simulated, but are instead imbedded in the throat transit-time distributions. The reduced flexibility and accuracy are due to the loss of control in reassigning incoming particles to accurate cross-sectional locations of the outlet throats of a pore. This means that once an outlet throat is chosen (based on some probability), a transit-time distribution is used to determine the time the particle requires to exit the throat. But the transit-time distribution is typically stationary in time and is derived with uniform inlet conditions for concentration. Therefore, there is no way to specify the cross-sectional location an incoming particle should be assigned to in an outlet throat of a pore. This is particularly important in simulating dispersion in ordered media, for which CPT will not yield the expected (e.g., Edwards et al. 1991) D_L ~ Pe₂^d scaling, where D_L is the longitudinal dispersion coefficient, and Pe_d is the Péclet number defined as the ratio of advection to diffusion (Mehmani 2014).

CPT methods can be further subdivided into those that use deterministic transit-time distributions (e.g., Sorbie and Clifford 1991; Picard and Frey 2007), and those that use ensemble-averaged transit-time distributions (e.g., Bijeljic and Blunt 2006; Rhodes et al. 2009). Bijeljic and Blunt (2006) derived an ensemble-averaged transit-time distribution for Berea sandstone by fitting a truncated power-law distribution to their simulation results. They provided physically meaningful interpretation of the distribution variables and fitted their simulated data with a single adjustable parameter. On the other hand, deterministic transit-time distributions are often derived based on similar mathematics and assumptions as those already discussed in the context of Laplace-domain Eulerian network models (e.g., Rhodes and Blunt 2006; Picard and Frey 2007). These assumptions include perfect mixing within volumeless pores and plug flow within throats. Deterministic transit-time distributions are often numerically inverted from the Laplace domain to the time domain in order to draw particle transit times within the throats. An exception to this is Sorbie and Clifford (1991) who derived deterministic transit-time distributions based on rigorous single throat simulations assuming non-uniform velocity profiles.

The most common sources of ambiguity and error in both Eulerian and Lagrangian network models of single-phase solute transport are in the descriptions of partial mixing conditions within pores (Fig. 5a), and shear dispersion (i.e., spreading of the solute due to non-uniform velocity profiles) within throats (Fig. 5b). Accurate modeling of these fundamental transport physics could have a significant impact on quantitative macroscopic predictions of solute dispersion and effective reaction rates. Almost all Eulerian network models developed in the literature seem to assume perfect mixing within pores and neglect shear dispersion within throats. Models assuming Taylor–Aris shear dispersion coefficients within throats are unrealistic because throat lengths are typically not long enough for an asymptotic regime to be reached (Mehmani 2014). On the other hand, the incorporation of shear dispersion into Lagrangian network models (e.g., PT) is quite straightforward (see the discussion of these models above), but accurate description of pore-level partial mixing is still a source of ambiguity in these models.

In reality, perfect mixing occurs when diffusive forces are much more dominant than either advective or reactive forces. In the absence of reactions, increasing the Péclet number changes pore-level mixing conditions from perfect mixing to partial mixing. Figure 5a shows direct numerical simulations of this scenario in a 2-D pore (Mehmani et al. 2014). When reactions are present, increasing the Damköhler number (= reaction/diffusion), while keeping the Péclet number constant, causes a high-concentration-gradient boundary layer to develop near the fluid–solid interface. In pore-network models, the effect of reaction-induced boundary layers may be implicitly taken into account in the calculation of pore concentrations, although this scenario unfortunately has not been explored enough in the literature. In the following, we shall focus on shear dispersion and partial mixing in the absence of chemical reactions, and defer further discussion on the effects and modeling approaches of reactions until the very end of this section.

The ambiguity in accurately describing partial mixing within pores stems from the difficulty in approximating flow streamlines within pores and the extent of diffusive mixing occurring therein (Fig. 5a). Sahimi et al. (1986) and Bruderer and Bernabé (2001) developed simple intuitive rules for the redistribution of solute particles from the inlet to the outlet throats of a 2-D cross-shaped volumeless pore. Both ignored diffusive mixing within pores, in the sense that the particles could not randomly “hop” between outlet throats. Sorbie and Clifford (1991) introduced general heuristic, and admittedly approximate, rules for redistributing particles arriving at a pore among its outlet throats. These rules, and variants thereof, were subsequently applied and analyzed in later publications (e.g., Acharya et al. 2004, 2007a,b; Bijeljic et al. 2004; Bijeljic and Blunt 2007). At high Péclet numbers, these rules reduce to redistribution of solute particles based on flow-rate-averaged probabilities; at low Péclet numbers, they reduce to redistribution based on throat cross-sectional-area-averaged (multiplied by other corrective parameters) probabilities. These rules, however, are quite inaccurate at moderate-to-high Péclet regimes as they completely ignore the underlying streamlines that map the inlets to the outlets of a given pore (Mehmani et al. 2014). Jha et al. (2011) proposed a deterministic mapping of incoming particles to the outlet throats, intendedly along the streamlines, in an attempt to generalize the rules developed by Bruderer and Bernabé (2001). The approach was limited to pores with a coordination number of four or less, and ignored mixing due to molecular diffusion. Despite being attractive due to its simplicity, the mapping neglected spatial orientations of the inlet throats, and consequently the interaction between inflowing streams, and was therefore less predictive (Mehmani et al. 2014).

The mixing problem within pores has additionally received quite a bit of attention from the field of fracture-network modeling (e.g., Berkowitz et al. 1994, Park and Lee 1999; Park et al. 2001a,b; Johnson et al. 2006). Although in this context, the “pores” correspond to the juncture at which two fractures intersect. Park and Lee (1999) derived physically sound, efficient, and accurate transition probabilities for mapping particles from the inlets to the outlets of fracture junctions. In doing so, they conceptualized partial mixing at fracture junctions as a one-dimensional problem amenable to analytical treatment (Fig. 4iii). While accurate for fracture networks, the mixing equations of Park and Lee (1999) are not applicable to 3-D networks of porous media. This is because fracture junctions are essentially 2-D cross-shaped “pores”, which are geometrically simpler than the 3-D pores with non-planar throat orientations found in porous media.

In an effort to accurately capture partial mixing within pores and shear dispersion within throats, Mehmani et al. (2014) and Mehmani (2014) developed two pore-network models, respectively: the streamline splitting method (SSM) and superposing transport method (STM). These models were developed with the aim of exploring the possibility of capturing the pore-scale physics discussed above under an Eulerian framework. The insistence on the models being Eulerian was because, in the context of pore-network modeling, such a framework is generally computationally more efficient than a Lagrangian counterpart, albeit possibly less accurate. Therefore, the goal was to find an appropriate balance between acceptable predictive accuracy and computational efficiency. To satisfy the latter, both models used MCM as the starting point for their development. In MCM, a single concentration value is assigned to every pore, which implicitly assumes perfect mixing within them, independent of the local Péclet regime. Species-balance equations are then written in every pore, and solved for the pore concentrations. To circumvent the perfect mixing assumption, SSM divides the volume of the pores into smaller compartments (or pockets) and assigns different concentration values to each compartment. The number of compartments within a pore is equal to the number of inlets of that pore, as each inlet may carry with it fluids with different concentrations that need to be kept separate. However, mass transfer (or mixing) between two adjacent compartments within the same pore is allowed, which takes place purely due to molecular diffusion (perpendicular to the streamlines). Unlike Park and Lee (1999), pore-wall effects on the mixing process were also taken into account. The outlets of each compartment are determined by “splitting” the inflowing “streams” between the outlet throats, using a constrained optimization algorithm; which was shown to be also applicable in Lagrangian network models (Mehmani et al. 2014). The compartments are then regarded as control volumes and species balance is enforced in all of them. Figure 6 provides a conceptual picture underlying SSM and contrasts it to that of MCM. It was shown that SSM predictions were in very good agreement with direct CFD simulations as well as micromodel experiments, and computational costs were only marginally higher than MCM.

On the other hand, the focus in STM is placed on shear dispersion due to the non-uniform velocity profiles within throats. This important pore-scale process is ignored in MCM as it assumes throats to be volumeless, which causes the solute to instantaneously transport from one pore to the next. In STM, elementary species rate expressions are constructed for each throat, as a function of their aspect ratio (i.e., diameter over length). These elementary rate expressions are then used to perform space-time superpositions across the pore network. This is accomplished by recording pore concentrations as they dynamically evolve during a simulation, which makes it similar to, but much more general than, the model developed by Martins et al. (2009) (the latter ignores shear dispersion and molecular diffusion). The superpositions in STM have the effect of performing convolutions of the elementary rate expressions across the pore network. This makes STM useful in other fields such as signal transmission in electrical engineering. The obvious limitation of STM is that it is applicable to linear transport problems only. This excludes scenarios with non-linear fluid–mineral reactions such as those encountered in CO₂ sequestration. But, for linear problems, the method was verified against direct CFD simulations and simple convolution integral expressions. It was additionally validated against longitudinal dispersion experiments for unconsolidated bead packs from the literature.

Yet a final question remains to be answered: what is the actual magnitude of the improvements imparted to the prediction of a parameter of interest by incorporating these additional pore-scale details into our network models? After all, modeling efforts should ideally seek ultimate simplicity while disposing of unnecessary complexity. Mehmani et al. (2014) showed that the impact of pore-level mixing assumptions, partial vs. perfect, is very significant in ordered media (e.g., micromodels), but comparatively insubstantial in 3-D disordered granular media (e.g., sandstones). In fact, the average difference between the concentration fields obtained via SSM and the relatively less accurate MCM for disordered granular media was ~6% of the maximum concentration value (although the impact on upscaled transverse dispersion coefficients is yet to be studied). A similar conclusion was drawn by Park et al. (2001a,b), who used PT to study the impact of various mixing assumptions (perfect mixing vs. no mixing) at fracture junctions on overall macroscopic transport behavior in ordered and disordered (2-D) fracture networks. They found that mixing assumptions have a larger impact in ordered media compared to random media (i.e., less than 5% of fracture junctions were affected in random media), and they attributed this to the lower effective coordination number and higher inlet flux ratios, at fracture junctions, in random networks. Despite the various differences between the two media (porous vs. fractured), the conclusions appear to be the same: for disordered (granular/fractured) media the model with the simpler description of pore-level mixing is the one that is preferred (e.g., MCM).

On the other hand, Mehmani (2014) showed that the incorporation of shear dispersion within throats in network models of disordered granular media is rather important at high Péclet numbers (>100). Namely, longitudinal dispersion coefficients predicted by STM undergo a ~2.5-fold improvement compared to MCM. This result corroborates earlier works that used PT to simulate longitudinal dispersion (Acharya et al. 2007a; Jha et al. 2011). A more important conclusion drawn using STM was that all Eulerian network models, including STM, are inherently limited at sufficiently high Péclet numbers when applied to ordered media (e.g., micromodels). This was demonstrated by simply conceptualizing a straight tube under purely advective transport, as an equivalent 1D “network” (or string) of shorter tubes (Fig. 7a). The concentration profiles obtained from STM and MCM were shown to be representative of neither of the analytically known flux-averaged or cross-sectional-averaged concentration profiles (Fig. 7b). Even worse was the fact that simulated profiles were shown to converge towards a Gaussian distribution with an increase in the number of tube segments traveled. This is known to be impossible under a purely advective regime (Fig. 7c). It is known from Taylor–Aris theory that in the presence of molecular diffusion (no matter how small its contribution) a D_L ~ Pe_d² asymptotic scaling is expected from this system (D_L is the longitudinal dispersion coefficient, and Pe_d is the Péclet number); similar to ordered media (Edwards et al. 1991). Instead, an Eulerian network model (e.g., STM) would yield a scaling that is closer to D_L ~ Pe_d at high Pe_d, since Gaussian profiles are reached not because of molecular diffusion, but because of cross-sectional smearing of concentrations at the “pores” (denoted by dashed lines in Fig. 7a). In other words, solute particle “memories” are effectively erased upon their arrival at the pores (i.e., Markovian process), whereas in reality they are retained over a given distance that depends on the local Péclet regime. This distance is larger than the length of a single throat at sufficiently high Péclet numbers. The foregoing problem is not just specific to Eulerian network models. Lagrangian methods that draw throat transit times from the same distribution, that is stationary in time, are equally limited. An example includes the CPT model developed by Sorbie and Clifford (1991). The implication is that the “simplest” model, capable of accurately capturing shear dispersion in ordered media, must preserve particle memories in passing them from the inlet to the outlet throats of a pore. This means that it is not just sufficient to know how long a particle takes to transit a given throat, but also the cross-sectional position of that throat whence it exits. Such detailed information seems only resolvable in a pore network, with adequate simplicity, under a Lagrangian framework.

Comprehensive reviews on modeling reactive transport can be found in Steefel and MacQuarrie (1996), Steefel et al. (2005, 2013). Reactions at the pore-scale can be divided into those occurring within the fluid bulk, referred to as homogeneous reactions, and those that take place at the fluid-mineral interface, referred to as heterogeneous reactions. Homogeneous reactions themselves consist of zero-order decay processes independent of species concentrations, and higher-order concentration-dependent reactions that occur due to the mixing of two or more solute species. Homogeneous reactions are commonly modeled as instantaneous processes at equilibrium, while heterogeneous reactions as kinetically controlled processes. In direct pore-scale models, the former is described via source terms in the solute balance equation (Eqn. 2a), and the latter via appropriate boundary conditions at the fluid–mineral interface (Eqn. 2b). In Eulerian network models, where pores are viewed as continuously stirred tank reactors, all reaction types (as well as adsorption) are almost invariably described via source terms in the mesoscale solute balance equation (Acharya et al. 2005, 2007b; Li et al. 2006, 2007a,b; Kim et al. 2011, 2012; Mehmani et al. 2012; Raoof et al. 2013; Nogues et al. 2013). This effectively ignores concentration gradients and transport-limited effects on reaction rates at the scale of individual pores. Depending on the application, this may or may not be a good assumption. Li et al. (2008) showed that for typical flow and chemical/mineralogical conditions relevant to geologic CO₂ sequestration, pore-level perfect mixing appears, for all practical purposes, appropriate. Their conclusion was based on molecular diffusion coefficient and activity of the hydrogen ion among other assumptions. In problems where these parameters assume very different values, transport-limited effects can become important and the perfect mixing assumption inappropriate.

The network models developed by Algive et al. (2010, 2012) and Varloteaux et al. (2013) are the only known exception, in which concentration gradients due to fluid–mineral reactions are implicitly taken into account by calculating effective transport parameters for the pores and throats. However, their model is limited to long-time asymptotic regimes, where changes in the concentration field are slow with time. Moreover, none of the above mentioned models accounts for diffusion-limited mixing reactions that occur when two fluids with differing chemical compositions come into contact within a pore. In SSM, discussed in Difficulties and sources of error, the rate of mass transfer between two adjacent compartments within a pore (pk1 and pk2 in Fig. 6b) is computed via solving a local bounded Riemann problem. Mehmani et al. (2014) proposed the possibility of modifying this local problem (via appropriate initial/boundary conditions) to derive mass transfer rates in the presence of mixing-induced reactions. For reactions at the fluid–solid interface within the throats, one may derive modified elementary rate expressions that can be used in the STM model discussed in Difficulties and sources of error (although limitations for ordered media still persist). This would complement the long-time asymptotic method of Algive and coworkers, for short-time pre-asymptotic regimes. Finally, particle tracking (PT) can be used to describe all foregoing reactions types at the pore scale in high detail. This requires the development of non-trivial stochastic particle interaction criteria and collision probabilities between fluid–fluid and solid–fluid species (Gillespie 1977; Benson and Meerschaert 2008; Dentz et al. 2011). However, the adoption and implementation of PT in pore networks for accurate simulations of geochemical reactions have unfortunately been scarce.

Flow and reactive transport phenomena in the subsurface often span a wide range of spatial scales (nanometer to kilometer), which renders the development of predictive models capable of accurately bridging them a formidable task. Macroscopic parameters (e.g., permeability) or closure relations (e.g., capillary pressure vs. saturation) are commonly extracted from pore-scale modeling or experiments on representative samples of the real medium, followed by their direct implementation into continuum-scale simulators. While such an approach is appropriate for situations in which a clear separation between scales exists, it may lead to large errors otherwise (Kechagia et al. 2002; Tartakovsky et al. 2008b; Battiato and Tartakovsky 2011). Mixing induced precipitation ensuing form a bimolecular reaction is one such example (Battiato et al. 2009). In these cases a dynamic communication between the pore and the continuum scales seems to be required.

Li et al. (2006, 2007a,b) used a pore-network model to study reaction kinetics of kaolinite and anorthite in the context of geologic carbon sequestration. They demonstrated that reaction rates obtained from continuum-scale representations of transport and/or using volume/flux-averaged concentrations in reaction rate expressions leads to large errors, sometimes even wrongly predicting the direction of the reactions (i.e., precipitation vs. dissolution). Although their study was qualitative in nature, using 3-D regular lattice networks, it provided an explanation for the commonly reported discrepancy between reaction rates observed at the field scale and those obtained from well-mixed batch experiments on crushed samples. Namely, transport limitations at the pore scale control overall reaction rates, which are non-existent under well-mixed conditions in batch experiments. Effects of flow rate and reactive cluster size/ abundance were also studied, and it was concluded that the higher the degree of incomplete mixing (i.e., spatial variability of concentration) the higher the scaling error. Incomplete mixing was found to be strongest at medium flow rates (or Péclet numbers). Kim et al. (2011) and Kim and Lindquist (2012) extended the work of Li and coworkers using networks extracted from XMT images of real sandstones. They were able to determine surface mineral distributions from the images allowing for better quantitative analysis. Similar conclusions were drawn regarding the “lab–field discrepancy”, and an approximately power-law scaling of reaction rates vs. flow rate was reported for anorthite, while a more complex scaling emerged for kaolinite. In the context of filtration combustion in porous media, Lu and Yortsos (2005) similarly observed that spatially averaged macroscopic reaction rates were very different (discrepancies of a factor of two or higher) than those determined from using averaged variables in microscopic reaction rate expressions. They attributed this to the strong influence of microscopic heterogeneities on macro scale behavior. Recently, Molins et al. (2012) conducted sophisticated direct pore-scale simulations of calcite dissolution, and showed that pore-scale heterogeneities can result in an underestimation of reaction rates, due to mass transport limitations, even when total reactive surface area and porosity are held constant between samples.

Kechagia et al. (2002) demonstrated that for reactive transport with fast/finite kinetics, homogenization of microscopic equations via volume averaging does not hold except at the limit of macroscopic equilibrium. They further showed that even under these circumstances an eigenvalue problem enforces a coupling between the micro and the macro scales. A systematic study was performed by Battiato and Tartakovsky (2011) to identify transport regimes (characterized by the Damköhler (Da) and Péclet (Pe) numbers), in which a purely continuum description of advection-diffusion with nonlinear surface reactions breaks down. They used multiple-scale expansions to upscale the pore-scale equations, and presented their results in the form of a Pe–Da phase diagram. The results were in agreement with an earlier work that used volume averaging (Battiato et al. 2009), which substantiated their independence from the specific upscaling method employed. A recent work by Boso and Battiato (2013) extended this analysis to three-component systems undergoing two homogeneous and one heterogeneous reactions (all reversible). It is interesting to note, that Molins et al. (2012) located their simulations on the aforementioned Pe–Da phase diagram and determined their correspondence to a case in which spatial scales were coupled. While all foregoing examples focused on geochemical reactions, a coupling between scales is also observed in problems involving viscous and density-driven instabilities in multiphase and miscible displacement among others (Tartakovsky et al. 2008c).

The implications of the foregoing studies (among others) have given rise to a new class of modeling approaches referred to as “hybrid multiscale methods” (Scheibe et al. 2007), in which micro- and macro-scale simulations are simultaneously performed on the same computational domain. Balhoff et al. (2007) were one of the first to couple a pore-scale model to the continuum. However, their iterative coupling strategy had strong limitations in terms of flexibility and efficiency. These limitations were later lifted by Balhoff et al. (2008) through the introduction of mortars. Mortars are finite-element function spaces that ensure the (weak) continuity of flux at the interface between two coupled models (e.g., pore-scale and continuum-scale), and are the essential component of a highly flexible, efficient, and accurate non-overlapping domain decomposition method (Bernardi et al. 1994; Arbogast et al. 2000; Peszynska et al. 2002). Subsequently, Sun et al. (2012a) showed that mortars can be used as accurate upscaling tools for pore-scale models in obtaining macroscopic properties (e.g., permeability). They demonstrated that a large heterogeneous pore-scale domain can be decomposed along structural discontinuities and coupled via mortars to closely approximate the true permeability. Sun et al. (2012b) developed a single-phase reservoir simulator, in which Darcy grids in the near-well region were substituted with pore-scale models (Fig. 8ii). The study focused on upscaling strategies for the permeability field of the near-well region. In the foregoing studies, application of mortars was limited to linear, single-phase Newtonian flow without species transport and computational aspects were left unexamined. These issues were later addressed by Mehmani et al. (2012) and Mehmani and Balhoff (2014), who extended the use of mortars to non-linear (power-law) flow and (passive/reactive) solute transport and demonstrated computational efficiency and parallel scalability of their new algorithms (Fig. 8iii). The application of the above hybrid mortar methods is most appropriate when a tight coupling between spatial/temporal scales exists and a characterization of the pore-scale subdomain is available. This scenario is most likely to occur in a “skin-deep" region around wellbores (e.g., matrix acidization, and CO₂ leakage through wellbore cement), in which fluids are furthest away from equilibrium and direct access to pore-scale samples for geometric characterization is available.

Tartakovsky et al. (2008a) non-iteratively coupled the pore scale and the continuum in a diffusion-reaction system using an SPH formulation for both scales. They were able to show that their hybrid method produced reliable predictions, compared to modeling purely at the pore scale, of mixing-induced precipitation, which is typically not amenable to a purely continuum description as the reactions occur over the length of a few grain diameters (Fig. 8i). While the non-iterative nature of the method makes it very attractive from a computational standpoint, the method lacks flexibility in the sense that it requires both scales to be formulated using SPH. Advection was additionally ignored in that work. A novel overlapping method for coupling pore-scale inclusions to the surrounding continua was developed by Battiato et al. (2011) and was successfully verified for the case of Taylor dispersion in a reactive fracture. While the method is more general compared to the previous work (it included advection and was not limited to SPH), it appears to be limited to rather small pore-scale point inclusions that are dependent on the underlying macro-grid structure of the domain. The specific coupling algorithm used can also be rather computationally expensive, since the differential-algebraic system was solved by iterating between the differential and algebraic parts, rather than solving the system as a whole. Roubinet and Tartakovsky (2013) built upon the work of Battiato et al. (2011) and developed a non-overlapping hybrid method in 1D, in which the computational performance was enhanced by solving a global (instead of a sequential) differential-algebraic system. The method used finite volumes to discretize both the continuum- and the pore-scale subdomains, and introduced additional unknowns for each grid on either side of the interface. This dependence of interface unknowns on the number of interface grids can render the approach quite expensive for large problems. It is also noteworthy, that the approach seems to be a special case of the Global Jacobian methods of Ganis et al. (2012) and Mehmani and Balhoff (2014).

Chu et al. (2012) proposed an approach based on the heterogeneous multiscale method (HMM), in which macroscopic conservation equations are written assuming unavailability of constitutive flow equations (e.g., Darcy’s law) at the macro scale. The missing data were instead supplied from locally sampled pore-network simulations across the domain. Specifically, pore-scale models were used as providers of accurate in/out-fluxes at the boundaries of macroscopic control volumes. Chu et al. (2013) extended this work to two-phase flow problems, where pore-network models were used to track macroscopic two-phase drainage fronts (Fig. 8iv). While the method is a superior alternative to pure continuum-scale models in cases where local effects are dominant, its extension to problems with strong fluid–mineral reactions (where a macroscopic source/sink term is unknown and highly coupled to transport itself) is not obvious. Sheng and Thompson (2013) developed a dynamic coupling algorithm for two-phase flow, in which network models were embedded at the center of 1D macroscopic control volumes. The network models provided relative permeability values to the continuum model, while the continuum model provided fractional flow rates as boundary conditions to the network model. The network model was specially designed to handle these non-trivial boundary conditions. The large time-scale discrepancy between the pore scale and the continuum scale was addressed by performing successive steady-state simulations at the pore scale. The work showed that macroscopic and pore-scale saturations do not necessarily agree, and highlighted the difficulty in enforcing such an agreement. Tomin and Lunati (2013) developed a hybrid algorithm based on the multiscale finite volume (MsFV) method. In this approach the macroscopic domain is discretized into primal control volumes, upon which a complementary dual grid is superposed. Local pore-scale problems are then solved on the dual grid, in order to construct a set of bases as well as correction functions; these can be thought of as pore-scale closures for the continuum model. Global balance equations are then solved by forming a linear combination of these bases/correction functions and computing their corresponding multipliers. The continuum solution is then used to construct a conservative velocity field at the pore scale for advancing tracer concentrations and/or phase saturations. The hybrid method was successfully verified against fully pore-scale simulations for tracer transport and stable/unstable drainage problems.

The essence of all hybrid methods discussed is a “two-way communication” between the pore scale and the continuum. Hybrid methods are a relatively recent development compared to single-scale modeling strategies such as molecular dynamics (MD), pore-scale modeling, and reservoir simulation. Naturally, their development has been based upon the vast diversity of single-scale methods and various combinations thereof. Consequently, they themselves are very diverse in the ways they approach the problem, which might raise the obvious question as to which hybrid method is most appropriate for a given problem. Scheibe et al. (2015) have attempted to provide a general road map for choosing an appropriate hybrid strategy, referred to as the Multiscale Analysis Platform (MAP), depending on the degree of complexity of the hydrological problem at hand (i.e., degree of separability between temporal/spatial scales). MAP classifies various hybrid methods into separate “motifs”, and leads its user (with a specific application in mind) towards a suitable motif by asking a series of questions regarding the nature of the problem to be solved (Fig. 9). Such a classification additionally provides a useful context, within which the ever-increasing multiscale methods developed in the literature can be categorized.

Recent advances in rock imaging technology (e.g., XMT) have made it possible to resolve the complex pore-space geometry in sufficient detail, to the point that direct numerical simulations of various pore-scale processes are becoming commonplace and an area of active research. In addition, novel micromodel fabrication and measurement techniques (Werth et al. 2010; Karadimitriou and Hassanizadeh 2012) are allowing for an unprecedented possibility of quantitative comparison between modeling and experiments. However, our recent capability to resolve much of the pore-scale details with high accuracy, both experimentally and numerically, does not detract from the importance of mesoscale (or pore-network) models and the role they play in bridging the pore scale to the continuum. On the contrary, we can now use observations from such detailed experimental and direct numerical results to increase the accuracy of traditional network models and remove some of their ambiguity. The significance for doing so is threefold: a) it allows us to translate our detailed observations to an intermediate scale that is simple enough to understand, b) forces us to determine what features of the problem are the most essential so we can discard unnecessary details, and c) regardless of advances in computational performance, mesoscale models will always be one scale ahead of direct pore-scale methods. While improving network models based on detailed numerical/experimental observations is often successful, there are however limitations to this process. In the first part of the chapter, we tried to emphasize this point by providing two examples in the context of single-phase solute transport.

In the first example, we demonstrate that modifications to a traditional Eulerian network model allows it to capture partial mixing within pores with higher accuracy, while preserving computational efficiency and relative simplicity. In the same example, we also show that the need for accurately capturing partial mixing within pores depends on the porous medium. If the porous medium is ordered, partial mixing is important. But in disordered media (e.g., sandstones), partial mixing does not seem to be as significant, which means that we do not need to add this additional complexity to our mesoscale model. In the second example, we similarly show that modifications of a traditional Eulerian network allows it to capture shear dispersion within throats, and improves predictions of macroscopic dispersion in disordered (granular) media. However, when the porous medium is ordered, no amount of modification appears to prevent it from producing the wrong D_L vs. Pe_d scaling due to the Eulerian nature of the network model. Thus, a Lagrangian network model seems to be the next simplest mesoscale description. Similar conclusions may or may not hold for problems involving multiphase flow, geochemical reactions, adsorption, etc. However, the current state of computational power and experimental advances are allowing for the possibility of investigation into these questions more than ever.

On the other hand, the practical need to perform large-scale predictions requires pore-scale information to be translated to the continuum scale. While traditional upscaling approaches are successful in various scenarios, their applicability is rather limited when insufficient separation between scales exists. Overwhelming evidence in the literature hints towards this being quite common in the presence of strong fluid-mineral reactions; rendering averaged reaction rates transport-limited at the pore scale. Another example involves viscous and density-driven instabilities in multiphase and miscible displacement scenarios (Tartakovsky et al. 2008c). In recent years, various “hybrid” models have been developed as a means of simulating such scenarios, with the typical premise that the break-down of macroscopic continuum equations occurs locally. The essence of all such methods is a “two-way communication” between the pore and the macro scales, accomplished by incorporating models of both scales into the same computational domain. A useful classification given by Scheibe et al. (2015) facilitates one in choosing a suitable hybrid method, from a large selection developed in the literature, depending on the particular application at hand. As developments in both pore-scale and hybrid models continue into the future, the hope is for all these efforts to converge towards a deeper understanding of the scale transition problem and arrive at a predictive description (modeling and/or theory) of the relevant processes. Such an understanding would undoubtedly be indispensable in tackling currently pressing and highly challenging problems such as the safe sequestration and storage of anthropogenic CO₂, increasing hydrocarbon recovery from mature formations, and clean-up and remediation of nuclear waste repositories.

This material is based upon work supported as part of the Center for Frontiers of Subsurface Energy Security, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, and Office of Basic Energy Sciences under Award Number DE-SC0001114.