Computational Microfluidics for Geosciences

Computational microfluidics for geosciences is the third leg of the scientific strategy that includes microfluidic experiments and high-resolution imaging for deciphering coupled processes in geological porous media. This modeling approach solves the fundamental equations of continuum mechanics in the exact geometry of porous materials. Computational microfluidics intends to complement and augment laboratory experiments. Although the field is still in its infancy, the recent progress in modeling multiphase flow and reactive transport at the pore-scale has shed new light on the coupled mechanisms occurring in geological porous media already. In this paper, we review the state-of-the-art computational microfluidics for geosciences, the open challenges, and the future trends.


INTRODUCTION
Since their appearance in the early eighties, microfluidic experiments have revolutionized the knowledge in porous media research. They brought new insights into the debates that drive the community for decades including the controversial use of Darcy's law for modeling multiphase flow and the understanding of the complex feedback between hydro-bio-chemical systems. The most famous example is the work of Lenormand et al. (1988) who identified different flow instabilities in unsaturated porous systems according to the fluid pair and the flow rates using micromodels. These microfluidic devices enable direct visualization of the flow, transport, and chemical processes in rock microstructure replica sandwiched between two parallel plates separated by a small distance (Jahanbakhsh et al., 2020;Morais et al., 2020). Classically, a representation of the pore space is etched on a substrate made of silicon, glass or polymeric materials and covered by a transparent material. Recently, 3D printing has been used for fast and robust micromodel fabrication (Watson et al., 2019), allowing to conduct many flow experiments at lower cost for 2D to complex 3D geometry, although for now resolution is not as accurate as etched 2D micromodels. Coupled with high-resolution imaging both in space and time, microfluidic experiments have shed new light on the key mechanisms controlling subsurface flow. For example, particle image velocimetry has highlighted internal fluid eddies within droplets trapped by capillary forces whose large scale impact is not included in multiphase Darcy's law (Kazemifar et al., 2016;Roman et al., 2016Roman et al., , 2020Zarikos et al., 2018). Microfluidics also led to substantial advancements in the understanding of reactive systems as the most recent microchips involve reactive materials to investigate dissolution and precipitation of solid minerals in fractured porous media (Zhang et al., 2013;Song et al., 2014;Porter et al., 2015;Osselin et al., 2016;Soulaine et al., 2017Soulaine et al., , 2018Yoon et al., 2019;Agrawal et al., 2020;Poonoosamy et al., 2020;Yun et al., 2020). In just a few decades, microfluidics has become an indispensable tool in geosciences to decipher the complex mechanisms that occur in porous systems and the discipline keeps improving.
Computational microfluidics is the digital counterpart of microfluidic experiments. Like any modeling approach, it is a natural companion of laboratory experiments as it can be used to perform sensitivity analysis for identifying the key parameters of the underlying processes or to explore ranges of conditions including pressure and temperature that are difficult to reach in the laboratory without dedicated equipment. Moreover, computational microfluidics can also assist in the design of micromodels. Unlike pore-network modeling-a widely used approach for describing pore-scale physics (Fatt, 1956)-that relies on an approximation of the pore geometries and on pore-averaged physical laws, computational microfluidics solves the fundamental equations of continuum mechanics in the exact microstructure geometry. Historically, this kind of simulations was considered computationally expensive and almost impossible to run on complex and large domains. The enormous progress in high-resolution numerical simulation of complex flow and high-performance computing techniques, however, makes it possible to perform routinely computational microfluidics for single-phase flows on domains that reach the size of a Representative Elementary Volume . Although it heavily relies on advanced and modern approaches for simulating fluid dynamics, computational microfluidics for geosciences has its own challenges related to the complexity of natural systems including the heterogeneity of geological porous media, the evolving porous microstructure along with geochemical processes, and the complex description of interfacial phenomena at the mineral surfaces (Meakin and Tartakovsky, 2009).
The recent improvements in computational microfluidics for geosciences enable the modeling of multiphase flow and reactive transport. Current capabilities capture many of the mechanisms observed experimentally. For example, two-phase flow simulations in micromodels describe the complex interplay between viscous, capillary, and gravity forces leading to viscous fingerings if a fluid drains another fluid (Ferrari and Lunati, 2013;Alpak et al., 2016;Chen et al., 2019;Zhao et al., 2019), or to the trapping of carbon dioxide droplets in the pore space (Chen et al., 2015;Hu et al., 2017). Computational microfluidics allows the understanding of subsurface processes by identifying key mechanisms. Hence, the integrated usage of computational and experimental microfluidics led  to highlight that snap-off events, i.e., the break-up of the invading fluid at the entrance of a pore-throat, can be followed by a reconnection to the flowing phase. They observed and characterized cycles of snap-off and reconnection events that contradict the criteria that are generally used to assess the storage capacity by capillary trapping. Another example is the study of Ferrari and Lunati (2014) which establishes that inertial forces arise during pore invasion-a statement also supported by glass bead experiments (Moebius and Or, 2012) and interface tracking in micromodels )suggesting that creeping flow assumption for the derivation of multiphase Darcy's law is erroneous. Moreover, as computational microfluidics provides a mapping of pressure, velocity, species concentration, and mineral distribution, the characterization of constitutive laws using volume averaging is straightforward. For instance, Soulaine et al. (2017) used computational microfluidics to compute the accessible reactive surface area of a rock sample. A more comprehensive review of the state-of-the-art including the authors' contributions, as well as the remaining challenges that need to be solved for being truly predictive, are presented in section 3.
Our definition of computational microfluidics goes beyond the modeling of flow and geochemical processes in micromodels as it can also be applied to complex 3D structures. In this sense, computational microfluidics for geosciences is intimately linked to Digital Rock Physics, an emerging technology concerned with the estimation of effective upscaled parameters (e.g., permeability, dispersion tensor) in three-dimensional images of rock samples . In Digital Rock Physics, these properties can be calculated using different techniques, such as pore network analysis (Andre et al., 2013a,b) or closure problem approaches (Whitaker, 1999;Soulaine et al., 2013). Alternatively, computational microfluidics can be used to directly solve the flow and transport equations within the pore-space and compute the desired upscaled properties. Although it is often restricted to small micro-CT images (<300 3 voxels) or to single-phase flow due to high computing cost, computational microfluidics is also well-suited to investigate pore-to-pore physics (Ferrari and Lunati, 2014;Pavuluri et al., 2020) with the objective of improving upscaling techniques (e.g., pore entry pressure for pore network modeling). A major asset of microfluidics is the ability to design well-controlled experiments and obtain high-resolution measurements-both in space and time-of the velocity fields , phase distribution, film thickness , aqueous species concentrations (Chang et al., 2017), and mapping in mineral changes (Poonoosamy et al., 2020). Microfluidics is inestimable to verify the numerical models using direct comparison with microfluidic experiments (Willingham et al., 2008;Yoon et al., 2012;Chapman et al., 2013;Oostrom et al., 2016;Roman et al., 2017;Soulaine et al., 2018). Therefore, the combination of computational and experimental microfluidics promotes the development of robust and validated numerical models for porescale modeling.
Exposing simulation results to experimental data, however, is not always an easy task as it raises questions regarding the measurement resolution and also questions on the modeling spatial description. Micromodels are often referred to as "twodimensional" objects by opposition to real geological objects or their reproduction using 3D printing. This qualification is somewhat imprecise as the micromodels are made of the assembly of two parallel plates and that important physical phenomena occur in the confined space between the plates. This confusion leads many scientists to use two-dimensional models for representing their microfluidic experiments numerically. In most cases, however, 2D simulations are not able to reproduce the experiments. Even in cases for which the micromodel depth is much larger than the smallest pore-throat diameter the twodimensional hypothesis breaks down because physical events are not uniform in the thickness. For example, reactive particles can nucleate in the mid-plane or at the solid walls (Beuvier et al., 2015) and instabilities in the liquid films that wet the top and bottom plates can lead to the formation of bubbles in the bulk (Hansen and Toong, 1971;Zhao et al., 2018). The systematic usage of 3D models, however, is computationally intensive. To reduce the computational efforts, 2D depth-integrated models are an interesting alternative for solving flow in microfluidic chips. These approaches, however, can have strong limitations. In section 2, we present depth-integrated models for single-phase, two-phase, and species transport, and we discuss their range of validity.
This paper aims at reviewing the modeling approaches in computational microfluidics for geosciences and their open challenges, as well as the authors' contributions to this emerging scientific discipline. It is organized as follows. First, we introduce the fundamental equations used in computational microfluidics (section 2.1) and the common numerical approaches to solve them (section 2.2). We also present models that solve directly for the depth-averaged variables at a reduced simulation cost (section 2.3). Then, we summarize the main advances in modeling multiphase flow (section 3.1) and reactive transport (section 3.2), and we discuss future lines of research.

MATHEMATICAL MODELS AND NUMERICAL SOLUTIONS
In this section, we introduce the fundamental equations used in computational microfluidics for solving flow and transport. We also discuss the standard numerical approaches for solving these equations.

Fundamental Governing Equations
Computational microfluidics relies on the fundamental equations of continuum mechanics for fluid flow, namely the Navier-Stokes equations, solved in the exact geometry of the pore space. They are composed of a set of mass and momentum balance equations combined with boundary conditions at the solid surface.

Flow Equations
For single-phase flow, the conservation of mass reads, where ρ and v are the fluid density and velocity, respectively. For liquids, it is common to assume that the fluid is incompressible, and the mass balance equation (Equation 1) reduces to the divergence-free velocity, In an Eulerian frame, the momentum balance writes, where p is the pressure field, g is the gravity, and µ is the fluid viscosity. The left-hand side corresponds to the inertial effects. As mass flow rates in the soils and the subsurface are usually low, it is common to neglect this term and to use Stokes instead of the Navier-Stokes momentum equation. On the right-hand side, the first term is the pressure gradient, the second is the gravity acceleration and the third term is the viscous dissipation. The latter characterizes the internal friction of the fluid particles with each other. The higher the viscosity, the larger the pressure gradient has to be to initiate the fluid movement. At the porescale, it is also common to neglect the gravity term. For multiphase flow, the mass and momentum equations, Equations (1)-(3) apply to each fluid phase. Alternative forms of the Navier-Stokes equations including non-Newtonian and compressible fluids as well as descriptions of the equations in a Lagrangian frame are found in any fluid mechanics textbook (e.g., Landau and Lifshitz, 1987;Caltagirone, 2013).

Boundary Conditions at the Solid Walls
The flow equations are supplemented with boundary conditions at the solid surface that arise from particle/solid molecular interactions. Because the fluid molecules cannot penetrate into the impermeable boundary, the mass flux at the solid surface is null and the normal component of the velocity is zero, i.e., n.v = 0. For viscous flow, it is common to consider that the fluid molecules adhere to the surface due to molecular interactions including van der Waals or Coulombic forces. Because of this adherence, the fluid particles at the solid wall cannot slide on the surface and the tangential component of the velocity is null, leading to the no-slip condition, v = 0. (4) The no-slip boundary condition is valid for most of the Newtonian liquids flowing through porous media with porethroats larger than 100 nm. Although widely used in fluid dynamics, the no-slip condition is not always valid in geological systems. In particular, in complex fluids involving colloids or polymers, the fluid elements adjacent to the surface can travel along the surface by rolling or sliding. In such cases, the fluid in contact with the solid surface slide over it with a finite value. Moreover, the no-slip condition poses a problem in viscous flow theory at the place where an interface between two fluids meets a solid boundary. Instead, the so-called partial slip condition, proposed by Navier (1823) and later revisited by Maxwell (1879), relates the tangential velocity to the strain rate. In this relation, n is the normal vector at the solid surface and β is called the slip length and is typically of the order of ten nanometers. The description of boundary conditions at the solid walls for multiphase flow is discussed in section 3.1.

Interfacial Conditions Between Two Fluids
The interface between two fluid phases (labeled fluid 1 and fluid 2, respectively) is a discontinuity on which boundary conditions are applied. The interface deforms with the fluid motion and moves at the interface velocity, w. The latter is obtained using the conservation of mass at the vicinity of the interface, where n 12 is the normal at the interface separating the two fluids, and ρ i and v i are the density and velocity of phase i, respectively. If there is no phase change, the interface velocity and the fluid velocity balance each other at the interface and n 12 .v 1 = n 12 .v 2 . The normal stress balance at the interface gives (Landau and Lifshitz, 1987), where p i is the fluid pressure in fluid i, σ is the surface tension, and κ is the interface curvature. At equilibrium, this condition becomes Laplace law. An overview of the simulation techniques used in computational microfluidics to solve Navier-Stokes equations along with the interfacial conditions is given in section 3.1.

Solute Transport
In its simplest form, the transport of a species i in a fluid phase is modeled by an advection-diffusion equation: where the first term is the accumulation term, the second term is the advection term and the right-hand side corresponds to diffusion effects where D i is the diffusion coefficient of species i into the fluid phase. More complex forms of the solute transport includes multi-component diffusion, thermal diffusion and electrostatic gradients (Bird et al., 2002). As for flow equations, the solute balance equation is supplemented by boundary conditions at the solid walls. For a passive tracer the no flux condition writes, If the species reacts with the solid surface the reactive flux condition is, where r m is a source or sink term that describes kinetic reactions or adsorption phenomena (Molins et al., 2014). It is also common to consider that the species satisfies thermodynamics equilibrium condition with the solid wall using, where C eq i is a fixed value imposed by thermodynamics equilibrium.

Numerical Engines for Solving Navier-Stokes Equations
Solving Navier-Stokes equations in complex geometries is not trivial, and except in very few situations, there are no analytical solutions. Computational microfluidics uses numerical approximations to solve the flow within the pore space. Complexities include the pressure-velocity coupling and the nonlinearity due to the inertial term in the momentum equation. In this section, we introduce some of the standard approaches for solving Navier-Stokes equations. Specific techniques for multiphase flow and reactive transport are introduced and discussed in section 3.

Computational Fluid Dynamics
In Computational Fluid Dynamics (CFD), the Navier-Stokes equations are solved by discretizing the spatial differential operators on an Eulerian grid using techniques such as the finite difference method (FDM), the finite volume method (FVM), or the finite element method (FEM) (Ferziger and Peric, 2002). In principle, the three families of techniques yield the same solution if the computational grid is fine enough. For fluid flow, however, FVM is usually preferred as it is locally conservative and deals with unstructured grids by construction. In FVM, the domain is subdivided into a finite number of contiguous control volumes whereby the integral form of the conservation laws is applied. Modern gridders generate meshes that follow the complex boundaries of the microstructure (see Figure 1). Nowadays, all the commercial or open-source CFD simulators embed very efficient Navier-Stokes solvers that can deal with very large grids and complex geometries using multigrid strategies and High-Performance Computing. Solution algorithms usually rely on predictor-corrector methods to solve for the pressure-velocity coupling. The Pressure-Implicit with Splitting of Operators (PISO) algorithm is one of the most popular for solving transient flow (Issa, 1985). Algorithms such as Semi-Implicit Method for Pressure Linked Equations (SIMPLE) solve directly for steadystate flow (Patankar, 1980).

Lattice Boltzmann Method
The lattice Boltzmann method (LBM) is an alternative CFD technique to solve fluid flow at the pore-scale. The popularity of LBM comes from the ease to program massively parallel codes. Unlike traditional CFD approaches that solve for Navier-Stokes equations, LBM is based on the Boltzmann equation that describes the probability of finding a given particle at a given position with a given velocity. This is a statistical approach and the averaged behavior approximates the continuum mechanics described by the Navier-Stokes equations. In LBM, a discrete Boltzmann equation is solved on a lattice where the particles can only move along a finite number of directions. The degree of freedom is determined by the nature of the lattice classified with the nomenclature DnQm where n corresponds to the dimension of the lattice (n = 2 for two-dimensional simulations, n = 3 for three-dimensional simulations) and m stands for the number of directions. The Boltzmann equation involves a collision operator that describes the change in the probability distribution function induced by the collisions between particles. In the Lattice Boltzmann method, the boundary conditions are described in terms of probability distribution functions. The description of the wall conditions corresponds actually to the physical model of boundary events proposed by Maxwell (1879). For a no-slip, the bounce-back condition is used where the particles are reflected at the wall. Hence, the net average flow rate at the solid surface is zero. Slip conditions can be obtained using specular reflections (Succi, 2002). Partial slip is modeled with a mix of bounceback and specular reflections. The implementation of boundary conditions in LBM can become very complex when the solid surface is not aligned with the lattice.

Smoothed-Particle Hydrodynamics
Smoothed-Particle Hydrodynamics (SPH) is a mesh-free modeling technique that describes fluids as a set of discrete moving particles in a Lagrangian framework. The method was initially developed to model the interaction of stars in astrophysics (Gingold and Monaghan, 1977;Lucy, 1977) and is now being used increasingly in real-time animation and video games. The particles have a spatial distance over which their properties are smoothed by a kernel function centered on the particles. This distance defines the sphere of influence for each particle and is used to reconstruct continuous variables such as fluid density. The SPH particles are then used as interpolation points to discretize and solve the governing equations. SPH has several benefits over traditional grid-based methods. First, it is meshless and there is no need to go through complex gridding processes. Second, it explicitly conserves mass and momentum. Because the method is Lagrangian, there is no non-linear term in the momentum equation and that allows an accurate description of the advection terms when dealing with high flow rates. Third, the SPH Lagrangian framework allows for the modeling of moving and deformable boundaries without using complex tracking algorithms. SPH has also important drawbacks. Although the SPH algorithms due to their simplicity and structure are easy to program, the implementation of boundary conditions can be tricky. Moreover, for the same accuracy, SPH is more computationally intensive than FVM or FEM because it requires much more neighboring particles than the tensile of grid-based methods. Nevertheless, because SPH can be easily paralleled, the rise of parallel computing including Graphical Processing Units makes SPH a promising technique (Vacondio et al., 2020).

Depth-Averaged Models: 2D vs. 3D
Most of the time, measurements in microfluidic experiments relate to depth-averaged values. For example, even though micro-Particle Image Velocimetry that uses confocal microscopy or laser pretends to measure the velocity profiles in a cross-section, these techniques do measure values averaged in a thickness of the order of few micrometers that corresponds to the laser thickness or the close surrounding of the focal plan . Likewise, microfluidic experiments using optical luminescence for measuring a concentration profile (Watson et al., 2019) or mapping the water distribution (Zhao et al., 2016;Roman et al., 2017) provides values that are averaged along the micromodel depth, h, without information of the depth profile of the species concentration or water saturation. The depth-averaged variables are defined as the vertical integration of the 3D fields, In this section, we present depth-averaged models that directly provide the solution for the depth-averaged velocity, v 2D , pressure, p 2D , and concentration, C i 2D , fields without solving 3D problems. These models result from the integration of the fundamental equations described in section 2.1 along the depth axis and can significantly reduce the simulation costs. They have, however, ranges of validity that are discussed below.

Single-Phase Flow
The integration of Navier-Stokes momentum equation over the micromodel depth along with the assumption that a parabolic profile approximates the flow in the thickness gives, where the last term of the right-hand side, µ 12 h 2 v 2D , is the Hele-Shaw correction (Lamb, 1906), i.e., a drag force related to the friction of the fluid along the micromodel top and bottom plates. Numerical implementation of Equation (13) has been reported using CFD  and LBM (Venturoli and Boek, 2006;Laleian et al., 2015). The Hele-Shaw correction is similar to a Darcy law for which the permeability is k = h 2 12 and the porosity is φ = 1. This simple depth-integration was historically the first step toward the derivation of Darcy's law from Stokes equation (Vasil'ev, 2009). Actually, for micromodels with a low aspect ratio between the cell thickness and the pore-throat size, h/d ≪ 1, the drag force is the dominant force and Equation 13 can be replaced by Darcy's law for Hele-Shaw cells, Theoretically, the hypothesis that flow sandwiched in between two parallel plates obey a parabolic profile is only valid for infinite plates-their dimensions are much larger than the space FIGURE 2 | Single-phase flow simulations using 3D, 2D, and 2D-depth-averaged models for different aspect ratio between the micromodel thickness, h and the smallest pore-throat, d. The computational domain is gridded with 1.4 × 10 5 cells for the 2D cases and 11.4 × 10 5 × 10 for the 3D cases. The typical pore-throat size is d typ = 4 × 10 −3 m, the largest is d large = 12 × 10 −3 m and the smallest is d = 0.5 × 10 −3 m. We used two different depth thickness, namely h = 10 −2 m and h = 10 −4 m. Simulations are performed with OpenFOAM. For each simulation, the velocity profile is normalized by the maximum value, and we indicate the maximum value normalized by the inlet velocity. If h d > 1, 2D Navier-Stokes simulations tend toward the mid-plane 3D Navier-Stokes solution and 2D Darcy is far from the reference solution. If h d < 1, 2D Darcy captures the 3D reference solution whereas 2D Navier-Stokes breaks down. The 2D depth-integrated Navier-Stokes model captures the 3D flow behavior regardless the micromodel thickness.
between them. Micromodels mimicking rock porous structures, however, contain pillars that violate this assumption with threedimensional flowlines at the vicinity of the corner formed by the pillars and the top and bottom plates. Nonetheless, in practice, Equation (13) is a very good approximation of the flow behavior in Hele-Shaw cells. In Roman et al. (2016), the very good agreement in the comparison of the velocity profiles obtained experimentally with micro-PIV and numerically solving Equation (13) illustrates that the depth-averaged flow model is relevant for micromodels. This statement is verified numerically in the simulation results presented in Figure 2 for which the velocity profile in a micromodel was computed based on 2D, 2D depth-integrated, or 3D simulations for different micromodel thickness. For 3D simulations, the 2/3 factor in the velocity magnitude of the 3D depth-integrated profile with respect to the mid-plane value is representative of a parabolic profile in the third dimension. The 2D Navier-Stokes model only tends to the 3D solution if h/d > 1 and corresponds to the 3D mid-plane profile. In this case, the scales separation hypothesis leading to Darcy's law is not satisfied and the 2D Darcy solution cannot describe the flow through the micromodel. If h/d < 1, the friction of the fluid along the top and bottom plates are the dominant resistance to the flow. In such a case, the 3D flow tends to the 2D Darcy solution whereas the 2D Navier-Stokes solution-that corresponds to an infinite thicknessis completely off the reference data. The 2D depth-averaged simulations remain consistent with the 3D solution regardless of the aspect ratio.

Single-Phase Scalar Transport
The integration of an advection-diffusion equation for the transport of a scalar-e.g., species concentration or temperatureacross the micromodel thickness gives, where v 2D is obtained using Equation (13), and D(t) is a time-dependent dispersion tensor. Because of the parabolic velocity profile in the micromodel depth, there are hydrodynamic dispersion effects that must be considered even if the aspect ratio between the pore-throat size and micromodel thickness is lower than one, h/d > 1. For transport in between two parallel plates, the dispersion tensor exactly writes (Dentz and Carrera, 2007), where the first term is the molecular diffusion and the second term corresponds to Taylor-Aris dispersion in the direction of the flow (Taylor, 1953;Aris, 1956)-notes that v 2D v 2D is a tensor. The FIGURE 3 | Transport of a scalar in a micromodel with aspect ratio h/d > 1. We use 3D, 2D, and 2D-depth-averaged models for different Péclet number by varying the diffusion. The results correspond to isocontour of the concentration for C = 0.2 (blue), C = 0.5 (green), C = 0.8 (orange) at t = 42 s. 2D models without dispersion break down for Pe h > 10.
time-dependent infinite series accounts for the pre-asymptotic regime at early time and tends toward zero at long time. In practice, only the two first terms of the series are necessary. Hydrodynamic dispersion effects are characterized by the Peclet number, Pe h = hU D , that compares the transport by advection and the transport by diffusion. They are illustrated in the simulation results of Figure 3 that solve the transport of a scalar using 3D, 2D, and 2D depth-averaged models for Peclet numbers ranging from 0.1 to 10 4 . The deviations between midplane and depth-averaged values of the 3D simulations illustrate the hydrodynamic dispersion when advection is the dominant transport mechanism. Experimentally, the measurement of species concentration using optical index corresponds to depthaveraged values (Watson et al., 2019). Confusion regarding the nature of measured value may lead to misinterpretations of the comparison between experiments and 3D simulations. For transport dominated by diffusion, Pe h < 1, the concentration profile is flat over the micromodel depth. Therefore, there are no dispersion effects and the 2D models are in good agreement with 3D simulations. If the advection takes over, Pe h > 1, hydrodynamic dispersion becomes important and 2D model breaks down. The two-dimensional depth-averaged simulations using the time-dependant dispersion tensor approximate well the depth-integrated solution of the 3D simulations, regardless the Peclet number. Figure 4 highlights the importance of considering the pre-asymptotic regime. Neglecting the dispersion (i.e., D(t) = D i I) or considering only the asymptotic Taylor-Aris dispersion (i.e., D(t) = D i I + h 2 v 2D v 2D 210D i ) tends to overdisperse the solute. Likewise, it is crucial to consider the product v 2D v 2D as a 210Di ] tends to overdisperse the solute. Isotropic dispersion cannot capture the transport of scalars in micromodels.
tensor and not as a scalar v 2D 2 because Taylor-Aris dispersion in microfluidic devices is not isotropic and occurs in the direction of the flow.

Two-Phase Flow
Simulation of two-phase flow processes in microfluidic cells using 2D depth-averaged models is challenging as the fluidfluid interface has two curvature radii: one in the micromodel plane, another in the micromodel thickness. Horgue et al. (2013) proposes 2D depth-averaged equations in which they integrated the interface curvature along the micromodel thickness in addition to the Hele-Shaw correction. Their model, however, is limited to cases for which no liquid films are flowing on the top and bottom plates, and cannot capture phenomena such as snapoff. Alternatively, Cueto-Felgueroso and Juanes (2014) proposes a Darcy-like model to simulate two-phase flow in Hele-Shaw cells that accounts for water films flowing on the top and bottom plates using relative permeabilities. The applicability of their model, however, is restricted to aspect ratios between the micromodel thickness and the pore-throat size smaller than one, h/d ≪ 1.
In absence of reliable 2D depth-integrated models for two-phase flow, it is recommended to use 3D rather than 2D simulations since the former can capture the impact of wetting films along the top and bottom plates as well as corner flows (Zhao et al., 2019). Moreover, phenomena such as snap-off can only be captured if the two curvature radii are considered .

Multi-Scale Approach
Multi-modal pore-size distributions are quite common in geological porous media. This special feature appears, for example, in porous fractured media for which the width of a fracture is orders of magnitude larger than the typical pore throat diameter in the porous matrix. In such a case, a full Navier-Stokes modeling approach requires to resolve all the porosity explicitly and leads to tremendous computational grid size. Alternative approaches including micro-continuum models propose to model only the flow in the fracture using Navier-Stokes while the flow in the matrix is described by Darcy's law . Micro-continuum approach relies on a hybrid-scale flow model based on the Darcy-Brinkman-Stokes equation (Brinkman, 1947;Neale and Nader, 1974) discretized with the Finite-Volume Method. All the geometric structures below a threshold that corresponds to the spatial resolution of the simulation is filtered and modeled as a porous medium. The momentum equation reads where the overline notation denotes filtered variables, φ is porosity field, and k is the local permeability of the porous regions. Recent studies show that the micro-continuum framework is well-suited to solve problems involving two characteristic length-scales (e.g., fractured media, micro-porous rocks) (Arns et al., 2005;Scheibe et al., 2015;Soulaine et al., , 2021Carrillo et al., 2020;Menke et al., 2020;Poonoosamy et al., 2020) and also to describe the movement of fluid/solid interfaces caused by dissolution/erosion and precipitation/deposition Soulaine et al., 2017;Molins et al., 2020) or by solid deformation due to swelling or fracturing (Carrillo and Bourg, 2019). Microcontinuum models are also able to simulate two-phase flow processes (Soulaine et al., 2018(Soulaine et al., , 2019Carrillo et al., 2020).

APPLICATION TO GEOSCIENCES AND OPEN CHALLENGES
In this section, we review the main advances in computational microfluidics involving multiphase flow and reactive transport including the authors' own contributions. We also discuss the open-challenges and the current and future lines of research in this area.

Modeling of Immiscible Two-Phase Flow
Modeling immiscible multi-phase flow is one the biggest challenge in the porous media community. Despite a very large usage in subsurface engineering, many studies have emphasized the limits of Darcy's law for two or more fluids (Danis and Quintard, 1984;Rose, 2000;Cinar et al., 2009). An accurate description of the flow mechanisms in porous media is challenged by the complex interplay between capillary, viscous, and gravitational forces that may lead to highly unstable flow and trapping mechanisms. The heterogeneous nature of geological formations further complicates the development of predictive models.
Computational microfluidics for multiphase flow offers an appealing framework to probe the underlying mechanisms in unsaturated porous media. Despite enormous improvements, however, the numerical modeling of multi-phase flow remains challenging. The difficulties are twofold: on the one hand, numeric must be improved (for calculating the interface curvature dynamically, for tracking the composition of the fluids, for speeding-up the simulation times); on the other hand, the physical description of interfacial properties is still not well-established (e.g., wettability conditions, thin films). Current research intends to address these open-challenges by developing reliable and efficient multi-phase flow simulators at the pore-scale.

Efficient Numerical Solvers for Capillary-Dominated Flow
The standard approaches for tracking the interface movement of two immiscible fluids rely on solving single-field equations for the single-field variables, i.e., variables that includes the contribution of each phase. Whether the modeling is gridbased or particle-based, a phase indicator function maps the distribution of the fluids within the pore-space and evolves as the fluid-fluid interface propagates along with hydrodynamic forces. This function denotes either the fluid volume fraction in every grid cells for the Volume-of-Fluid (VOF, Hirt and Nichols, 1981) and phase-field (PF, Jacqmin, 1999) approaches, a distance to the fluid interface for Level-Set (LS, Sussman et al., 1994) approaches, or a color function for Smooth Particles Hydrodynamics (SPH, Tartakovsky and Meakin, 2005) and Lattice Boltzmann Methods (LBM, Gunstensen et al., 1991). The literature also reports on modeling based on the Density Functional Theory for hydrodynamics (Demianov et al., 2014). A review of approaches is found in Wörner (2012). The phase indicator function is used to define single-field variables by weighting the fluid properties (e.g., viscosity, density) according to the distribution of fluid phases. Its gradient is also used to calculate the interface curvature. Today, none of the existing approaches can solve two-phase flow problems accurately and efficiently for a wide range of flow regimes (Abadie et al., 2015). One of the bottlenecks is the existence of spurious velocities near the interface attributed to inaccuracies in the calculation of the interface curvature, regardless of the computational approach. This problem acknowledged for more than 30 years (Brackbill et al., 1992), is particularly limiting for capillary-dominated flow (Scardovelli and Zaleski, 1999) that are the dominant regimes in many subsurface systems (Cinar and Riaz, 2014). An important part of the current research in computational microfluidics for two-phase flow consists in developing robust solvers that minimize these numerical issues (Pavuluri et al., 2018). One of the popular approach in computational microfluidics is the algebraic VOF approach. This is the technique used in the simulations that illustrate the application section. In this technique, the phase indicator functions, α l and α g , described the volume fraction of each fluid within each cell of the computational domain and the governing equations are obtained by volume averaging the Navier-Stokes equations over a control volume (Maes and Soulaine, 2020). For incompressible fluids, the continuity equation writes, the mass balance equation for the liquid phase is and the momentum equation writes where v r and F c are sub-grid models that describe the fluid-fluid interactions in cells that contain the interface. They represent the relative velocity and the surface tension force due to the curvature of the interface, respectively. Although other methods such as CFD with level-set (Sussman et al., 1994;Abu AlSaud et al., 2017 or combinations of LBM and level-set (McClure et al., 2016) can provide a more accurate description of the sharp interface, the VOF method is attractive due to its flexibility, robustness in terms of mass conservation, and adaptability to more complex physic (e.g., multi-component mass transfer). Because the derivation of the governing equations relies on volume averaging principles, algebraic VOF methods are the cornerstone of multi-scale solvers for two-phase flow using micro-continuum (Soulaine et al., 2018(Soulaine et al., , 2019Carrillo et al., 2020).

Modeling of Wettability and Interfacial Properties
Wettability plays a key role in multiphase flow in porous media (De Gennes, 1985;Singh et al., 2018). Accurate prediction of multiphase flow in the soils and the subsurface for geoenvironmental issues including Carbon Capture and geological Storage, water resources remediation, or Enhanced Oil Recovery requires an in-depth understanding at the pore-scale of the key mechanisms that influence the fluid displacements and dynamics including viscous and capillary forces but also the solid surface wettability. Recent studies have shown that the spatial distribution of wettability play a crucial role (AlRatrout et al., 2018;Yesufu-Rufai et al., 2020). The rock wettability and the resulting multiphase flow can be modified by a change in pore water composition (pH, salinity) and interfacial physicochemical properties. For example, an oil droplet trapped within a pore may be more easily remobilized using low salinity water (Jackson et al., 2016). Accurate modeling of wettability conditions at the mineral surface is therefore crucial to develop efficient and reliable computational microfluidics for two-phase flow. Traditionally, wettability is described considering that the fluid-fluid interface and the solid surface form a contact angle (see Figure 5). The value of the contact angle is provided from static measurements that described the contact angle at the nanoscale. In computational microfluidics, the cell size near the solid surface is usually orders of magnitude larger than nanometers. Therefore, upscaled laws including the Cox-Voinov model estimates an apparent dynamic contact angle that depends on the interface velocity near the solid surface and accounts for viscous bending (Voinov, 1976;Cox, 1986). The alteration of interfacial and solid surface properties due to a change of pH and salinity requires constitutive laws for modifying the contact angle accordingly (Maes and Geiger, 2018). Moreover, flow models based on contact angles break down in the presence of thin water films on the mineral surface (Starov et al., 2007).
Alternative flow models including lubrication theory offer an appealing framework to account for physicochemical mechanisms controlling the wettability of rocks according to pore water chemistry and surface roughness for instance (Abu AlSaud et al., 2017. Instead of using a contact angle, lubrication theory provides a partial differential equation governing the evolution of the film thickness, δ, on the solid surface including molecular forces . It reads, where µ w is the water viscosity, σ is the surface tension, κ is the interface curvature, (δ) is the disjoining pressure term, and τ is the shear-stress exerted on the film surface. The disjoiningor Derjaguin-pressure includes attractive (e.g., van der Waals) and repulsive (e.g., electrostatic) forces that change with the ionic strength, i.e., pH and salinity. The film equation is solved on the solid surface and is coupled with the two-phase Navier-Stokes solver in the volume of the pore. By getting rid of the concept of contact angle, it is believed that simulators using lubrication models are more realistic and mechanistic for the description of multiphase systems in geological formations for energy and environmental issues. Applications include the modeling of wettability alteration measured in supercritical CO 2 stored in deep saline aquifers (Wan et al., 2014) and the detachment of crude oil from mineral surface observed when changing the brine from high to low salinity (Berg et al., 2010). The state-of-the-art computational microfluidics using lubrication theory, however, is limited to flat solid boundaries (Abu AlSaud et al., 2017. Recent work considers curvilinear surfaces on cylindrical solids (Qin et al., 2020). However, pore-scale multiphase flow solvers with lubrication models on generic solid geometries are still under development.

Modeling of Species Transfer Across Fluid-Fluid Interfaces
In many subsurface processes including the injection and sequestration of supercritical carbon dioxide (sCO 2 ) into deep geological formations, there is a transfer of species across the fluid-fluid interface. Moreover, the ability to track the composition of fluids in two-phase flow is crucial for modeling efficiently the change of wettability conditions with pH evolution (Maes and Geiger, 2018). The Continuous Species Transfer (CST) approach allows tracking multicomponent multiphase flow at the pore-scale (Haroun et al., 2010;Marschall et al., 2012). In this approach, the thermodynamic equilibrium at the interface is described by a partitioning relationship such as Henry's law (Henry, 1803) and the continuity of mass fluxes is guaranteed while the interface evolves in the pore-space due to fluid displacements including drainage and imbibition (see Figure 6).
CST is a single-field approach along the same lines as the Volume-Of-Fluid technique used for simulating immiscible twophase flow, i.e., a unique partial differential equation solves for a single-field concentration regardless of the content of a grid block. It reads, where D * A is the weighted-average of the diffusion coefficients in both fluids and A is the CST flux enforcing the concentration discontinuity due to Henry's law. Graveleau et al. (2017) derived a boundary condition on a wetted wall for the CST and used the modeling approach to upscale the rate of mass transfer in unsaturated porous media as a function of the flow rates. Yang et al. (2017), however, pointed out that when convection dominates diffusion locally near the interface, the CST method generates large numerical errors. They show that in order to capture the species concentration discontinuity accurately, the mesh should be extremely refined at the vicinity of the interface. This workaround, however, increases dramatically the computational cost, especially for advection-dominated systems.
The CST derivation was carefully revisited using the method of volume averaging (Maes and Soulaine, 2018b) leading to a new approach without limitation on the transport regimes or the grid refinement level-referred to as Compressive-CST or simply C-CST-that is fully consistent with the phase advection scheme unlike the standard CST scheme that was causing unbalanced mass transfer. C-CST can simulate the mass transfer of a passive scalar in a dynamic system (see Figure 7), and therefore, enable the assessment of the mean rate of mass transfer in unsaturated porous media in terms of saturation, flow rates and interfacial area. C-CST can be used to evaluate the role played on the mass transfer by the internal re-circulation within the trapped droplets that was observed experimentally using microfludics Zarikos et al., 2018).
In more recent work, the C-CST has been extended for modeling cases for which the species transfer across the fluidfluid interface is no longer passive but leads to local volume changes (Maes and Soulaine, 2020). Such situations happen, for example, when a soda can is opened, the pressure in the headspace falls abruptly and Henry's law is no longer satisfied. To re-establish equilibrium, the carbon dioxide concentration in water decreases by forming gas bubbles. During their travel toward the surface, gas bubbles grow because of the diffusive flux from the carbon dioxide dissolved in the liquid (Power et al., 2009). The other way around occurs in carbon sequestration processes when sCO 2 is stored in deep saline aquifers by capillary effects and trapped sCO 2 droplets dissolve into the surrounding brine (Kim et al., 2012). As a consequence, the volume of droplets decreases, and the equilibrium of hydrodynamic forces that was maintaining the sCO 2 trapped in the pore-space is displaced, eventually leading to a mobilization of the droplets. Computational microfluidics using C-CST for phase change can estimate the volume of CO 2 trapped in the porous structure and assess the timescale of potential remobilization.

Multi-Scale Two-Phase Flow
The micro-continuum approach for two-phase flow enables the simultaneous modeling of multiphase flow at two different length scales: (i) a Darcy-scale where sub-voxel fluid-fluid and fluidsolid interactions within a porous medium are modeled through relative permeability and capillary constitutive models, (ii) a pore-scale (or Navier-Stokes scale) where the solid material is non-porous and fluid-fluid interactions are described through a continuum representation of the Young-Laplace equation. The model combines the Volume-of-Fluid approach for tracking the fluid-fluid interface in regions free of solid and the Darcy-Brinkman-Stokes equation for modeling flow in porous regions (Soulaine et al., 2019). It includes classic concepts such as saturation, relative permeability, and capillary pressure curves in the porous matrix (Carrillo et al., 2020). This multiscale approach offers new possibilities in computational microfluidics to investigate multiphase flow in fractured porous media including the fracture-matrix interactions and to study the role played by microporosity in processes such as imbibition and drainage.
For example, the two-phase micro-continuum model was used to simulate the expulsion of methane stored in microporous matrix into a soaked fracture (see Figure 8). The results highlight two completely different expulsion flow regimes whether the fracture boundaries are hydrophobic or hydrophilic. On the one hand, there is a continuous gas production, and on the other hand a bubbly regime because the surface tension effects at the wet fracture boundary generate an energetic barrier that tends to form gas bubbles as far as the pressure in the microporous regions is sufficient to inflate them. This kind of simulations enables the assessment of the rate of hydrocarbon recovery in shale gas according to the spatial distribution of the kerogen FIGURE 8 | Multi-scale simulations of the expulsion of methane stored in a microporous matrix into a soaked fracture. We investigated two scenarios: the fracture/matrix interface is non-wetting and the methane production is continuous, and the fracture/matrix interface is wetting and the gas is produced in a bubbly flow regime. The simulation setup is fully described in Soulaine et al. (2019).

Three-Phase Flow
Computational microfluidics for three-phase flow (e.g., gaswater-oil) are still scarce in the literature because the community focuses mostly on two-phase flow simulators for which the open challenges discussed in previous sections are still limiting the predictive capabilities. Such a tool is crucial to understand fully the underlying mechanisms of non-aqueous phase liquid (e.g., hydrocarbon, chlorinated solvents) trapped in the unsaturated zone by capillary effects and to propose innovative remediation strategies. In the literature, Bayestehparvin et al. (2015) use a Volume-of-Fluid approach, Helland et al. (2017) and Mohammadmoradi and Kantzas (2017) use quasi-static displacement using a Level-Set method, Panchenko (2016) use Smooth Particle Hydrodynamics, Dinariev andEvseev (2016) use a phase-field method and van Kats and Egberts (1999) and Jiang and Tsuji (2017) use a lattice-Boltzmann method to investigate three-phase Navier-Stokes flow at pore scale. The development of efficient three-phase flow simulators at the pore-scale is still under development and microfluidic experiments will help to assess the robustness of the computational models.

Reactive Transport Modeling at the Pore-Scale
Most of computational microfluidics developments so far have been devoted to solve the Navier-Stokes equations under single (Spanne et al., 1994;Bijeljic et al., 2013;Guibert et al., 2015; and two-phase flow conditions (Horgue et al., 2013;Raeini et al., 2014;Graveleau et al., 2017;Maes and Soulaine, 2018a;Pavuluri et al., 2020). Flow in geological porous media, however, has the distinctive feature to interact chemically with the solid walls. The consideration of surface reactions in computational microfluidics requires inclusion of comprehensive reaction networks along with flow and transport, an approach known as Reactive Transport Modeling (Steefel et al., 2005). Despite the growing investment in the development of RTM at the pore-scale-pioneer simulators date back to the late 90s (Békri et al., 1995(Békri et al., , 1997)-the field is still in its infancy. The main challenge consists in moving the fluid/solid boundary with respect to chemical reactions at the mineral surfaces. A large variety of methodologies has been proposed by the community to solve this problem and a comprehensive review can be found in Molins et al. (2020). They include Level-Set (Molins et al., 2014(Molins et al., , 2017, or phase-field Meakin, 2008, 2011) techniques in a Eulerian grid, LBM (Kang et al., 2003;Szymczak and Ladd, 2009;Prasianakis et al., 2013;Chen et al., 2014b;Yoon et al., 2015), SPH (Tartakovsky et al., 2007), and Arbitrary Lagrangian-Eulerian (Oltéan et al., 2013;Starchenko et al., 2016) frameworks. Another family of methods for displacing the mineral interface uses the micro-continuum approach Soulaine et al., 2017).

Modeling of Dissolution Processes
The modeling of pore-scale dissolution is the most advanced of the reactive transport processes using computational microfluidics. In this case, the solid mineral dissolves into the fluid phase according to thermodynamics and chemical reactions at the solid surface. High-resolution data of well-controlled microfluidic experiments have highlighted that the shape of a cylindrical crystal evolves into an elongated profile due to non-uniform distribution of the solute along the solid surface leading to faster dissolution rates upstream than downstream Dutka et al., 2020) as illustrated in Figure 9. The elongated profile of the crystal as it dissolved in the acid flowing solution illustrates well the complex coupling between streamlines, solute transport, and surface reaction. The stateof-the-art of computational microfluidics for reactive transport using different techniques including micro-continuum, Level-Set, LBM, moving grids with conformal mapping, and Vortex methods can reproduce the shape evolution of the dissolving crystal with great fidelity (Molins et al., 2020). This success leads to the strong conclusion that pore-scale technologies for moving fluid-solid boundaries along with geochemical reactions are now mature enough to be used as predictive tools. Current development uses computational microfluidics along with geochemical packages to model comprehensive reaction networks.
Computational microfluidics was used to investigate the dissolution of micro-model pore networks (Szymczak and Ladd, 2009;Kang et al., 2014;Molins et al., 2017;Soulaine et al., 2017). Five different dissolution regimes are observed according to the constant of reaction at the mineral surface and the injection mass flow rate as illustrated in Figure 10. For a given fluid, these two effects are characterized by the Damkohler number Da I (i.e., ratio of the characteristic timescale of transport and surface reaction) and the Peclet number Pe (i.e., ratio of advection and diffusion), respectively. These regimes are: compact dissolution (Pe < 10 −2 and Da I > 1), conical dissolution (1 > Pe > 10 −2 and Da I > 1), dominant wormhole (10 > Pe > 1 and Da I > 1), ramified wormholes (Pe > 10 and Da I > 1) and uniform dissolution (Da I < 1 and regardless the Peclet number). Because computational microfluidics provides a full mapping of the species and mineral evolution in the system, macroscopic properties can be obtained directly by averaging the simulation results, something that is barely achievable experimentally. For example, knowing the exact distribution of solute concentration including the local mass flux at the mineral surface  volume-averaged the pore-scale results and proposed a new law for quantifying the accessible reactive surface area, A, based on the transport and reactivity conditions. It reads where A 0 is the geometric surface area of the porous sample and n and m are model parameters. Similar relations that reduce the accessible surface area based on hydrodynamics conditions were used to reconcile experimental or numerical data at the larger scales (Wen and Li, 2017;Deng et al., 2018).

Modeling of Precipitation Processes
Although computational microfluidics for simulating precipitation processes are reported in the literature (Tartakovsky et al., 2007;Li et al., 2008;Yoon et al., 2012;Huber et al., 2014;Chen et al., 2015;Prasianakis et al., 2017;Ray et al., 2019), the field is still scarce and more development and verification are required to reach predictive capabilities. Pore-scale precipitation cannot be treated simply as the reverse problem of dissolution.
Indeed, the creation of the solid phase during the precipitation of minerals needs to consider the nucleation of mineral at the solid surface. Moreover, the development of depth-averaged model for precipitation is challenging because in microfluidic experiments it is possible that precipitate preferentially form either on the top and bottom plates or in the bulk (Beuvier et al., 2015;Yoon et al., 2019). Model verification will use the ability of microfluidics to design well-controlled experiments along with high-resolution imaging and chemical characterization (Poonoosamy et al., 2020).

Multi-Scale Reactive Transport Models
Hybrid-scale models have been proposed to describe reactive systems that include multiple characteristic length-scales for which some regions are described using pore-scale modeling while others are modeled with continuum approaches (Liu and Ortoleva, 1996;Liu et al., 1997). Two kinds of approaches solve hybrid-scale problems. On the one hand, the domain FIGURE 11 | Numerical simulation of tertiary low-salinity flooding in a micromodel. The colormap corresponds to the ionic strength (mol/l) in the domain. The non-aqueous phase appears in green (I = 0), the high salinity water appears in dark blue (I = 3.6) and the low-salinity water appears in light blue (I = 0.9). The simulation parameters are fully described in Maes and Geiger (2018).
decomposition technique solves different physics on separate domains-one for Darcy flow, another for Stokes flow-linked together through appropriate boundary conditions (Tartakovsky et al., 2008;Molins et al., 2019). On the other hand, microcontinuum models use a single set of partial differential equations throughout the computational domain regardless of the content of a grid block (Steefel et al., 2015;Soulaine and Tchelepi, 2016). The latter approach is particularly well-suited to capture the dynamic displacement of the interface between the porous and solid-free regions without involving complex re-meshing strategies. For example, micro-continuum models have been used successfully to simulate the formation and growth of wormholes in acidic environments (Ormond and Ortoleva, 2000;Golfier et al., 2002;Soulaine and Tchelepi, 2016;Faris et al., 2020).

Two-Phase Flow in Reactive Environments
The dynamics modeling of multiphase flow phenomena in reactive environments is quite complex because it involves tracking multiple interfaces that evolve as a function of the details of the reactive transport: on the one hand, the fluid/fluid interface is subject to surface tension forces, on the other hand, the fluid/solid interface moves with chemical reactions at the surface of the solid minerals. The modeling is complicated by the dynamics of the contact line or thin films along the solid walls, and by the transfer of mass across interfaces. Few research works have been devoted to the modeling of multiphase flow with solid boundaries that evolve with chemical reactions. Parmigiani et al. (2011) have proposed an LBM approach to study the evolution of capillary fingers in a porous medium that dynamically evolves with the melting of the solid phase. Later, LBM has been extended to multicomponent multiphase fluid flow and applied to reactive transport with dissolution and precipitation (Chen et al., 2013(Chen et al., , 2015. LBM-based approaches, however, are unstable for highdensity ratios between the liquid and the gas (Chen et al., 2014a), which, therefore, limits the predictive aspects of such modeling frameworks. Alternatively, multicomponent multiphase reactive transport can be modeled using the Volume Of Fluid method to track the fluid/fluid interface and the C-CST approach to model the transport of reactive species in the system, including the effect of speciation and surface reactions . Maes and Geiger (2018) used this model to track pH, solid surface potential and contact angle alteration during injection of lowsalinity water in micromodels and were able to observe and characterize for the first time the effect of secondary and tertiary low-salinity flooding at the pore-scale (Figure 11). Moreover, Soulaine et al. (2018) proposes new computational microfluidics FIGURE 12 | Comparison between experimental and computational microfluidics. A calcite crystal is dissolved in a flow of hydrochloric acid at atmospheric conditions. Carbon dioxide gas bubbles nucleate at the mineral surface, grow, coalesce, and eventually detach. The micro-continuum model for two-phase flow in reactive environments can reproduce accurately the microfluidic experiments. The microfluidic experiment and simulation are fully described in Soulaine et al. (2018). to simulate mineral dissolution in unsaturated porous media. The model combines the reactive micro-continuum approach developed in Soulaine et al. (2017) with the VOF and C-CST methods. Chemical reactions and wettability conditions at the mineral surface are described as immersed boundaries that evolve with the dissolution of the solid phase. Microfluidic experiments that consist of a calcite crystal dissolving in a microchannel after the injection of an aqueous solution containing up to 1% of hydrochloric acid were used to verify the predictive aspect of the model (see Figure 12).
The model has been used to investigate the effect of carbon dioxide production during the acidizing of carbonate formations. Results show that under single-phase flow conditions, i.e., when the dissolution time of CO 2 into the aqueous phase is instantaneous, the solid structure dissolves forming wormholes. However, if the characteristic time scale of CO 2 dissolving in brine is much longer compared to the rate of gas production, then CO 2 bubbles grow, coalesce and form flow barriers that limit the transport of the acid in the domain (see Figure 13). These flow barriers prevent the emergence of wormholes and limit significantly the overall dissolution rate. This finding is of great interest in acid stimulation processes for which a gas phase can be produced (Prutton and Savage, 1945;Thompson and Gdanski, 1993). Hence, the presence of a second fluid phase under the same flow conditions may lead to very different dissolution regimes, and the common behavior diagrams that characterize the dissolution pattern based on the Peclet and Damkohler numbers (Golfier et al., 2002) have to be complemented, at a minimum, by a third dimension that quantifies the solubility of the CO 2 in the aqueous phase.

CONCLUSION AND PERSPECTIVES
Computational microfluidics for geosciences offers an appealing framework to investigate the coupled processes in geological porous media. The integrated usage of experimental and computational microfluidics is, therefore, a powerful approach to decipher the complex mechanisms that occur in the soils and the subsurface and to bring new insights into the open challenges associated with the continuum-scale modeling of flow and transport in geological systems. Microfluidic experiments provide valuable benchmark data sets, leading to more accurate and robust numerical solvers, while numerical simulations augment microfluidic experiments by providing high-resolution mapping of pressure and velocity profiles, solute concentration, and mineral distribution that are not easily measurable experimentally. Modern approaches can simulate pore-scale events over a wide range of conditions including multiphase flow and the microstructure evolution due to precipitation and dissolution. Combining experimental and computational microfluidics have led to new insights into complex subsurface processes, such as CO 2 geological storage, low salinity flooding and acid stimulation. The paper presents a snapshot of the current capabilities and discuss the validity of depth-averaged models.
The future of this emerging scientific discipline is guaranteed by the continuous improvements of numerical algorithms, modeling, and high-performance computing. The next generations of computational microfluidics will be able to simulate the transport of particles including the micro-dynamics of colloids, their aggregation and stability, the attachment/detachment from the solid walls, and pore-clogging mechanisms according to pH ad salinity conditions. This will necessitate handling the strong coupling between hydro-electrochemical effects. Important efforts will be dedicated to the consideration of biological processes including the growth of biofilms.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

ACKNOWLEDGMENTS
We acknowledge the first reviewer for suggesting timedependent dispersion. The research leading to these results has received support from the French Agency for Research (Agence Nationale de la Recherche, ANR) through the Equipex Planex ANR-11-EQPX-36, the Labex Voltaire ANR-10-LABX-100-01, the grant CATCH ANR-18-CE05-0035, and through the FraMatI project under contract ANR-19-CE05-0002. It has also received financial support from the CNRS through the MITI interdisciplinary programs, and from the UK EPSRC funded project on Direct Numerical Simulation for Additive Manufacturing in Porous Media (grant reference EP/P031307/1). Finally, SR and JM would like to thank ALLIANCE (Partenariat Hubert Curien) for the grant that enable the scholar exchange of scientists.