2.1 Modelling approach

Figure 1.

Schematic of the micro-continuum approach extended to multiple fluid phases. (a) The full Navier–Stokes representation of the physics (solid grains in black, liquid in blue and gas in orange); (b) map of the solid volume fraction after filtering; (c) map of the liquid saturation after filtering.

Numerical simulation of multiphase dissolution at the pore scale is challenging because it involves multiple interfaces (gas/liquid, fluid/solid) that evolve as a function of the details of the reactive transport. Instead of the full Navier–Stokes simulation approach (see figure 1 a), we employ a micro-continuum formulation. That is, all the information regarding the solid structure at a scale smaller than the computational cells (control volumes) is filtered and modelled as a fluid/solid aggregate. In the DBS approach, all the physical variables of the solid–fluids system are averaged over the cell’s control volume, $V$ . Cells that contain pore space (fluids) and those that contain solid are differentiated using the scalar field $\unicode[STIX]{x1D700}$ , which corresponds to the volume fraction of the pore space in the control volumes. It is a field with different values in the computational domain that range from 0 to 1. This filtering approach is illustrated in figure 1(b). The free zone is characterized by $\unicode[STIX]{x1D700}=1$ , and the flow is governed by the two-phase Navier–Stokes equations. For intermediate values, $0<\unicode[STIX]{x1D700}<1$ , there are solid ‘obstacles’ in the control volume, and the physics of flow and transport in porous media is applied. When $\unicode[STIX]{x1D700}=0$ , the control volume contains solid only, and there is no flow. We assume that the cells that are fully occupied by the solid always contain a tiny porosity, $\unicode[STIX]{x1D700}=0\equiv 0.001$ , and that allows for a single-field DBS formulation regardless of the make-up of a computational cell. Different from other micro-continuum models proposed to simulate core-scale phenomena (Liu et al. Reference Liu, Ormond, Bartko, Ying and Ortoleva1997; Ormond & Ortoleva Reference Ormond and Ortoleva2000; Golfier et al. Reference Golfier, Zarcone, Bazin, Lenormand, Lasseux and Quintard2002; Scheibe et al. Reference Scheibe, Perkins, Richmond, McKinley, Romero-Gomez, Oostrom, Wietsma, Serkowski and Zachara2015; Soulaine & Tchelepi Reference Soulaine and Tchelepi2016a ), this work uses a DBS formulation to simulate pore-scale physics. Thus, only the bounding values $\unicode[STIX]{x1D700}=1$ and $\unicode[STIX]{x1D700}=0$ are used to map the entire structure at the pore scale. The intermediate values of the pore space (or the solid) describe the interface between the fluid and the solid. This approach has been used successfully to track the displacement of a fluid/solid interface during dissolution in the presence of a single-phase fluid (Soulaine et al. Reference Soulaine, Roman, Kovscek and Tchelepi2017).

The micro-continuum approach is a physically based model, whereby the governing equations are locally averaged. This approach inherits all the theoretical developments of the physics of flow and transport in porous media using the method of volume averaging (Whitaker Reference Whitaker1999). The micro-continuum approach relies on the fact that for each variable of the system, a single conservation equation holds for both the porous and pore (fluids) domains: in the free zone (pore space), the equations tend asymptotically to the standard equations of fluid dynamics, while in the porous regions they reduce to the partial differential equations that govern flow and transport in porous media. For example, for single phase flow, the Darcy–Brinkman–Stokes (DBS) equation (Brinkman Reference Brinkman1947) can be used to simulate flow in a channel surrounded by a porous matrix (Neale & Nader Reference Neale and Nader1974; Soulaine et al. Reference Soulaine, Gjetvaj, Garing, Roman, Russian, Gouze and Tchelepi2016; Soulaine & Tchelepi Reference Soulaine and Tchelepi2016a ). In the channel, the drag force vanishes since there is no solid, and the DBS formalism turns to the Stokes equation. In the ‘porous’ matrix, the drag force is dominant with respect to the viscous term, and the momentum balance tends toward Darcy’s law. In this work, the solid phase is described by a low-permeability low-porosity matrix and the Darcy term is used to penalize the flow in this region.

When dissolution occurs in a control volume, the volume fraction of solid, $\unicode[STIX]{x1D700}_{s}$ , decreases. Consequently, the local porosity, $\unicode[STIX]{x1D700}$ , increases. This is because the relation

(2.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D700}_{s}=1-\unicode[STIX]{x1D700} & & \displaystyle\end{eqnarray}$$

is always valid. The variation of the solid volume fraction is used to track the evolution of the pore topology.

In the micro-continuum approach extended to multiphase system, the pore space is occupied by two fluids, namely gas and liquid denoted $g$ and $l$ , respectively. The distribution of fluids in the domain is described by the volume fraction of gas, $\unicode[STIX]{x1D700}_{g}$ , and the volume fraction of liquid, $\unicode[STIX]{x1D700}_{l}$ , in every control volumes (see figure 1 c). Note that by construction, the relations $\unicode[STIX]{x1D700}_{s}+\unicode[STIX]{x1D700}_{g}+\unicode[STIX]{x1D700}_{l}=1$ and $\unicode[STIX]{x1D700}=\unicode[STIX]{x1D700}_{g}+\unicode[STIX]{x1D700}_{l}$ are always valid. These quantities evolve with chemical reactions at the mineral surface and with the hydrodynamics of the system. The distribution of the fluids in the domain is tracked using the concept of saturation, $S_{f}$ with $f=l,g$ , that represents the ratio of the volume occupied by the $f$ -phase over the volume of the pore space,

(2.2) $$\begin{eqnarray}\displaystyle S_{f}\equiv \frac{\unicode[STIX]{x1D700}_{f}}{\unicode[STIX]{x1D700}}\quad \text{with}~f=l,g. & & \displaystyle\end{eqnarray}$$

It has the advantage to be normalized: $S_{f}$ varies in the range 0–1. The relation,

(2.3) $$\begin{eqnarray}\displaystyle S_{l}+S_{g}=1, & & \displaystyle\end{eqnarray}$$

is always satisfied. Thus, we only need to compute the liquid (or gas) saturation. In the pore space, i.e. $\unicode[STIX]{x1D700}=1$ , saturation is similar to the colour function used in the volume-of-fluid (VOF) approach to track fluid/fluid interfaces. Namely, in the control volumes occupied by liquid only, $S_{l}=1$ ; in regions occupied by gas only, $S_{l}=0$ ; and for gas/liquid interface, $0<S_{l}<1$ .

As for the volume fractions and the saturation, all the physical variables of the system can be defined as averaged quantities over the control volume. For each phase ( $f=l,g$ ), we define the phase averaged pressure $\bar{p}_{f}=1/V_{f}\int _{V_{f}}p_{f}\,\text{d}V$ , the phase averaged velocity $\bar{\mathsf{v}}_{f}=1/V_{f}\int _{V_{f}}\mathsf{v}_{f}\,\text{d}V$ and the phase averaged mass fraction of a chemical species $A$ in the mixture, $\bar{\unicode[STIX]{x1D714}}_{f,A}=1/V_{f}\int _{V_{f}}\unicode[STIX]{x1D714}_{f,A}\,\text{d}V$ where $V_{f}$ is the volume occupied by the fluid $f=l,g$ in the control volume.

In the pore-scale multiphase micro-continuum model, the problem is formulated in terms of single field (global variables) regardless of the content of cell. Hence, a single momentum equation is used to solve for the mixture velocity, $\bar{\mathsf{v}}$ , defined as,

(2.4) $$\begin{eqnarray}\displaystyle \bar{\mathsf{v}}=\unicode[STIX]{x1D700}_{l}\bar{\mathsf{v}}_{l}+\unicode[STIX]{x1D700}_{g}\bar{\mathsf{v}}_{g}=\unicode[STIX]{x1D700}(S_{l}\bar{\mathsf{v}}_{l}+(1-S_{l})\bar{\mathsf{v}}_{g}). & & \displaystyle\end{eqnarray}$$

On the one hand, in the free zone, $\bar{\mathsf{v}}$ corresponds to the velocity field used in the VOF method: if $S_{l}=1$ , it denotes the velocity of the liquid, if $S_{l}=0$ , it depicts the velocity of the gas phase and if $0<S_{l}<1$ , it corresponds to the velocity of the gas/liquid interface. In the solid–fluids aggregates domain, $\bar{\mathsf{v}}$ corresponds to Darcy’s mixture velocity (Wang & Beckermann Reference Wang and Beckermann1993).

The single-field pressure is defined as

(2.5) $$\begin{eqnarray}\displaystyle \bar{p}=S_{l}\bar{p}_{l}+(1-S_{l})\bar{p}_{g}, & & \displaystyle\end{eqnarray}$$

and the concentration field for species $A$ as,

(2.6) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D714}}_{A}=S_{l}\bar{\unicode[STIX]{x1D714}}_{l,A}+(1-S_{l})\bar{\unicode[STIX]{x1D714}}_{g,A}. & & \displaystyle\end{eqnarray}$$

Next, we introduce the governing equations for the evolution of the solid volume fraction $\unicode[STIX]{x1D700}_{s}$ , the saturation $S_{l}$ , the velocity field $\bar{\mathsf{v}}$ , the pressure $\bar{p}$ and the concentration field $\bar{\unicode[STIX]{x1D714}}_{A}$ .

2.2 Mass transfer

In this system, the pore geometry, and hence the porosity, permeability and flow behaviour evolve with the chemical reactions that take place at the solid boundaries. For acidization processes involving hydrochloric acid (HCl) in carbonate formations, mineral dissolution is dictated by the following reactions at the surface of the calcite crystal:

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle \text{CaCO}_{3(s)}+2\text{HCl}_{(aq)}\rightarrow \text{CaCl}_{2(aq)}+\text{CO}_{2(aq)}+\text{H}_{2}\text{O}_{(l)} & \displaystyle\end{eqnarray}$$

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle \text{CO}_{2(aq)}\rightarrow \text{CO}_{2(g)}. & \displaystyle\end{eqnarray}$$

The $\text{CO}_{2}$ forms in the aqueous phase and then nucleates as a gas phase on the surface of the grain. We neglect the $\text{CO}_{2(aq)}$ in the mass balance and we consider that the $\text{CO}_{2}$ from the aqueous phase instantaneously goes to the gas phase. Moreover, we assume that the characteristic time of the $\text{CO}_{2}$ mass transfer from the gas to the liquid phase is much longer compared with the rate of $\text{CO}_{2}$ production from the dissolution of the rock. Hence, HCl is the only species that needs to be tracked to compute the reaction rate at the solid surface, ${\dot{m}}_{l,HCl}$ . Then, it follows that the rates of the phase change for the solid, gas and liquid phases are directly described by the reaction rate of the acid (HCl). Consequently, we can write:

(2.9) $$\begin{eqnarray}\displaystyle \dot{m_{i}}=\unicode[STIX]{x1D6FE}_{i}{\dot{m}}_{l,HCl}\quad \text{with}~i=s,l,g, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}_{i}$ is based on molar mass and stoichiometric coefficients: $\unicode[STIX]{x1D6FE}_{s}=M_{CaCO_{3}}/2M_{HCl}$ , $\unicode[STIX]{x1D6FE}_{g}=-M_{CO_{2}}/2M_{HCl}$ and $\unicode[STIX]{x1D6FE}_{l}=(2M_{HCl}-M_{CaCl_{2}}-M_{H_{2}O})/2M_{HCl}$ .

2.3 Mass balance equation

The volume fractions of the three phases (solid, liquid, gas) that occupy the domain evolve in time as a function of the hydrodynamics and the chemical reactions at the solid surface. The volume fraction of solid in a control volume is an unknown (degree-of-freedom) of the system, and its evolution is governed by the mass balance equation:

(2.10) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D700}_{s}\unicode[STIX]{x1D70C}_{s}}{\unicode[STIX]{x2202}t}={\dot{m}}_{s}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}_{s}$ is the density of the solid, and ${\dot{m}}_{s}$ , is the rate of dissolution of solid into gas and liquid. The latter has a non-zero value only at the fluid/solid interface. The dissolved mass goes into the two fluid phases. Hence, the conservation laws for the liquid and gas are:

(2.11) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D700}_{l}\unicode[STIX]{x1D70C}_{l}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D700}_{l}\unicode[STIX]{x1D70C}_{l}\bar{\mathsf{v}}_{l})=\dot{m_{l}}, & & \displaystyle\end{eqnarray}$$

and,

(2.12) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D700}_{g}\unicode[STIX]{x1D70C}_{g}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D700}_{g}\unicode[STIX]{x1D70C}_{g}\bar{\mathsf{v}}_{g})={\dot{m}}_{g}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}_{l}$ and $\unicode[STIX]{x1D70C}_{g}$ stand for the density of the liquid and gas, respectively. The source terms $\dot{m_{l}}$ and $\dot{m_{g}}$ denote the rate of production of liquid and gas, respectively.

In our DBS framework, we recast this set of equations into a global mass balance equation and a saturation (transport) equation. The global mass conservation equation is obtained by summing (2.10)–(2.12). Assuming that both the liquid and the gas phases are incompressible, we have

(2.13) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\bar{\mathsf{v}}={\dot{m}}_{l,HCl}\left(\frac{\unicode[STIX]{x1D6FE}_{l}}{\unicode[STIX]{x1D70C}_{l}}+\frac{\unicode[STIX]{x1D6FE}_{g}}{\unicode[STIX]{x1D70C}_{g}}+\frac{\unicode[STIX]{x1D6FE}_{s}}{\unicode[STIX]{x1D70C}_{s}}\right). & & \displaystyle\end{eqnarray}$$

To track the liquid saturation, the mass conservation equation for the liquid, equation (2.11) is transformed into (see appendix A for the derivation)

(2.14) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D700}S_{l}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\bar{\mathsf{v}}S_{l})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D700}S_{l}(1-S_{l})\bar{\mathsf{v}}_{r})=\frac{\dot{m_{l}}}{\unicode[STIX]{x1D70C}_{l}}, & & \displaystyle\end{eqnarray}$$

where $\bar{\mathsf{v}}_{r}=(\bar{\mathsf{v}}_{l}-\bar{\mathsf{v}}_{g})$ is a compression velocity. In the free zone, $\unicode[STIX]{x1D700}=1$ , and this equation turns into the phase indicator equation used in the algebraic VOF solver implemented in OpenFOAM^® (Rusche Reference Rusche2003). The last term on the left-hand side is non-zero only in the gas/liquid transition region. It is used to sharpen the interface. The value of the relative velocity, $\bar{\mathsf{v}}_{r}$ , is estimated from the maximum velocity magnitude in this interfacial area (Rusche Reference Rusche2003).

2.4 Momentum balance equation and adhesion to the solid

The flow problem under consideration involves the evolution of multiple interfaces. The delineation between the free zone, i.e.  $\unicode[STIX]{x1D700}=1$ , and the solid (modelled as a porous medium) changes as the solid dissolves. Thus, the topology of the void space changes substantially. In the free zone, the gas and liquid evolve as two immiscible fluids separated by an interface with surface tension. When the three phases meet, the fluid/fluid interface forms a contact angle with the solid surface, which depends on the wettability and the hydrodynamic conditions.

To model the momentum balance equation, we adopt a combination of the Darcy–Brinkman–Stokes equation, which has proven to be quite effective for solving single-phase dissolution problems (Soulaine et al. Reference Soulaine, Roman, Kovscek and Tchelepi2017), and the classic VOF method used to track gas–liquid interfaces on a fixed Eulerian grid (Hirt & Nichols Reference Hirt and Nichols1981). In this model, a single equation holds for the liquid, gas and the porous region. The momentum equation is expressed as:

(2.15) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D700}}\left(\frac{\unicode[STIX]{x2202}\bar{\mathsf{v}}}{\unicode[STIX]{x2202}t}+\frac{\bar{\mathsf{v}}}{\unicode[STIX]{x1D700}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{\mathsf{v}}\right)=-\unicode[STIX]{x1D735}\bar{p}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D700}}(\unicode[STIX]{x1D735}\bar{\mathsf{v}}+^{t}\unicode[STIX]{x1D735}\bar{\mathsf{v}})\right)+\boldsymbol{F}_{c}-\unicode[STIX]{x1D707}k^{-1}\bar{\mathsf{v}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}$ is the density, $\unicode[STIX]{x1D707}$ is the dynamic viscosity, $\boldsymbol{F}_{c}$ is the surface-tension force due to curvature of the interface. The last term of the right-hand side, $\unicode[STIX]{x1D707}k^{-1}\bar{\mathsf{v}}$ , represents the momentum exchange term between the fluids and the solid phase. The density and viscosity depend on the saturation through the relations, $\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{l}S_{l}+\unicode[STIX]{x1D70C}_{g}(1-S_{l})$ and $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{l}S_{l}+\unicode[STIX]{x1D707}_{g}(1-S_{l})$ . This model is similar to the penalization technique applied to the VOF method proposed by Horgue, Prat & Quintard (Reference Horgue, Prat and Quintard2014). Here, there is the additional complexity that the solid structure evolves in time. Horgue et al. (Reference Horgue, Prat and Quintard2014) solve a momentum equation similar to (2.15) and use the flow resistance term, $\unicode[STIX]{x1D707}k^{-1}\bar{\mathsf{v}}$ , to force the velocity values to near zero in regions identified as solid, and therefore enforce a no-slip boundary condition at the immersed boundaries. The wall adhesion condition is ensured by local modification of the surface tension force, $\boldsymbol{F}_{c}$ , at the wall boundary to enforce the value of the contact angle formed by the gas–liquid interface and the solid surface.

As for the Darcy–Brinkman–Stokes model applied to the simulation of flow in large channels surrounded by porous media (Neale & Nader Reference Neale and Nader1974), the resistance term, $\unicode[STIX]{x1D707}k^{-1}\bar{\mathsf{v}}$ , vanishes in the free zone ( $\unicode[STIX]{x1D700}=1$ ), so that (2.15) reduces to the classic VOF momentum equation, and is dominant in the porous regions so that the momentum equation tends toward Darcy’s law in these regions. Because the fluids/solid interface moves with the dissolution process, the value of the local permeability evolves with the local porosity, $\unicode[STIX]{x1D700}$ .

We use a Kozeny–Carman relationship to compute $k$ from the local porosity (Voller Reference Voller2009):

(2.16) $$\begin{eqnarray}\displaystyle k^{-1}=k_{0}^{-1}\frac{(1-\unicode[STIX]{x1D700})^{2}}{\unicode[STIX]{x1D700}^{3}}, & & \displaystyle\end{eqnarray}$$

with $\sqrt{k_{0}}$ the characteristic length of the micro-structure. With this formulation, $k^{-1}=0~\text{m}^{-2}$ in the free zone and tends towards a very large value when the porosity is close to zero. In the latter case, the velocity in the porous region drops down to values near zero, which mimics a solid with no flow as in the fictitious domain approach (Angot, Bruneau & Fabrie Reference Angot, Bruneau and Fabrie1999; Khadra et al. Reference Khadra, Angot, Parneix and Caltagirone2000; Horgue et al. Reference Horgue, Prat and Quintard2014; Soulaine & Tchelepi Reference Soulaine and Tchelepi2016a ). With dissolution, the cell porosity evolves and that changes the local permeability until all the solid contained in a computational cell is fully consumed, and the model switches automatically from Darcy’s law to the VOF-based Navier–Stokes representation.

Figure 2.

Illustration of the wall adhesion condition. The normal to the gas/liquid interface, $\boldsymbol{n}_{lg}$ , forms an angle $\unicode[STIX]{x1D703}$ with the normal to the solid wall, $\boldsymbol{n}_{wall}$ . These two normal vectors are computed from the gradient of the volume fraction fields.

The capillary effects due to the surface tension between the two fluids and the wettability with respect to the solid minerals are described by a unique volumetric source term, $\boldsymbol{F}_{c}$ , defined throughout the computational grid. The nature of this term, however, is different according to the cell porosity. In the free zone, the two fluid phases are immiscible, and the curvature of the interface leads to a discontinuity in the stress to balance the pressure gradient according to Young–Laplace law. At the fluid/solid boundary, the liquid–gas interface meets the solid surface forming a contact angle that reflects the wettability of the system. To model $\boldsymbol{F}_{c}$ with a unique formulation that is valid everywhere in the computational domain, we adopt the strategy proposed by Horgue et al. (Reference Horgue, Prat and Quintard2014). Regardless of the region of interest, we have

(2.17) $$\begin{eqnarray}\displaystyle \boldsymbol{F}_{c}=\unicode[STIX]{x1D700}\unicode[STIX]{x1D70E}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\hat{\boldsymbol{n}}_{lg})\unicode[STIX]{x1D735}S_{l}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}$ is the surface tension between the gas and the liquid, and $\hat{\boldsymbol{n}}_{lg}$ is the normal to the gas/liquid interface locally modified to account for the wall adhesion condition. In the clear-fluid region, $\hat{\boldsymbol{n}}_{lg}$ corresponds to the normal to the gas–liquid interface $\boldsymbol{n}_{lg}=\unicode[STIX]{x1D735}S_{l}/\Vert \unicode[STIX]{x1D735}S_{l}\Vert$ and we have $\boldsymbol{F}_{c}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}S_{l}/\Vert \unicode[STIX]{x1D735}S_{l}\Vert )\unicode[STIX]{x1D735}S_{l}$ , which corresponds to the continuum surface force (CSF) introduced by Brackbill, Kothe & Zemach (Reference Brackbill, Kothe and Zemach1992). At the solid boundaries, $\hat{\boldsymbol{n}}_{lg}$ is modified to enforce the wall adhesion condition:

(2.18) $$\begin{eqnarray}\displaystyle \hat{\boldsymbol{n}}_{lg}=\cos \unicode[STIX]{x1D703}\boldsymbol{n}_{wall}+\sin \unicode[STIX]{x1D703}\boldsymbol{t}_{wall}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{n}_{wall}$ and $\boldsymbol{t}_{wall}$ are the normal and tangent vectors to the solid surface, respectively, and $\unicode[STIX]{x1D703}$ is the contact angle (see figure 2).

We use the correction strategy proposed by Horgue et al. (Reference Horgue, Prat and Quintard2014): after solving the saturation equation. The normal, $\boldsymbol{n}_{lg}$ , to the gas/liquid interface does not satisfy (2.18) but forms an angle $\unicode[STIX]{x1D703}_{I}$ with the vector normal to the wall, $\boldsymbol{n}_{wall}$ . That is, $\unicode[STIX]{x1D703}_{I}=\cos ^{-1}(\boldsymbol{n}_{wall}\boldsymbol{\cdot }\boldsymbol{n}_{lg})$ . The locally modified normal vector, $\hat{\boldsymbol{n}}_{lg}$ , used to compute the capillary forces in (2.17) – and defined over the entire computational domain – is corrected to enforce the correct value of $\unicode[STIX]{x1D703}$ by using the formula:

(2.19) $$\begin{eqnarray}\displaystyle \hat{\boldsymbol{n}}_{lg}=\frac{\cos \unicode[STIX]{x1D703}-\cos \unicode[STIX]{x1D703}_{I}\cos (\unicode[STIX]{x1D703}_{I}-\unicode[STIX]{x1D703})}{1-\cos ^{2}\unicode[STIX]{x1D703}_{I}}\boldsymbol{n}_{wall}+\frac{\cos (\unicode[STIX]{x1D703}_{I}-\unicode[STIX]{x1D703})-\cos \unicode[STIX]{x1D703}_{I}\cos \unicode[STIX]{x1D703}}{1-\cos ^{2}\unicode[STIX]{x1D703}_{I}}\boldsymbol{n}_{lg}. & & \displaystyle\end{eqnarray}$$

The derivation of this formula can be found in appendix B. The normal vectors to the solid wall and to the gas/liquid interface are computed using the volume fractions of the solid and liquid, respectively, with $\boldsymbol{n}_{wall}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D700}_{s}/\Vert \unicode[STIX]{x1D735}\unicode[STIX]{x1D700}_{s}\Vert$ and $\boldsymbol{n}_{lg}=\unicode[STIX]{x1D735}S_{l}/\Vert \unicode[STIX]{x1D735}S_{l}\Vert$ . With (2.19), both Brackbill’s CSF and the wall adhesion condition are satisfied using a unique formulation. In the free region, $\boldsymbol{n}_{wall}=\mathbf{0}$ and we have $\hat{\boldsymbol{n}}_{lg}/\Vert \hat{\boldsymbol{n}}_{lg}\Vert =\boldsymbol{n}_{lg}/\Vert \boldsymbol{n}_{lg}\Vert$ , and then $\boldsymbol{F}_{c}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}S_{l}/\Vert \unicode[STIX]{x1D735}S_{l}\Vert )\unicode[STIX]{x1D735}S_{l}$ . At the solid boundary, equation (2.19) automatically enforces the wall adhesion condition, $\hat{\boldsymbol{n}}_{lg}=\cos \unicode[STIX]{x1D703}\boldsymbol{n}_{wall}+\sin \unicode[STIX]{x1D703}\boldsymbol{t}_{wall}$ . Refer to Horgue et al. (Reference Horgue, Prat and Quintard2014) for details regarding the implementation and validation of this technique. In the porous regions, the surface-tension forces under confinement corresponds to the gradient of the capillary pressure, $\boldsymbol{F}_{c}=\unicode[STIX]{x1D735}p_{c}=(\unicode[STIX]{x2202}p_{c}/\unicode[STIX]{x2202}S_{l})\unicode[STIX]{x1D735}S_{l}$ . In the pore-scale simulations presented in this paper, the porous regions are initially saturated by the liquid at the thermodynamic equilibrium with the solid. Hence, there is no chemical reaction except at the boundaries of the porous regions. Therefore, in the matrix, we always have $\unicode[STIX]{x1D735}S_{l}=\mathbf{0}$ , and then $\boldsymbol{F}_{c}=\mathbf{0}$ .

2.5 Concentration equation

In a multiphase system, the chemical components are present in both fluid phases. For instance, there is $\text{CO}_{2}$ in the gas phase and dissolved $\text{CO}_{2}$ in the aqueous phase. The component concentration in each phase is governed by the advection–diffusion conservation equations coupled with the thermodynamic equilibrium condition at the gas/liquid interface. This latter condition is usually described using a partitioning relation, such as Henry’s laws (see figure 3):

(2.20) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{l,A}=H_{A}\unicode[STIX]{x1D714}_{g,A}, & & \displaystyle\end{eqnarray}$$

where $H_{A}$ is the partitioning coefficient or Henry’s constant. When this condition is not fulfilled, there is a mass transfer from one phase to one another in order to reach thermodynamic equilibrium.

Figure 3.

Illustration of the mass transfer between two immiscible phases. If Henry’s law is not satisfied at the gas/liquid interface, then there is a mass transfer across the interface in order to reach thermodynamic equilibrium.

In the micro-continuum approach, a single conservation equation governs the global concentration field, $\bar{\unicode[STIX]{x1D714}}_{A}$ . Two challenges arise to transport $\bar{\unicode[STIX]{x1D714}}_{A}$ using a single equation in a domain that contains two fluids. In the free zone, this equation has to be compatible with the VOF formulation to transport the species in a two-fluid system, where the thermodynamic equilibrium condition at the interface between the fluids, equation (2.20), is always satisfied. This equation should also capture the chemical reaction at the minerals surface. The model we propose for the evolution of $\bar{\unicode[STIX]{x1D714}}_{A}$ is a combination of the method proposed by Haroun, Legendre & Raynal (Reference Haroun, Legendre and Raynal2010) to transport a concentration field in a two-fluid system and the immersed boundary technique for micro-continuum models proposed by Soulaine et al. (Reference Soulaine, Roman, Kovscek and Tchelepi2017). Specifically, we write

(2.21) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D700}\unicode[STIX]{x1D70C}\bar{\unicode[STIX]{x1D714}}_{A}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}\bar{\mathsf{v}}\bar{\unicode[STIX]{x1D714}}_{A})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D700}\unicode[STIX]{x1D70C}D_{A}^{\ast }(\unicode[STIX]{x1D735}\bar{\unicode[STIX]{x1D714}}_{A}+\unicode[STIX]{x1D731}_{A}))-{\dot{m}}_{l,A}, & & \displaystyle\end{eqnarray}$$

where $D_{A}^{\ast }$ is the effective diffusion coefficient, $\unicode[STIX]{x1D731}_{A}$ is the continuum species transfer (CST) flux (Haroun et al. Reference Haroun, Legendre and Raynal2010; Marschall et al. Reference Marschall, Hinterberger, Schüler, Habla and Hinrichsen2012; Graveleau et al. Reference Graveleau, Soulaine and Tchelepi2017; Maes & Geiger Reference Maes and Geiger2017; Maes & Soulaine Reference Maes and Soulaine2018) and ${\dot{m}}_{l,A}$ is a reaction rate at the mineral surface. This equation is valid everywhere in the computational domain, regardless the cells content.

The effective diffusion coefficient, $D_{A}^{\ast }$ , is:

(2.22) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D700}D_{A}^{\ast }=\unicode[STIX]{x1D700}^{2}\frac{D_{l,A}D_{g,A}}{S_{l}D_{g,A}+(1-S_{l})D_{l,A}}, & & \displaystyle\end{eqnarray}$$

where $D_{g,A}$ and $D_{l,A}$ are the molecular diffusivity of the component $A$ in the gas and liquid phase, respectively. Haroun et al. (Reference Haroun, Legendre and Raynal2010) and Deising, Marschall & Bothe (Reference Deising, Marschall and Bothe2016) have demonstrated that this harmonic formulation is more robust than a simple mixing rule. In the free zone, $D_{A}^{\ast }=D_{g,A}$ in the gas phase and $D_{A}^{\ast }=D_{l,A}$ in the liquid phase. Similar to the single-phase micro-continuum model of Soulaine et al. (Reference Soulaine, Roman, Kovscek and Tchelepi2017), the exponent 2 allows the value of diffusion in the solid regions to drop to near zero.

The additional flux in (2.21), $\unicode[STIX]{x1D731}_{A}$ , describes the thermodynamic equilibrium at the gas/liquid interface. It transforms the solubility condition, equation (2.20), into a volumetric term under the framework of the VOF formulation (Haroun et al. Reference Haroun, Legendre and Raynal2010). We have:

(2.23) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D731}_{A}=\frac{\bar{\unicode[STIX]{x1D714}}_{A}(\unicode[STIX]{x1D70C}_{g}H_{A}-\unicode[STIX]{x1D70C}_{l})}{\unicode[STIX]{x1D70C}_{l}H_{A}S_{l}+\unicode[STIX]{x1D70C}_{g}(1-S_{l})}\unicode[STIX]{x1D735}S_{l}. & & \displaystyle\end{eqnarray}$$

The CST term is responsible for the concentration jump at the interface, while the flux continuity across the interface is always insured. This approach is similar of the continuum surface force (CSF) used for the modelling of the surface tension between two fluids (Brackbill et al. Reference Brackbill, Kothe and Zemach1992).

The reaction term, ${\dot{m}}_{l,A}$ , describes the chemical reaction of the mineral surface and the acid species included in the liquid phase. In (2.21), this term is negative because acid concentration decreases at the liquid/solid interface due to the chemical reaction. For a first-order reaction rate, we have

(2.24) $$\begin{eqnarray}\displaystyle {\dot{m}}_{l,A}=\unicode[STIX]{x1D70C}_{l}a_{v}^{l}r\bar{\unicode[STIX]{x1D714}}_{A}, & & \displaystyle\end{eqnarray}$$

where $r$ is the constant of reaction and $a_{v}^{l}$ is the wetted effective surface area in one grid block. The latter is a function of the specific area, $a_{v}$ and the liquid saturation. In the micro-continuum formulation for pore-scale simulation with a sharp interface between the non-obstructed regions and the low-permeability low-porosity regions, the surface area per control volume can be evaluated using the porosity gradient, $a_{v}=\Vert \unicode[STIX]{x1D735}\unicode[STIX]{x1D700}\Vert$ (Soulaine et al. Reference Soulaine, Roman, Kovscek and Tchelepi2017). To enforce the reactive boundary solution at the fluid/solid interface only, the specific area is computed as $a_{v}=\Vert \unicode[STIX]{x1D735}\unicode[STIX]{x1D700}\Vert \unicode[STIX]{x1D713}$ , where $\unicode[STIX]{x1D713}$ is a diffuse interface function. We use $\unicode[STIX]{x1D713}=4\unicode[STIX]{x1D700}(1-\unicode[STIX]{x1D700})$ (Luo et al. Reference Luo, Quintard, Debenest and Laouafa2012). Moreover, because the reaction occurs on wetted surface only, we simply write $a_{v}^{l}=a_{v}S_{l}$ .

2.6 Numerical implementation

The mathematical model formed by (2.13)–(2.15) and (2.21) has been implemented in the open-source computational fluid dynamics platform OpenFOAM^® (http://www.openfoam.org). The equations are first discretized using a finite-volume method and solved sequentially. The pressure–velocity coupling is handled by a predictor–corrector strategy based on the pressure-implicit with splitting of operators (PISO) algorithm (Issa Reference Issa1985). The full solution algorithm is described in appendix C.

3.1 Experiments: materials and method

A calcite post embedded in a single straight polydimethylsiloxane (PDMS) microchannel is used to visualize directly the dissolution dynamics. The calcite post is a regular octagonal prism with a width of 0.5 mm and a height of 0.25 mm. The microchannel is 37 mm long with a rectangular cross-section of $1.5~\text{mm}\times 0.2~\text{mm}$ . The well-known technique of soft lithography is used to make the microchannels (McDonald et al. Reference McDonald, Duffy, Anderson, Chiu, Wu, Schueller and Whitesides2000). Briefly, a glass-chrome photomask is used to create the positive relief of the microchannel etched in silicon (Buchgraber et al. Reference Buchgraber, Al-Dossary, Ross and Kovscek2012; Roman et al. Reference Roman, Soulaine, AlSaud, Kovscek and Tchelepi2016), that is then used as a master for the PDMS moulding process. A cast of the channel is obtained by pouring transparent liquid PDMS (Silicone Elastomer, Sylgard) onto the silicon master, baking it for 60 min at $65\,^{\circ }\text{C}$ to harden it, and peeling it from the master. Holes for the fluidics connections are made in the PDMS. Then, the PDMS is cleaned with soap and carefully dried. The calcite crystal is placed at the centre of the channel using a set of Teflon tweezers. Then, another cast of PDMS is prepared by pouring liquid PDMS into a flat mould and partially curing it (30–40 min, $55\,^{\circ }\text{C}$ ), (Eddings, Johnson & Gale Reference Eddings, Johnson and Gale2008; Roman et al. Reference Roman, Lorthois, Duru and Risso2012). The semi-cured PDMS had the ability to deform and create a seal around the calcite crystal. Eventually, the flat plate is placed on top of the calcite-embedded open channel and heated at $110\,^{\circ }\text{C}$ for 15 min until the two PDMS layers are securely bonded.

A polytetrafluoroethylene (PTFE) tube is connected on one side to the PDMS microchannel and on the other to a syringe (BD, 10 ml) filled with acid solution. The microchannel is initially saturated with air. A syringe pump (Harvard Apparatus Pump 11 Elite) provides a constant injection rate into the microchannel. The fluid injected is a mixture of deionized water, ethanol (EtOH) and hydrogen chloride (HCl). The PDMS channel walls are hydrophobic, thus ethanol is added to the acid mixture in order to wet the PDMS walls. Two sets of experimental conditions are used: (i) 0.5 % HCl concentration and $1.7\times 10^{-10}~\text{m}^{3}~\text{s}^{-1}$ flow rate ( $5.8\times 10^{-4}~\text{m}~\text{s}^{-1}$ mean velocity), (ii) 1 % HCl concentration and $6.9\times 10^{-10}~\text{m}^{3}~\text{s}^{-1}$ flow rate ( $2.3\times 10^{-3}~\text{m}~\text{s}^{-1}$ mean velocity).

The microfluidic device is placed under a Nikon ME600 microscope for direct visualization of the dissolution dynamics. The objective lens had a magnification of 4 and a numerical aperture of 0.13. Sequences of images are recorded using a pco.edge 5.5 camera (sCMOS sensor up to 50 f.p.s.).

Sequences of images of the dissolution process are analysed in order to measure the evolution of the size of the calcite post and to measure the evolution of $\text{CO}_{2}$ bubbles over time. The image processing workflow was developed using the image processing toolbox of MATLAB^®. The parameters measured for each time step are: size of the calcite post, number of $\text{CO}_{2}$ bubbles connected to the calcite and size of $\text{CO}_{2}$ bubbles connected to the calcite. The dissolution is assumed to be uniform in the depth of the channel (Soulaine et al. Reference Soulaine, Roman, Kovscek and Tchelepi2017) and only the top view of the dissolution process is visualized and measured for post-processing.

3.2 Experimental results

Figure 4.

Image sequence of the formation of $\text{CO}_{2}$ bubbles during the acidifying process of a calcite crystal. Experimental conditions: 1 % HCl concentration, $6.9\times 10^{-10}~\text{m}^{3}~\text{s}^{-1}$ flow rate.

Figure 4 shows images of the calcite dissolution with production of gas bubbles at different times during the injection at $6.9\times 10^{-10}~\text{m}^{3}~\text{s}^{-1}$ of a solution containing 1 % of HCl. Similar results are obtained for the injection of 0.5 % of HCl at a lower rate. First, the sequence focuses on the events that occur just after the passage of the acid front, i.e. from 0 to about 15 s after the acid front reached the calcite. First, many small $\text{CO}_{2}$ bubbles nucleate at the crystal surface (figure 4, $t=0.45$  s, $t=0.88$  s). The nucleation sites are arranged uniformly all around the calcite surface. Then, as the chemical reaction continues, the bubbles grow and merge with their closest neighbours. After about 5 s, the bubble size distribution changes from large numbers of small bubbles to only a few – three or four – large bubbles (figure 4, $t=5.00$  s, $t=9.55$  s, $t=13.65$  s). The spacing between the largest bubbles remains relatively uniform. These larger bubbles are fed by smaller bubbles produced at the nucleation sites that are still exposed to the acidified aqueous phase. Some of these large bubbles come off the calcite surface, and they are transported away by the flow. At some point, one, or several, large bubbles continue to grow while no small bubbles are visible (figure 4, $t=13.65$  s, $t=261$  s, $t=399$  s). This suggests that very small bubbles are generated at the surface of the mineral, but migrate to the large ones before they become visible. Eventually, the large bubbles coalesce, detach and the sequence is repeated. For these experimental conditions, the calcite post is fully dissolved after 13 min. A similar behaviour is observed for the injection of 0.5 % HCl at $1.7\times 10^{-10}~\text{m}^{3}~\text{s}^{-1}$ (figure 6): $\text{CO}_{2}$ small bubbles are generated for approximately 35 s, then 3 or 4 large bubbles are observed, that will eventually coalesce and detach. In this case the calcite post is fully dissolved after approximately 50 min.

Figure 5.

Distribution of bubble sizes over time for two experimental conditions.

Quantitative measurements of the bubble size distribution are presented in figure 5 for both experiments. For both cases, only small bubbles (blue stars and red squares) are generated in the beginning. Then, these bubbles grow and coalesce and larger bubbles appear, after 7.5 and 3.5 s, figure 5(a,b) respectively. For the greater acid concentration and greater flow rate, we observe large bubbles only after 10.5 s of injection (figure 5 b, green and purple markers) and after 40 s for the other case (figure 5 a). The same regimes of bubble generation, i.e. from a lot of very small bubbles to a couple of medium sized bubbles and to one or two very large bubbles, are observed.

Figure 6.

Mass of calcite normalized by the initial mass ( $t=0$  s) and mass of $\text{CO}_{2}$ normalized by the maximum mass of $\text{CO}_{2}$ that can be contained in the field of view, i.e. normalized mass of $\text{CO}_{2}=1$ when bubbles occupy the whole channel width. Experimental conditions: 0.5 % HCl concentration, $1.7\times 10^{-10}~\text{m}^{3}~\text{s}^{-1}$ flow rate.

The evolution of the mass of the solid and of the gas are plotted in figure 6 for 0.5 % HCl concentration and $1.74\times 10^{-10}~\text{m}^{3}~\text{s}^{-1}$ flow rate. The area of the calcite normalized by the initial area is measured and is equivalent to the normalized mass of calcite. The area occupied by $\text{CO}_{2}$ bubbles, that are connected to the calcite, is normalized by the area of the channel in the field of view, meaning that when the gas occupies the whole space, this value is equal to 1. In figure 6 we observe that the mass of calcite decreases at a rate of approximately $0.025\,\%~\text{s}^{-1}$ from 0 to 300 s. At this point a very large bubble is protecting the calcite crystal from being exposed to the acid solution, see picture at 277 s on figure 6. Notes that the lack of experimental points from approximately 70 to 220 s is due to memory limitation of the camera when recording sequences of high-resolution images.

We also observe that the crystal geometry does not evolve to the petal-like shape identified during the dissolution of solid minerals under single-phase flow conditions (Molins et al. Reference Molins, Trebotich, Miller and Steefel2017; Soulaine et al. Reference Soulaine, Roman, Kovscek and Tchelepi2017). Instead of an elongated profile in the main direction of the flow, the shape of the crystal when $\text{CO}_{2}$ is produced as a separate phase is difficult to characterize and depends on the complex history of the gas bubble formation, from the nucleation to the drying of the mineral surface which stops the chemical reaction at this location.

3.3 Simulation results

The multiphase DBS framework is used to model the dissolution of the first experimental dataset (0.5 % HCl at $1.7\times 10^{-10}~\text{m}^{3}~\text{s}^{-1}$ ). The geometry consists in a $1.5~\text{mm}\times 1.5~\text{mm}$ two-dimensional box that contains a cylindrical pillar in the middle. Two grids are considered: the domain is discretized by a $100\times 100$ Cartesian grid and by a finer $200\times 200$ grid. The cell porosity is set to $\unicode[STIX]{x1D700}=0.001$ for the cells that make up a cylinder of 0.5 mm in diameter centred in the middle of the domain. Elsewhere, $\unicode[STIX]{x1D700}=1$ . The top and bottom boundaries of the domain are set as impermeable fixed walls with no reaction. Initially, the domain is saturated with liquid water both in the clear-fluid region and in the porous region ( $S_{l}=1$ everywhere). At the pillar boundary, a constant contact angle $\unicode[STIX]{x1D703}=45^{\circ }$ is imposed and the constant of reaction is set to $r=5\times 10^{-3}~\text{m}~\text{s}^{-1}$ based on the value obtained by Soulaine et al. (Reference Soulaine, Roman, Kovscek and Tchelepi2017). A solution of 0.5 % of acid is injected in the aqueous phase from the left-hand side of the domain at a constant velocity $v_{0}=5.8\times 10^{-4}~\text{m}~\text{s}^{-1}$ for 300 s. The right-hand side is an outflow condition. The other properties are $\unicode[STIX]{x1D70C}_{s}=2165~\text{kg}~~\text{m}^{-3}$ , $\unicode[STIX]{x1D70C}_{l}=1000~\text{kg}~~\text{m}^{-3}$ , $\unicode[STIX]{x1D70C}_{g}=1~\text{kg}~\text{m}^{-3}$ , $\unicode[STIX]{x1D707}_{l}=10^{-3}~\text{Pa}~\text{s}$ , $\unicode[STIX]{x1D707}_{g}=10^{-3}~\text{Pa}~\text{s}$ , $\unicode[STIX]{x1D70E}=0.03~\text{kg}~\text{s}^{-2}$ , $D_{l,A}=D_{g,A}=5\times 10^{-9}~\text{m}^{2}~\text{s}^{-1}$ , $H_{A}=1$ and $k_{0}=10^{-15}~\text{m}^{2}$ . To investigate the effect of the gas phase, simulations are also run under single-phase flow conditions assuming that the $\text{CO}_{2}$ produced at the mineral surface is entirely soluble in the aqueous phase.

Figure 7.

Simulation results for the evolution of a calcite grain during the injection of acid under single-phase and two-phase conditions.

Snapshots of the dissolution process are presented in figure 7 for the two grids and the single-phase dissolution. For the single-phase case, the solid grain evolves from a cylindrical shape to an oval with only one axis of symmetry. This evolution is due to the coupling with hydrodynamics and a faster dissolution rate upstream than downstream (Molins et al. Reference Molins, Trebotich, Miller and Steefel2017; Soulaine et al. Reference Soulaine, Roman, Kovscek and Tchelepi2017). For the multiphase simulations, the volume fraction of gas around the grain increases during the injection of acid. When enough $\text{CO}_{2}$ is produced, the surface tension forms small gas bubbles. Note that the interface between gas and liquid is resolved in the images when $S_{l}<0.5$ . The number of bubbles increases with the resolution of the computational grid. Similar to the sequence of events observed in the experiments, the bubbles grow and coalesce until three or four main bubbles occupy the domain with a uniform distribution around the calcite. The larger bubbles are fed by smaller gas bubbles that are not always visible in the image but are described with saturation between and eventually detach.

The evolution of the mass of solid and gas in the domain are plotted in figure 6 and compared with the experimental data. The prediction for the evolution of the mass of solid is in very good agreement with the experiment. Because the presence of the gas as a separate phase creates a barrier to the acid, the mean dissolution rate under two-phase condition is less important than under single-phase conditions. The rate of production of gas is also well captured by the multiphase DBS framework for the first 70 s. The sudden drops in the evolution of the mass of gas after that show, however, that the dynamics of the detachment is faster in the simulations than in the experiments. The fact that the simulations are two-dimensional while the experiment is actually three-dimensional can explain the premature detachments in the simulation. Indeed, the three-dimensional effects due to the confinement of the fluids between the top and bottom plates can play an important role in the multiphase dynamics. The modelling of the wall adhesion condition can also explain this difference and more accurate models for the contact angle could improve the prediction.

The successful comparison between the experimental measurements and the simulation results builds up confidence in the multiphase DBS predictions. In particular, the immersed boundaries used for the wettability conditions and the chemical reactions at the mineral surfaces capture the dynamics of the mineral dissolution.