Dissolution Mass Transfer of Trapped Phase in Porous Media

Abstract

Dissolution mass transfer of trapped phase (TP) to flowing phase (FP) in porous media plays significant roles in hydrogeology, e.g., groundwater contamination by non-aqueous phase liquids, groundwater in-situ bioremediation, and geological carbon sequestration. In this chapter, this phenomenon is described. First, the physical and mathematical models are given. Afterwards, various conditions affecting this process, i.e., porous media characteristics, capillary trapping characteristics, flow bypassing, TP characteristics, and FP velocity, are discussed. These various conditions are described based on three parameters affecting the dissolution mass transfer: TP interfacial area (A), TP dissolution ratio (ξ), and mass transfer coefficient (k). Eventually, models to predict the mass transfer are formulated based on non-dimensional model. All of the data in this chapter are based on the experiments obtained by using micro-tomography and a series of image processing techniques from our latest works.

1. Introduction

1.1 Applications of dissolution mass transfer in porous media

Dissolution mass transfer of trapped phase (TP) to flowing phase (FP) in porous media plays significant roles in hydrogeology (Figure 1), e.g., groundwater contamination by non-aqueous phase liquid (NAPL) [1, 2, 3], groundwater in-situ bioremediation [4, 5], and geological carbon sequestration (GCS) [6].

Figure 1.

Dissolution mass transfer in porous media phenomena in hydrogeology (modified from Patmonoaji et al. [7] with permission).

In groundwater contamination by NAPL, NAPL could leak from industrial plant and contaminate the groundwater stream [1, 2, 3]. When the NAPL enters the groundwater, it displaces the existing groundwater, and then it was trapped due to capillary trapping. As it is trapped, mass transfer from the trapped NAPL to the groundwater flow occurs, transferring the dissolved NAPL to the groundwater. As a result, groundwater is contaminated, and the period of contamination depends on the mass transfer rate of trapped NAPL to the groundwater. Low mass transfer rate results in a prolong contamination, and vice versa.

In-situ bioremediation is one of the decontamination method of NAPL contamination [4, 5]. This method harnesses the ability of native bacteria in the groundwater to biodegrade the contaminants. To enhance the biodegradation process, O₂, CH₄, NH₃, CO₃ and N₂ gases are injected at the upstream of the contamination area to supply nutrients and stimulate the bacteria. After the injection, these gases are trapped due to capillary trapping and thus mass transfer occurs. The mass transfer rate is important because it affects the efficiency of the gas transport. Higher mass transfer rate means faster delivery but faster depletion of trapped gases. Therefore, it should be controlled to match the bacteria needs.

In GCS, the captured CO₂ from the point source, e.g., coal power plant, steel factory, and cement factory, is injected into deep saline aquifers for sequestration. After the injection, the CO₂ is trapped in the groundwater due to capillary trapping [6, 8]. Although the CO₂ can be sequestered this way, the reservoir pressure increases due to the injection. However, by the time, the trapped CO₂ will dissolve, reducing the reservoir pressure. As a result, the mass transfer rate is critical to approximate the pressure reduction rate for safer CO₂ sequestration process.

Given the scale of the phenomena is massive, the observation scale can be categorized into pore-scale, core-scale, and field-scale (Figure 1). Pore-scale is in the order of μm, and core-scale is in the order of cm. Field-scale, however, is in the order of m to km. Each of the scale is used to explain the scales above it. Pore-scale and core-scale are performed in a laboratory, whereas field-scale is performed in the field. Only pore-scale and core-scale will be discussed in this chapter.

Given the opaqueness of porous media system, direct observation of dissolution process is difficult to be performed. To evaluate the mass transfer process in detail, the interfacial area ($\mathit{A}$) of TP needs to be measured. With the information of interfacial area of TP, the mass transfer coefficient ($\mathit{k}$) can be calculated (described later in Eq. (5)). In addition, the distribution of the TP clusters is also need to be monitored to assess the dissolution mass transfer process.

Earlier studies mainly relied only on core-scale studies by the measurement of dissolved phase concentration from the effluent and amount of trapped phase through gravimetrical measurement without knowing the trapped phase characteristics inside [1, 2, 3, 4, 5]. However, because the trapped phase characteristics are required to measure the $\mathit{A}$ for the mass transfer, and its characteristics are strongly controlled by the porous media characteristics, the studies become phenomenological because $\mathit{k}$ cannot be calculated. In addition, the overall dissolution process cannot be observed, and thus flow bypassing cannot be monitored although it also affects the mass transfer process.

In this chapter, dissolution mass transfer in porous media is discussed based on experimental observation on the dissolution mass transfer of various trapped gas into flowing water in porous media by using micro-tomography and image processing [7, 9, 10, 11, 12]. By using this methods, $\mathit{A}$ and distribution of TP clusters can be monitored throughout dissolution process, and thus $\mathit{k}$ can be calculated. In addition, by monitoring dissolution process, flow bypassing was investigated, and thus its effect to $\mathit{k}$ was found. Moreover, unique dissolution process of TP with high dissolution ratio ($\mathit{\xi }$), which is the ratio of solubility limit (${\mathit{C}}_{\mathit{sol}\mathit{,}\mathit{i}}$) and density (${\mathit{\rho }}_{\mathit{i}}$) was observed. Its effect to $\mathit{k}$ was also found. Therefore, all of three parameters effecting dissolution mass transfer ($\mathit{A}$,$\mathit{\xi }$, and $\mathit{k}$), were elucidated.

The contents of this chapter are described as follow. In Section 2, the physical and mathematical model are given. In Section 3, the general description about the micro-tomography and image processing methods are described. In Section 4, various parameters affecting this process, i.e., porous media characteristics, capillary trapping characteristics, flow bypassing, TP properties, and FP velocity, are discussed. All of these condition will be investigated based on the point of view of three parameters affecting dissolution mass transfer: $\mathit{A}$,$\mathit{\xi }$, and $\mathit{k}$. Porous media together with capillary trapping characteristics affect the amount of$\mathit{A}$ for mass transfer process. Flow bypassing affects the uniformity of mass transfer process and eventually $\mathit{k}$. TP properties affect the $\mathit{\xi }$ through its ${\mathit{C}}_{\mathit{sol}\mathit{,}\mathit{i}}$ and ${\mathit{\rho }}_{\mathit{i}}$. Lastly, FP velocity affects the advection rate of the solute, influencing $\mathit{k}$. Eventually, the mass transfer model based on non-dimensional number is also formulated at the end of the chapter.

2. Physical and mathematical model

2.1 What is dissolution mass transfer?

At first, we can start from the basic by describing mass transfer. Mass transfer is a process of a movement of mass from one medium to another medium. The examples can be found easily in the surrounding of our daily lives. For example, the mass transfer of camphor in a bathroom into the surrounding air to improve the smells of the room (Figure 2a). Another example is the mass transfer of oxygen from the air into a lake (Figure 2b). Without the supply of oxygen through mass transfer, the fishes living in the lake cannot breathe the oxygen and stay alive.

Figure 2.

(a) Mass transfer of camphor to surrounding air, (b) mass transfer of oxygen from air to water in lake, and (c) illustration of mass transfer process.

Similar with other transport phenomena process, i.e., heat transfer and momentum transfer, mass transfer occurs because of non-equilibrium condition [13, 14, 15]. In case of mass transfer, equilibrium can be achieved once the system possesses equal chemical potential ($\mathit{\eta }$) of the considered component ($\mathit{i}$):

$${\mathit{\eta }}_{\mathit{i},\mathbf{1}}={\mathit{\eta }}_{\mathit{i},\mathbf{2}}$$E1

with subscripts 1 and 2 represent the phases. This equation can also be described with partition coefficient or Henry’s coefficient (Figure 2c) as follow:

$${\mathit{C}}_{\mathit{i},\mathbf{1}}=\frac{{\mathit{C}}_{\mathit{i},\mathbf{2}}}{{\mathit{H}}_{\mathit{i}}^{\mathit{cc}}}$$E2

with ${\mathit{C}}_{\mathit{i}}$ as concentration and ${\mathit{H}}_{\mathit{i}}^{\mathit{cc}}$ as dimensionless Henry coefficient of component $\mathit{i}$, which is constant with pressure but varies with temperature.

The concept of mass transfer nowadays is mainly based on the works of Thomas Graham (1805–1869) and Adolf Fick (1829–1901). Thomas Graham was the pioneer in the study of mass transfer through his works in the diffusion of H₂ gas in a tube [16]. By using analogies with Fourier’s work [17] in heat conduction, Adolf Fick formulated the based equation for mass transfer that is the Fick’s law of diffusion mass transfer [18] defined as follow:

$$\mathit{j}=-{\mathit{D}}_{\mathit{i}}\mathbf{\nabla }{\mathit{C}}_{\mathit{i}}$$E3

with $\mathit{j}$ as mass flux and ${\mathit{D}}_{\mathit{i}}$ as diffusion mass transfer coefficient of compound $\mathit{i}$ in the medium. Another form is based on Newton’s law of cooling as follow:

$$\mathit{j}=\mathit{k}\mathbf{\Delta }{\mathit{C}}_{\mathit{i}}$$E4

with $\mathit{k}$ as mass transfer coefficient. In this approach, $\mathit{k}$ and $\mathbf{\Delta }{\mathit{C}}_{\mathit{i}}$ is used instead of ${\mathit{D}}_{\mathit{i}}$ and $\mathbf{\nabla }{\mathit{C}}_{\mathit{i}}$. As a result, the mass transfer is judged based on concentration difference without knowing the concentration gradient. Because the concentration gradient in this system is difficult to obtain in the system, this form of mass transfer equation is more preferable in mass transfer from an interface into bulk solution, such as in the hydrogeology processes mentioned in Section 1.

In case of dissolution mass transfer, an additional process present. Dissolution here refers to the decrease in volume due to mass transfer. As a result, the Eq. (4) can be modified to volume perspective as follow:

$${\mathit{\rho }}_{\mathit{i}}\frac{\mathit{dV}}{\mathit{dt}}=\mathit{kA}\left({\mathit{C}}_{\mathit{sol},\mathit{i}}-{\mathit{C}}_{\mathit{i}}\right)$$E5

with ${\mathit{\rho }}_{\mathit{i}}$ as density of compound $\mathit{i}$, $\mathit{V}$ as volume, and $\mathit{t}$ as time. Therefore, high mass transfer rate also results in high volume decrease rate or dissolution rate. The example of camphor in Figure 2a provides an example of dissolution mass transfer, whereas oxygen mass transfer in a pool of water in Figure 2b provides an example of mass transfer without dissolution.

2.2 Physical model of mass transfer from an interface into a bulk solution

There are several mechanisms models explaining mass transfer from an interface into a bulk solution [14]. Due to the analogy between mass transfer and heat transfer, similar models with heat transfer are also used. However, because most of heat transfer application is in fluid–solid system, whereas mass transfer applications occur in fluid–fluid system, not all model can be used for both.

The simplest physical model is the stagnant film theory (Figure 3a) [19] that assumes the existence of a stagnant film near the interface of mass transfer process. This stagnant film is also referred as an unstirred layer because no fluid motion exists in this film. However, this film is hypothetical because fluid motions still occur beyond the no-slip condition at interface. The effects of velocity and other parameters are then represented by the thickness of this stagnant film. Higher mass transfer rate results in thinner δ, and vice versa.

Figure 3.

(a) Illustration of stagnant film model, (b) penetration theory model at low and high velocity, and (c) surface renewal theory model.

Second mechanisms model is penetration theory (Figure 3b) [20]. The key assumption is that a thick moving film is continuously generated in the bulk fluid, which the mass transfer occurs by diffusion across this film. The effect of velocity to the mass transfer is represented by the concentration distribution. As the solute diffuse into the bulk liquid, it will be flushed away by the flow. Higher velocity results in faster removal of the solute. As a result, higher velocity results in steeper concentration distribution, leading to higher mass transfer rate, and vice versa.

Third mechanisms model is surface renewal theory (Figure 3c) [21], which is the further development of penetration theory. The bulk liquid is depicted into two regions, a region near the interface that is renewed quickly and a region of well-mixed bulk fluid. The renewed region is renewed by the flow and similar with penetration theory. However, it is constantly exchanged with new elements from a second bulk region, making this theory closer to the realistic physical situation.

Among these models, however, no theory is completely satisfactory. Nevertheless, these models provide valuable insight to the mechanisms in mass transfer. For more detail, readers are advised to check the articles in the references [19, 20, 21].

For more empirical modeling of the mass transfer process, an empirical equation based on dimensionless number is usually used. The experimental data of mass transfer is represented by power-law relation of dimensionless number as follow:

$$\mathit{Sh}={\mathit{ϵ}}_{\mathbf{1}}{\mathit{Re}}^{{\mathit{ϵ}}_{\mathbf{2}}}$$E6

with ${\mathit{ϵ}}_{\mathbf{1}}$, ${\mathit{ϵ}}_{\mathbf{2}}$, and ${\mathit{ϵ}}_{\mathbf{3}}$ as constants, $\mathit{Sh}$ as Sherwood number representing the mass transfer rate and $\mathit{Re}$ as Reynolds number representing the bulk fluid velocity and geometry of the system. These non-dimensional number are defined as follow:

$$\mathit{Sh}=\frac{\mathit{kL}}{{\mathit{D}}_{\mathit{i}}}$$E7

$$\mathit{Re}=\frac{{\mathit{\rho }}_{\mathit{FP}}\mathit{UL}}{{\mathit{\mu }}_{\mathit{FP}}}$$E8

with ${\mathit{\rho }}_{\mathit{FP}}$ as FP density, ${\mathit{\mu }}_{\mathit{FP}}$ as FP absolute viscosity, $\mathit{U}$ as FP velocity, and $\mathit{L}$ as characteristic length.

In case of porous media, the physical mechanism of mass transfer is actually similar. In an open space, when dissolution occurs, the solute can be flushed away easily. In porous media, on the other hand, these processes are governed by the solute transport in porous media. In addition, in the application of hydrogeology, the TP clusters are distributed throughout the porous media (Figure 4), resulting in continuous mass transfer process along the porous media. As a result, a special attention to the solute transport throughout the porous media is required before assessing the mass transfer process with Eq. (5).

Figure 4.

Illustration of dissolution mass transfer in porous media and the predicted ${\mathit{C}}_{\mathit{i}\mathit{,}\mathit{x}\mathit{,}\mathit{t}}$ based on equation 16.

2.3 Fluid flow and solute transport in porous media

The fluid flow in hydrogeology application is relatively slow and can be classified as Stokes flow or creeping flow, which the ratio between fluid velocity and viscosities is relatively low. This condition is represented by $\mathit{Re}$ lower than one. Steady-state Stokes equation is given as follow:

$${\mathit{\mu }}_{\mathit{FP}}{\mathbf{\nabla }}^{\mathbf{2}}\mathit{U}=\mathbf{\nabla }\mathit{P}-{\mathit{\rho }}_{\mathit{FP}}\mathit{g}$$E9

with $\mathit{P}$ as pressure and $\mathit{g}$ as gravitational acceleration. This equation is then further derived into Darcy law, which is introduced by Henry Darcy [22] as a relationship to describe fluid flux ($\mathit{q}$) in sand filters, which was first proposed by Nutting [23] and Wyckoff et al. [24], described as follow:

$$\mathit{q}=\frac{\mathit{K}}{{\mathit{\mu }}_{\mathit{FP}}}\left(\mathbf{\nabla }\mathit{P}-{\mathit{\rho }}_{\mathit{FP}}\mathit{g}\right)$$E10

with $\mathit{K}$ as permeability. Given its original development, Darcy equation has been used mainly in the field of fluid flow in porous media. Similarly, with Stokes flow, to satisfy Darcy regime, the $\mathit{Re}$ needs to be low. Experimental tests demonstrated that it can still be achieved when the $\mathit{Re}$ is lower than ten [25].

For the solute transport in porous media, it is affected by the tortuosity and velocity profile in the pore networks. Nevertheless, under Darcy regime, it can be predicted by using advection-dispersion equation [26] as follow:

$$\frac{\mathit{\partial }{\mathit{C}}_{\mathit{i}}}{\mathit{\partial t}}={\mathit{D}}_{\mathit{hdl}}{\mathbf{\nabla }}^{\mathbf{2}}{\mathit{C}}_{\mathit{i}}-\mathrm{q}\mathbf{\nabla }{\mathit{C}}_{\mathit{i}}-\mathit{R}$$E11

with, ${\mathit{D}}_{\mathit{hdl}}$ as hydrodynamic dispersion coefficient and $\mathit{R}$ as external mass transfer source or sink. With this equation, solute transport along the porous media can be predicted, and it can be used to assess the mass transfer with Eq. (5).

2.4 Calculation of dissolution mass transfer in porous media

To evaluate the dissolution mass transfer in porous media, combination of modified version of Eqs. (5) and (11) are used [1, 2, 3, 4, 5, 9]. To calculate the mass transfer coefficient in porous media, Eq. (5) is slightly modified as follow:

$${\mathit{\rho }}_{\mathit{i}}\mathit{\varphi }\frac{\mathit{d}{\mathit{S}}_{\mathit{TP}}}{\mathit{dt}}=\mathit{k}{\mathit{a}}_{\mathit{TP}}\left({\mathit{C}}_{\mathit{sol},\mathit{i}}-{\mathit{C}}_{\mathit{i}}\right)$$E12

with $\mathit{\varphi }$ as porosity, ${\mathit{S}}_{\mathit{TP}}$ as TP saturation, and ${\mathit{a}}_{\mathit{TP}}$ as specific interfacial area, which is the amount of $\mathit{A}$ per porous media volume. However, the ${\mathit{C}}_{\mathit{i}}$, which is the driving force of the mass transfer, is unknown. To approximate it, ${\mathit{C}}_{\mathit{i}}$ can be calculated using Eq. (11). Because the experiment can be assumed as a one-dimensional column, one-dimensional advection-dispersion equation as follow:

$$\mathit{\varphi }\left(\mathbf{1}-{\mathit{S}}_{\mathrm{T}\mathit{P},\mathit{cs}}\right)\frac{\mathit{\partial }{\mathit{C}}_{\mathit{i}}}{\mathit{\partial t}}=\frac{\mathit{\partial }}{\mathit{\partial x}}\left[\mathit{\varphi }\left(\mathbf{1}-{\mathit{S}}_{\mathit{TP},\mathit{cs}}\right){\mathit{D}}_{\mathit{hdl}}\frac{\mathit{\partial }{\mathit{C}}_{\mathit{i}}}{\mathit{\partial x}}\right]-\mathit{q}\frac{\mathit{\partial }{\mathit{C}}_{\mathit{i}}}{\mathit{\partial x}}-{\mathit{\rho }}_{\mathit{i}}\mathit{\varphi }\frac{\mathit{\partial }{\mathit{S}}_{\mathit{TP},\mathit{cs}}}{\mathit{\partial t}}$$E13

with ${\mathit{S}}_{\mathit{TP}\mathit{,}\mathit{cs}}$ as cross-sectional average of TP saturation and $\mathit{x}$ as axial direction, can be used. In addition, it can be simplified with some assumptions [1, 2, 3, 9]. First, because dissolution mass transfer is a prolong process, the problem can be approached as a pseudo-steady state, neglecting the time accumulation term at left-hand side of the equation. Second, the dispersion term, which is the first term at the right-hand side of the equation, can be neglected due to the dominance in advection. It can be assessed by calculating the Peclet number ($\mathit{Pe}$) as follow:

$$\mathit{Pe}=\frac{\mathit{U}{\mathit{d}}_{\mathit{par},\mathbf{50}\%}}{{\mathit{D}}_{\mathit{i}}}$$E14

with ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$ as median particle size. $\mathit{Pe}$ represents the ratio of advection and diffusion of solute transport. High $\mathit{Pe}$ number represents advection dominant regime, whereas low $\mathit{Pe}$ represents diffusion dominant regime. When the value of $\mathit{Pe}$ is higher than five, it represent advection dominant regime [27]. Thus, the dispersion term can be neglected. This was also proven with order of magnitude analysis [3, 9]. The order of magnitudes of accumulation and dispersion terms are significantly lower than the others. Hence, it can be simplified to:

$$\mathit{q}\frac{\mathit{\partial }{\mathit{C}}_{\mathit{i}}}{\mathit{\partial x}}=-{\mathit{\rho }}_{\mathit{i}}\mathit{\varphi }\frac{\mathit{\partial }{\mathit{S}}_{\mathit{TP},\mathit{cs}}}{\mathit{\partial t}}$$E15

By discretizing $\mathit{\partial x}$ to $\mathit{\Delta x}$, ${\mathit{C}}_{\mathit{i}}$ at the location of $\mathit{x}\mathit{+}\mathit{\Delta x}$ and at the given time of $\mathit{t}$ (${\mathit{C}}_{\mathit{i}\mathit{,}\mathit{x}\mathit{+}\mathit{\Delta x}\mathit{,}\mathit{t}}$) as shown in Figure 4, can be calculated as follow:

$${\mathit{C}}_{\mathit{i},\mathit{x}+\mathit{\Delta x},\mathit{t}}={\mathit{C}}_{\mathit{i},\mathit{x},\mathit{t}}-\frac{{\mathit{\rho }}_{\mathit{i}}\mathit{\varphi }}{\mathit{q}}\frac{\mathit{\Delta }{\mathit{S}}_{\mathit{TP},\mathit{cs}}}{\mathit{\Delta t}}\mathit{\Delta x}$$E16

The average of concentration in between $\mathit{x}$ and $\mathit{x}\mathit{+}\mathit{\Delta x}$ at the time of $\mathit{t}$ (${\mathit{C}}_{\mathit{i}\mathit{,}{\mathit{x}}_{\mathit{avg}}\mathit{,}\mathit{t}}$) (Figure 4) then can be calculated by averaging the ${\mathit{C}}_{\mathit{i}\mathit{,}\mathit{x}\mathit{,}\mathit{t}}$ and ${\mathit{C}}_{\mathit{i}\mathit{,}\mathit{x}\mathit{+}\mathit{\Delta x}\mathit{,}\mathit{t}}$ as follow:

$${\mathit{C}}_{\mathit{i},{\mathit{x}}_{\mathit{avg}},\mathit{t}}=\frac{{\mathit{C}}_{\mathit{i},\mathit{x},\mathit{t}}+{\mathit{C}}_{\mathit{i},\mathit{x}+\mathit{\Delta x},\mathit{t}}}{\mathbf{2}}$$E17

This ${\mathit{C}}_{\mathit{i}\mathit{,}{\mathit{x}}_{\mathit{avg}}\mathit{,}\mathit{t}}$ then can be used to calculate the mass transfer coefficient as follow:

$$\mathit{k}=\frac{-{\mathit{\rho }}_{\mathit{i}}\mathit{\varphi }\frac{∆{\mathit{S}}_{\mathit{TP},\mathit{cs}}}{∆\mathit{t}}}{{\mathit{a}}_{\mathit{TP},\mathit{cs},{\mathit{t}}_{\mathit{avg}}}\left({\mathit{C}}_{\mathit{i},\mathit{sol}}-{\mathit{C}}_{\mathit{i},{\mathit{x}}_{\mathit{avg}},\mathit{t}}\right)}$$E18

with ${\mathit{a}}_{\mathit{TP}\mathit{,}\mathit{cs}\mathit{,}{\mathit{t}}_{\mathit{avg}}}$ as the cross-sectional average of ${\mathit{a}}_{\mathit{TP}}$ that is averaged in between the time of $\mathit{t}\mathit{-}\mathit{∆}\mathit{t}$ (${\mathit{a}}_{\mathit{TP}\mathit{,}\mathit{cs}\mathit{,}\mathit{t}\mathit{-}\mathit{∆}\mathit{t}}$) and $\mathit{t}$ (${\mathit{a}}_{\mathit{TP}\mathit{,}\mathit{cs}\mathit{,}\mathit{t}}$) calculated as follow:

$${\mathit{a}}_{\mathit{TP},\mathit{cs},{\mathit{t}}_{\mathit{avg}}}=\frac{{\mathit{a}}_{\mathit{TP},\mathit{cs},\mathit{t}}+{\mathit{a}}_{\mathit{TP},\mathit{cs},\mathit{t}-∆\mathit{t}}}{\mathbf{2}}$$E19

For the modeling of mass transfer, $\mathit{Re}$ and $\mathit{Sh}$ in Eqs. (7) and (8) were used with $\mathit{L}$ equal to ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$.

3. Experimental procedure

3.1 Experimental setup

As mentioned in Section 1, this chapter is based on the experimental observation by using micro-tomography and image processing [7, 9, 10, 11, 12]. For the micro-tomography, an X-ray microtomography scanner (micro-CT) was used. This micro-CT is able to generate 3D images in the form of 992 cross section slice image with the dimension of 992 × 992 voxels with the voxel size of 16.427 μm. As a result, observation area with the size of 16.34 × 16.34 × 16.34 mm³ was generated.

To generate the porous media, cylindrical container with an inner diameter of 10 mm was used as the casing of porous media. Inside this container, various particles will be packed inside to generate various porous media. All of the packing was performed under vibration to generate closed random packing classification [28]. As a result, the quality of the packing can be controlled. Additionally, sintered glass discs were placed at inlet and outlet to generate uniform flow.

To generate capillary trapping and similar condition with hydrogeology conditions in Section 1, the porous media will be sequentially saturated with degassed deionized (DDI) water and the gas. First, the gas was injected into the porous medium to saturate the porous media and remove all of other gases. Second, DDI water was injected upward to saturate the porous media, generating a condition similar with groundwater. Third, the experimental gas was injected again to remove most of the DDI water, resulting in similar condition with gas injection into groundwater. Finally, the DDI water was injected upward at the desired experimental flow rate. With this method, the generated trapped gas is uniform along the sample. Once the DDI water had breakthrough, and no gas flow was observed in the outlet tube, creating a quasi-steady state, the experiment began, and the dissolution process was monitored with micro-CT. The schematic about the experimental setup is given in Figure 5.

Figure 5.

Schematics of experimental setup with micro-tomography (modified from Patmonoaji et al. [9] with permission).

3.2 Image processing

For the data processing, image processing techniques were carried out to quantify the data from the micro-CT image. The main method is to differentiate the gas from the surroundings by performing image binarization (Figure 6). Although the contrast between water and particles is indistinguishable, the contrast of the gas among other phases is easily distinguishable. Therefore, a direct thresholding method was enough to separate gas from other phases. With this binarized images, then quantitative data, i.e., amount of trapped gas, and surface area of the trapped gas, can be simply derived with image processing to assess the dissolution mass transfer. Furthermore, by identifying separated trapped gas clusters, each of the trapped gas clusters can be characterized based on its volume and surface area. In addition, by performing similar method to porous media saturated with gas, porosity of the porous media can be obtained. Moreover, by applying watershed (WS) segmentation method [11], the throats and pores of the porous media can be identified and measured.

Figure 6.

Binarization image processing methods: (a) cross sectional image of the porous media, (b) histogram of the gray value, (c) binarized image with black as the gas, and (d) 3D binarized image of gas (modified from Patmonoaji et al. [1] with permission).

4. Factors affecting dissolution mass transfer in porous media

4.1 Porous media characteristics

A porous medium can be described as a solid containing void networks inside it. Structure-wise, the void network can be seen as a series of pores and throats (Figure 7a) [28, 29]. Pores are connected by a throat, which is a constriction in the network. Because capillary trapping mainly occupies pores, pore characteristics affect the morphology of the capillary trapping that eventually controls the amount $\mathit{A}$ for the mass transfer process (Figure 7b).

Figure 7.

(a) Illustration of pore and throat in pore network and (b) effect of porous media characteristics to TP morphology.

These pore and throat in porous media are mainly governed by the shape, size, and distribution of the granular particles forming the porous media. Shown in Table 1 are the properties of nine granular particles used to generate porous media [10, 11]. Those properties ${\mathit{C}}_{\mathit{u}}$ as particle distribution uniformity coefficient, ${\mathit{C}}_{\mathit{c}}$ as particle distribution gradation coefficient, and $\mathit{\Phi }$ as porosity. Three particle shapes, which are plastic beads (PB), glass beads (GB), and Toyoura sands (TS), are used. PB represents angular shape particle like quartz sands, whereas GB represents sphere particles. TS, on the other hand, represents natural silica sands.

Type	Size range (μm)	$d_{{par},50\%}$ (μm)	$C_u$	$C_c$	$\Phi$
PB	125–300	184	1.5	0.9	0.43
	125–710	212	1.8	0.8	0.41
	180–212	196	1.1	1.0	0.44
	355–425	390	1.1	1.0	0.45
	600–710	655	1.1	1.0	0.44
GB	180–212	196	1.1	1.0	0.33
	355–425	390	1.1	1.0	0.34
	600–710	655	1.1	1.0	0.35
TS	180–212	196	1.1	1.0	0.32

Table 1.

Properties of nine porous media generated from different shape, size, and distribution of granular particles.

Figure 8a shows the pore size distribution (PSD) of the porous media generated from those particles. The PSD is shown as probability density function (PDF) of pore equivalent diameter (${\mathit{d}}_{\mathit{\text{pore}}}$). To obtain this PSD, combinations of micro-tomography and WS method were performed [11]. GB tends to generate larger ${\mathit{d}}_{\mathit{\text{pore}}}$ than PB, whereas TS generate similar ${\mathit{d}}_{\mathit{\text{pore}}}$ with GB. For particle size, larger particles generate larger ${\mathit{d}}_{\mathit{\text{pore}}}$. However, in case of particle size distribution as compared by PB 180–212 μm, PB 125–300 μm, and PB 125–710 μm, similar PSD was found. Although the range of particle size distribution is different, the ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$ shown in Table 1, is actually similar. Pores are voids generated due to contact of granular particles (Figure 7b). Contact of large particles generates large voids. However, if the distribution of particles is wide. Large voids that is generated by large particles can be filled with the smaller particles. As a result, ${\mathit{d}}_{\mathit{\text{pore}}}$ tends to be controlled by the smaller particle size, particularly ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$. This is also consistent with the $\mathit{\Phi }$. Wider particle range generates lower $\mathit{\Phi }$ because the large pores generated from the contact of large particles are filled by smaller particles.

Figure 8.

(a) PSD and (b) ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$-normalized PSD.

Proportionality of the PSD is also found by normalizing PSD with the ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$ (Figure 8b). By normalizing the PSD with ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$, the PSD tend to converge into a single curve but not exactly. The differences are due to particle shape variation. Nevertheless, it demonstrates proportionality to ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$ for each particle shape.

The pore shape can also be evaluated through its surface area (${\mathit{A}}_{\mathit{\text{pore}}}$) and volume (${\mathit{V}}_{\mathit{\text{pore}}}$) as shown in Figure 9. Each point is the result obtained from averaging ${\mathit{A}}_{\mathit{\text{pore}}}$ of all pores in a log-scale range of ${\mathit{V}}_{\mathit{\text{pore}}}$. Typical pores shape detected in porous media with PB, GB, and TS is also given as yellow, blue, and green objects, respectively. As shown in normal-scale, pores in smaller ${\mathit{d}}_{\mathit{\text{pore}}\mathit{,}\mathbf{50}\%}$ generates larger ${\mathit{A}}_{\mathit{\text{pore}}}$. Related with particle shape, the ${\mathit{A}}_{\mathit{\text{pore}}}$ of PB tends to be higher than GB and TS. GB generates pore with the smallest ${\mathit{A}}_{\mathit{\text{pore}}}$, whereas TS generate pore with ${\mathit{A}}_{\mathit{\text{pore}}}$ in between PB and GB. This higher ${\mathit{A}}_{\mathit{\text{pore}}}$ is due to the angular shapes of PB.

Figure 9.

${\mathit{V}}_{\mathit{\text{pore}}}$ and ${\mathit{A}}_{\mathit{\text{pore}}}$ in the porous media and the typical pore shapes in each particle shape (modified from Patmonoaji et al. [11]).

4.2 Capillary trapping characteristics

Capillary, buoyancy, and viscous forces play important roles in the formation of capillary trapping (Figure 10a). These forces can be represented with capillary number ($\mathit{Ca}$) and Bond number ($\mathit{Bo})$ as follow:

Figure 10.

(a) Illustration of forces exerting on TP in porous media and morphology of (b) trapped NWP and (c) WP in porous media.

$$\mathit{Ca}=\frac{\mathit{U}{\mathit{\mu }}_{\mathit{WP}}}{\mathit{\gamma }}$$E20

$$\mathit{Bo}=\frac{∆\mathit{\rho g}{{\mathit{d}}_{\mathit{par},\mathbf{50}\%}}^2}{\mathit{\gamma }}$$E21

$\mathit{Ca}$ represents the comparison between viscous and capillary forces, whereas $\mathit{Bo}$ represents the comparison between buoyancy and capillary forces. Buoyancy and viscous forces exert on the TP to escape from the trapped location, whereas capillary force resist these forces to keep the TP in the pore. If $\mathit{Bo}$ and $\mathit{Ca}$ is higher than 0.35 and 10⁻³, respectively, no capillary trapping occurs due to excessive buoyancy and viscous forces [30, 31]. In hydrogeology applications, the $\mathit{Bo}$ and $\mathit{Ca}$ are much lower, resulting in the possibility of capillary trapping generation.

Wettability of the porous media also affects the TP clusters morphology (Figure 10b and c). In case the TP is non-wetting phase (NWP), the TP is located in the center of the pore. On the other hand, in case the TP is wetting phase (WP), the TP tends to surround the solid. In this work, we mainly focus on the former condition.

The trapped phase size distribution (TPSD) in the porous media described in Table 1 is given in Figure 11 as PDF of TP equivalent diameter (${\mathit{d}}_{\mathit{TP}}$). As shown in Figure 11a, with the increase in ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$, ${\mathit{d}}_{\mathit{TP}}$ increases. Compared with the PSD, the TPSD tends to be larger due to the possibility to form multi-pores bubbles. By normalizing with ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$ as shown in Figure 11b, the TPSD also converge demonstrating proportionality to ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$. This proportionality of TPSD with ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$ can be traced back to the proportionality of PSD to ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$ in Section 4.1. Because capillary trapping mainly occupies pore, the proportionality of PSD to ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$ also represents the proportionality of TPSD to ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$.

Figure 11.

(a) TPSD and (b) ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$-normalized TPSD (modified from Patmonoaji et al. [10]).

The TP cluster morphology can also be characterized through its volume (${\mathit{V}}_{\mathit{TP}\mathit{,}\mathit{\text{cltr}}}$) and surface area (${\mathit{A}}_{\mathit{TP}\mathit{,}\mathit{\text{cltr}}}$) as shown in Figure 12. Each point is the result obtained from averaging the ${\mathit{A}}_{\mathit{TP}\mathit{,}\mathit{\text{cltr}}}$ of all TP clusters in a log-scale range of ${\mathit{V}}_{\mathit{TP}\mathit{,}\mathit{\text{cltr}}}$. The examples of the TP clusters shape are also given in Figure 12. At smaller size, the markers from all porous media coincides one another. However, as the ${\mathit{V}}_{\mathit{TP}\mathit{,}\mathit{\text{cltr}}}$ increase, they branched out with TP clusters from smaller particle size exhibit larger ${\mathit{A}}_{\mathit{TP}\mathit{,}\mathit{\text{cltr}}}$. The point of the branch out is approximately at the size of ${\mathit{d}}_{\mathit{\text{pore}}\mathit{,}\mathbf{50}\%}$. When the TP cluster size is smaller than ${\mathit{d}}_{\mathit{\text{pore}}}$ of the porous media, it tends to occupy only one pore with a sphere-like shape. However, for larger cluster, it has to occupy more than one pore. As a result, it extends to the neighboring pores forming a channel of multi-pores TP cluster. These single-pore and multi-pores clusters can be modeled with a sphere and an elongated capsule model as follow:

Figure 12.

The average $A_{\mathit{TP},\mathit{\text{cltr}}}$ at the given $V_{\mathit{TP},\mathit{\text{cltr}}}$ and examples TP clusters (Patmonoaji et al. [10]).

$$A_{\mathit{TP},\mathit{\text{cltr}},\mathit{SP}}=4.836{V_{\mathit{TP},\mathit{\text{cltr}}}}^{\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}$$E22

$$A_{\mathit{TP},\mathit{\text{cltr}},\mathit{MP}}=1.047{d_{\mathit{\text{pore}},50\%}}^2+\frac{4}{d_{\mathit{\text{pore}},50\%}}{V}_{\mathrm{T}P,\mathit{\text{cltr}}}$$E23

For the ${\mathit{a}}_{\mathit{TP}}$, it is given in Figure 13a. The line represents the cross-sectional average ${\mathit{a}}_{\mathit{TP}}$ (${\mathit{a}}_{\mathit{TP}\mathit{,}\mathit{cs}}$), whereas the marker represents the overall average ${\mathit{a}}_{\mathit{TP}}$ (${\mathit{a}}_{\mathit{TP}\mathit{,}\mathit{ovr}}$). Figure 13a shows that smaller particles generate larger ${\mathit{a}}_{\mathit{TP}}$. As described previously, smaller particles generate smaller PSD and TPSD. The single-pore TP clusters are sphere with smaller diameter, and the multi-pores TP clusters are elongated capsule with smaller diameter. As a result, the TP clusters generate higher ${\mathit{a}}_{\mathit{TP}}$. For the effect of particle size distribution range, similarly with PSD and TPSD, it is mainly controlled by the ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$. For the particle shape, PB generates the highest ${\mathit{a}}_{\mathit{TP}}$, whereas GB and TS generate similar ${\mathit{a}}_{\mathit{TP}}$. Given the PSD and TPSD, PB generates smaller pore and TP clusters, and thus higher ${\mathit{a}}_{\mathit{TP}}$ was produced. GB and TS generate similar PSD but smaller than PB. Therefore, similar ${\mathit{a}}_{\mathit{TP}}$ was produced between GB and TS but lower than PB. In addition, as shown in Figure 9, pores of PB is more irregular, and thus it generates TP clusters with more surface area. On the other hand, the pores of GB and TS are much less irregular than PB, and look relatively similar. Therefore, based on these two reasons, PB generate higher ${\mathit{a}}_{\mathit{TP}}$, whereas GB and TS generate similar ${\mathit{a}}_{\mathit{TP}}$ but lower than PB. In addition, by normalizing ${\mathit{a}}_{\mathit{TP}}$ with ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$, the ${\mathit{a}}_{\mathit{TP}}$ converges into two lines (Figure 13b). All ${\mathit{a}}_{\mathit{TP}}$ of PB converge into a single line, whereas TS and GB tend to converge into another line. Hence, two models were generated:

Figure 13.

(a) ${\mathit{a}}_{\mathit{TP}}$ in all porous media and (b) ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$ normalized-${\mathit{a}}_{\mathit{TP}}$ with lines and markers represent ${\mathit{a}}_{\mathit{TP}\mathit{,}\mathit{cs}}$ and ${\mathit{a}}_{\mathit{TP}\mathit{,}\mathit{ovr}}$, respectively (taken with permission from Patmonoaji et al. [10]).

$$a_{\mathit{TP}}=2.884\frac{S_{\mathit{TP}}}{d_{\mathit{par},50\%}}$$E24

$$a_{\mathit{TP}}=1.889\frac{S_{\mathit{TP}}}{d_{\mathit{par},50\%}}$$E25

For the linearity between ${\mathit{a}}_{\mathit{TP}}$ and ${\mathit{S}}_{\mathit{TP}}$, it can be explained with these relations:

$$\frac{a_{\mathit{TP}}}{S_{\mathit{TP}}}=\Phi \frac{A_{\mathit{TP}}}{V_{\mathit{TP}}}$$E26

$$A_{\mathit{TP}}=\sum A_{\mathit{TP},\mathit{\text{cltr}},\mathit{MP}}$$E27

$$A_{\mathit{TP}}=\sum A_{\mathit{TP},\mathit{\text{cltr}},\mathit{SP}}=n_{\mathit{TP},\mathit{\text{cltr}}}{A}_{\mathit{TP},\mathit{\text{cltr}},\mathit{SP}}$$E28

$$n_{\mathit{TP},\mathit{\text{cltr}}}=\frac{V_{\mathit{TP}}}{V_{\mathit{TP},\mathit{\text{cltr}},\mathit{SP}}}$$E29

with ${\mathit{A}}_{\mathit{TP}}$ as trapped phase interfacial area, ${\mathit{V}}_{\mathit{TP}}$ as trapped phase volume, and ${\mathit{n}}_{\mathit{TP}\mathit{,}\mathit{\text{cltr}}}$ as the number of TP clusters. When the trapped gas is mainly multi-pores TP clusters, which the function of ${\mathit{A}}_{\mathit{TP}\mathit{,}\mathit{\text{cltr}}}$ is linear with ${\mathit{V}}_{\mathit{TP}\mathit{,}\mathit{\text{cltr}}}$ (Eq. (23)), linear function between ${\mathit{a}}_{\mathit{TP}}$ and ${\mathit{S}}_{\mathit{TP}}$ will be generated. In contrast, if it consists of large number of single TP cluster, the relation will be controlled by ${\mathit{n}}_{\mathit{TP}\mathit{,}\mathit{\text{cltr}}}$. Because ${\mathit{n}}_{\mathit{TP}\mathit{,}\mathit{\text{cltr}}}$ is a function of ${\mathit{V}}_{\mathit{TP}}$, the relation of ${\mathit{a}}_{\mathit{TP}}$ and ${\mathit{S}}_{\mathit{TP}}$ will be linear also.

As described in this section, the ${\mathit{a}}_{\mathit{TP}}$ generated in porous media can be predicted by using Eqs. (24) and (25), which is based on ${\mathit{S}}_{\mathit{TP}}$ and ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$. Given these two parameters only, the ${\mathit{a}}_{\mathit{TP}}$ can be predicted. Higher ${\mathit{S}}_{\mathit{TP}}$ generates higher ${\mathit{a}}_{\mathit{TP}}$. Smaller ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$ also generates higher ${\mathit{a}}_{\mathit{TP}}$. Porous media generated from more angular particles will also generates higher ${\mathit{a}}_{\mathit{TP}}$. As described in Eq. (12), higher ${\mathit{a}}_{\mathit{TP}}$ results in faster dissolution mass transfer. Therefore, porous media generated from angular and smaller particles leads to faster dissolution, especially at the condition of high ${\mathit{S}}_{\mathit{TP}}$.

4.3 Flow bypassing

In general, flow bypassing can be divided into two categories. First is bypassing due to heterogeneity of porous media [32, 33], and second is bypassing induced by dissolution fingering in homogeneous porous media [34, 35, 36]. Both of them direct water flow to preferential flow path, resulting in non-uniform dissolution mass transfer. As a result, it will affect the $\mathit{k}$. In this section, bypassing induced by dissolution fingering is examined through dissolution experiment performed in porous media from PB with various sizes and distributions [7] shown in Table 2. For the TP, N₂ gas was used.

Size range (μm)	125–150	180–212	355–425	600–710	125–300	250–425	125–710
$k$ category	low	low	high	high	low	high	high
$\lambda$ (mm)	4.806	6.308	12.935	21.68	5.745	12.001	7.582

Table 2.

Characteristics of particle, porous medium, flow, and mass transfer inside each porous medium sample.

In general, the mechanism of dissolution fingering development can be explained as follow. Initially, capillary trapping is relatively homogeneous, and thus the permeability tends to be homogeneous as well. However, a slight difference in permeability remains present. This slight difference generates preferential water flow, directing the water to flow to region with higher permeability. At this higher permeability region, water flow slightly faster, and thus the TP clusters dissolve slightly faster as well. As a consequence, it gradually increases the permeability of this slightly permeable region, enhancing the permeability differences gradually. This results in a positive feedback, and thus the region of this preferential path grows with the dissolution time, resulting in the development of dissolution finger.

In short, the possibility of dissolution fingering development can be evaluated by measuring the approximated width of dissolution finger ($\mathit{\lambda }$) and compare it with the given size of porous medium diameter sample (${\mathit{L}}_{\mathit{y}}$), which in this work is 10 mm. By using similar analysis for mineral dissolution instabilities based on linear stability analysis [37, 38], we can approximate the size of the dissolution finger by measuring the dimensionless growth rate ($\overline{\mathit{\omega }}$) using the calculation from Zhao et al. [36] as follow:

$$\begin{array}{c}\overline{\omega }=\frac{{\overline{p}}_{\mathit{axf}}^{\prime }}{\left(1+\beta \right)S_{N_2}}\left\{\frac{\left(1+\beta \right){\overline{p}}_{\mathit{axf}}^{\prime }}{1-{\overline{\alpha }}_L{\overline{p}}_{\mathit{axf}}^{\prime }}+\left[\left(1+\beta \right)+\left(\frac{1}{1-{\overline{\alpha }}_L{\overline{p}}_{\mathit{axf}}^{\prime }}-{\overline{\alpha }}_L\overline{m}\right)\left(1-\beta \right)\right]\sigma \right\}\\ -\frac{{\overline{p}}_{\mathit{axf}}^{\prime }}{\left(1+\beta \right)S_{N_2}}\left[\left(\overline{m}-\frac{{\overline{p}}_{\mathit{axf}}^{\prime }}{1-{\overline{\alpha }}_L{\overline{p}}_{\mathit{axf}}^{\prime }}-{\overline{\alpha }}_L{\overline{m}}^2\left(1-{\overline{\alpha }}_L{\overline{p}}_{\mathit{axf}}^{\prime }\right)\right)\frac{\left(1-\beta \right)}{1-{\overline{\alpha }}_L{\overline{p}}_{\mathit{axf}}^{\prime }}\right]\end{array}$$E30

$$\sigma =\frac{\sqrt{{\left(\frac{{\overline{p}}_{\mathit{axf}}^{\prime }}{1-{\overline{\alpha }}_L{\overline{p}}_{\mathit{axf}}^{\prime }}\right)}^2+\left(\frac{4{\overline{m}}^2\left(1-{\overline{\alpha }}_T{\overline{p}}_{\mathit{axf}}^{\prime }\right)}{1-{\overline{\alpha }}_L{\overline{p}}_{\mathit{axf}}^{\prime }}\right)}-\frac{{\overline{p}}_{\mathit{axf}}^{\prime }}{1-{\overline{\alpha }}_L{\overline{p}}_{\mathit{axf}}^{\prime }}}{2}$$E31

$$\beta ={\left(\frac{1-S_{N_2,0}-S_{\mathit{wi}}}{1-S_{\mathit{wi}}}\right)}^3$$E32

$${\overline{p}}_{\mathit{axf}}^{\prime }=\frac{U}{\sqrt{\mathit{\varphi \tau }{D}_{N_2}}}\sqrt{\frac{1}{{\beta }_0}}$$E33

$${\beta }_0=\frac{29.224\mathit{\varphi Re}{D}_{N_2}}{{d_{50}}^2}$$E34

$${\overline{\alpha }}_L=\frac{{\alpha }_L}{L^{\ast }}$$E35

$${\overline{\alpha }}_T=\frac{{\alpha }_T}{L^{\ast }}$$E36

$$L^{\ast }=\sqrt{\tau D_{N_2}{t}^{\ast }}$$E37

$$t^{\ast }=\frac{\varphi }{{\beta }_0}$$E38

with $\overline{\mathit{m}}$ as dimensionless wave number, ${\mathit{S}}_{\mathit{wi}}$ as irreducible water saturation, ${\mathit{S}}_{{\mathit{N}}_{\mathbf{2}}\mathit{,}\mathbf{0}}$ as initial gas saturation, $\mathit{\tau }$ as tortuosity, ${\mathit{\alpha }}_{\mathit{L}}$ as longitudinal dispersivity, ${\mathit{\alpha }}_{\mathit{T}}$ as transversal dispersivity. By calculating the $\overline{\mathit{\omega }}$ at the given value of $\overline{\mathit{m}}$, the graph of $\overline{\mathit{\omega }}$ at the given $\overline{\mathit{m}}$ can be generated as shown in Figure 14. At the value of $\overline{\mathit{\omega }}$ lower than zero, the finger will stabilize, whereas at the value $\overline{\mathit{\omega }}$ larger than zero, the finger will develop. The peak value of $\overline{\mathit{\omega }}$ represents the value of critical dimensionless wave number (${\overline{\mathit{m}}}_{\mathit{\text{critical}}}$), which the wave will grow fastest and most likely to be observed. Therefore, by finding the ${\overline{\mathit{m}}}_{\mathit{\text{critical}}}$ in each experiment cases, $\mathit{\lambda }$ in each experiment cases can be calculated as follow:

Figure 14.

The function of $\overline{\mathit{\omega }}$ and $\overline{\mathit{m}}$ for the prediction of dissolution finger width calculated for PB 600–710 μm.

$$\mathbf{\lambda }=\frac{\mathbf{2}\mathit{\pi }{\mathit{L}}^{\ast }}{{\overline{\mathit{m}}}_{\mathit{c}\mathrm{r}\mathit{\text{itical}}}}$$E39

The dissolution process inside PB 600–710 μm and PB 180–212 μm are given in Figure 15a as 3D images. In PB 600–710 μm, dissolution process occurs uniformly, whereas in PB 180–212 μm, dissolution occurs non-uniformly. The $\mathit{k}$ throughout the dissolution process is given in Figure 15b with the function to${\mathit{S}}_{\mathit{TP}\mathit{,}\mathit{cs}}$. Two groups of $\mathit{k}$ are found. First group is higher $\mathit{k}$, that occurs in PB 125–710 μm, PB 250–425 μm, PB 355–425 μm, and PB 600–710 μm. Second group is lower $\mathit{k}$ that occur in PB 125–300 μm, PB 125–150 μm, and PB 180–212 μm. From the predicted $\mathit{\lambda }$, the first group belong to the value of $\mathit{\lambda }$ larger than ${\mathit{L}}_{\mathit{y}}$, whereas the second group belong to the value of $\mathit{\lambda }$ smaller than ${\mathit{L}}_{\mathit{y}}$. When the size of $\mathit{\lambda }$ is smaller than ${\mathit{L}}_{\mathit{y}}$, dissolution fingers could fully develop, resulting in bypassing and lower $\mathit{k}$. This trend occurs except in PB 125–710 μm that was influenced by the particle size distribution. The classification of $\mathit{k}$ and predicted $\mathbf{\lambda }$ is given in Table 2.

Figure 15.

(a) Dissolution progress of (top) particle 600-710 μm and (bottom) particle 180-212 μm and (b) $\mathit{k}$ at the given saturation in all porous media (modified from Patmonoaji et al. [7]).

As explained in this section, in the presence of bypassing induced by dissolution fingering, $\mathit{k}$ is lower due to non-uniform dissolution process. Because this is also occur in homogeneous porous media, and the limiting condition is just the ratio of the container to the predicted size of dissolution fingering, it is expected to occurs in field-scale of hydrogeology processes.

4.4 Trapped phase properties

Another parameter that affect the dissolution mass transfer rate is $\mathit{\xi }$, which is the ratio of ${\mathit{C}}_{\mathit{i}\mathit{,}\mathit{sol}}$ and ${\mathit{\rho }}_{\mathit{i}}$. ${\mathit{\rho }}_{\mathit{i}}$ affects the decrease rate in volume of dissolution process, and ${\mathit{C}}_{\mathit{i}\mathit{,}\mathit{sol}}$ serves as the driving force of the mass transfer process. Higher $\mathit{\xi }$ results in faster dissolution mass transfer process. As comparison, $\mathit{\xi }$ varies for about 10⁻⁹–10⁻³ in mineral dissolution, 10⁻⁴–10⁻³ in NAPL dissolution, and 10⁻³–10⁻¹ in gas dissolution [36]. Timewise, this implicates that the dissolution of minerals could take months and years, whereas dissolution of NAPL could take days to weeks. In contrast, dissolution of gases could only take minutes to hours.

Dissolution in porous media is not only mass transfer process, but also solute transport. When mass transfer occurs from TP, the dissolve solute need to be transported away from the surrounding TP by the FP. Therefore, high value in $\mathit{\xi }$ could generate problems in the balance of mass transfer and solute transport.

To compare the effect $\mathit{\xi }$, CO₂ and N₂ gases were used. The ${\mathit{\rho }}_{\mathit{i}}$ of N₂ and CO₂ gases are 1.842 and 1.165 kg/m³, whereas the ${\mathit{C}}_{\mathit{i}\mathit{,}\mathit{sol}}$ are 1.449 and 0.0175 kg/m³, resulting in $\mathit{\xi }$ of 0.787 and 0.015, respectively. The differences in $\mathit{\xi }$ is more than fifty times. The experiments were performed with the same porous media and FP velocity, resulting in the same $\mathit{Re}$.

The axial average TP saturation throughout the dissolution process is given in Figure 16. The dissolution of N₂, and CO₂ gases behave differently. In case of N₂ gas, as the time proceeds, ${\mathit{S}}_{\mathit{TP}\mathit{,}\mathit{cs}}$ decrease from the initial value to zero. For CO₂, on the other hand, the ${\mathit{S}}_{\mathit{TP}\mathit{,}\mathit{cs}}$ does not dissolve completely but stop at the value about 0.03 before it dissolves completely, resulting in two dissolution stages. First dissolution stage occurs as the dissolution of ${\mathit{S}}_{\mathit{TP}\mathit{,}\mathit{cs}}$ to about 0.03. Second dissolution stage occurs as the complete dissolution of the remaining ${\mathit{S}}_{\mathit{TP}\mathit{,}\mathit{cs}}$. By varying $\mathit{Re}$, the same dissolution pattern of CO₂ is still observed [9].

Figure 16.

${\mathit{S}}_{\mathit{TP}\mathit{,}\mathit{cs}}$ of N₂ and CO₂ gas during dissolution process under Re of 0.04.

By observing the dissolution process from the 3D images of segmented TP clusters as shown in Figure 17, the different dissolution behavior can be observed in more detail. In case of N₂, TP clusters dissolve completely from multi-pores size. However, in case of CO₂, at first, the multi-pores TP clusters dissolve to single-pore TP clusters. Afterwards, second stage occurs, dissolving the remaining clusters.

Figure 17.

TP cluster of N2 and CO₂ observed at different dissolution time. Different colors indicate different clusters (modified from Patmonoaji et al. [9]).

For the $\mathit{k}$, additionally, there is a significant different between CO₂ and N₂ gas as shown in Figure 18. Given the same $\mathit{Re}$, the $\mathit{k}$ of CO₂ is more than one order of magnitude lower than N₂ gas. In addition, $\mathit{k}$ shows a dipping point in around 0.03, which is consistent with the two stages dissolution of CO₂. When the ${\mathit{S}}_{\mathit{TP}\mathit{,}\mathit{cs}}$ reaches 0.03, the dissolution process slowed down, and thus $\mathit{k}$ decreases. Then second dissolution stages occurs, dissolving the rest ${\mathit{S}}_{\mathit{TP}\mathit{,}\mathit{cs}}$, resulting in the increase of $\mathit{k}$.

Figure 18.

Average $\mathit{k}$ of all experiments with different gases and $\mathit{Re}$.

Compared with other gases, CO₂ exhibit about fifty times larger. When TP cluster dissolves, the solute should be transported away by fresh incoming FP. If the solute is not transported, the solute concentration surrounding the TP cluster accumulates, slowing down the mass transfer. During the dissolution of CO₂, this rapid dissolution is believed to generate area with high CO₂ solute concentration due to much faster mass transfer rate than advection rate. In addition, as the trapped gas dissolves rapidly, it quickly retracts to the center of pores, increasing the difficulty of the solute to be removed by FP. As the dissolution occurs rapidly, the incoming FP cannot compete with the rapid dissolution of TP. As a result, two dissolution stages occur. First stage occurs before the accumulation of solute, whereas the second stage occurs after the accumulated solute is advected away from the surrounding of TP.

As discussed in this section, although higher $\mathit{\xi }$ results in higher dissolution mass transfer rate, when it becomes too high, it could result in lower $\mathit{k}$. This condition can be represented by two dissolution stages that occurs due to competition between mass transfer rate and advection rate.

4.5 Flowing phase velocity

To investigate the effect of FP velocity, the experiments were also performed under different FP velocity ($\mathit{Re}$ between 0.0016 and 0.16) with N₂ and CO₂ gases. They were performed with porous media from PB 250–425 μm to avoid bypassing and focus the discussion to the effect of FP velocity. By using the $\mathit{Sh}$ and $\mathit{Re}$, the mass transfer can be modeled with Eq. (6) as follow (Figure 19):

Figure 19.

Comparison of predicted-$\mathit{Sh}$ and measured-$\mathit{Sh}$ from the model.

$$\mathit{Sh}=\mathbf{12.142}\mathit{Re}$$E40

$$\mathit{Sh}=\mathbf{0.386}{\mathit{Re}}^{\mathbf{0.685}}$$E41

The Eq. (40) is for N₂ gas, whereas Eq. (41) is for CO₂ gas. The difference is due to the effect of high $\mathit{\xi }$ in CO₂, resulting in two dissolution stages. Because the $\mathit{k}$ is also lower, the non-dimensional number model result in lower $\mathit{Sh}$ than N₂ gas.

In this section, the effect of FP velocity was shown. With the increase in FP velocity, $\mathit{k}$ increase. This relation can also be represented in non-dimensional number, as represented by $\mathit{Re}$ and $\mathit{Sh}$ in Eqs. (40) and (41). The effect of high $\mathit{\xi }$ is also observable as lower $\mathit{Sh}$ due to lower $\mathit{k}$.

5. Conclusion

Throughout this work, the physical mechanism and mathematical model of dissolution mass transfer of trapped phase in porous media was described. Various condition affecting it were also discussed from experiment data obtained from micro-tomography and a series of image processing techniques. The effects of these conditions were represented in three parameters: $\mathit{A}$,$\mathit{\xi }$, and $\mathit{k}\mathit{.}$

In case of porous media characteristics, the characteristics of porous media generated from various particle shape, size, and distribution were investigated. Smaller ${\mathit{d}}_{\mathit{par}\mathit{,}\mathbf{50}\%}$ results in higher ${\mathit{a}}_{\mathit{TP}}$. Porous media generated from more angular particles will also generates higher ${\mathit{a}}_{\mathit{TP}}$. In addition, ${\mathit{S}}_{\mathit{TP}}$ linearly effect the amount of ${\mathit{a}}_{\mathit{TP}}$. Therefore, porous media generated from angular and smaller particles leads to faster dissolution, especially at the condition of high ${\mathit{S}}_{\mathit{TP}}$.

Bypassing induced by dissolution fingering was also described. When the predicted $\mathit{\lambda }$ is smaller than ${\mathit{L}}_{\mathit{y}}$, dissolution finger can fully develop, resulting in bypassing and decrease in $\mathit{k}$. Because this occurs in homogeneous porous media, it is expected to occurs in field-scale of hydrogeology processes.

In addition, the effect of $\mathit{\xi }$ was also investigated. If the $\mathit{\xi }$ of the TP is high, different dissolution process, which is represented by two stage dissolutions, occurs, leading to decrease in $\mathit{k}$. The competition between mass transfer and advection rate is believed to be the reason of this different dissolution process.

Eventually, the effect of FP velocity based on $\mathit{Sh}$ and $\mathit{Re}$ was formulated. With the increase of FP velocity, $\mathit{k}$ increase as represented by the increase in $\mathit{Sh}$ at the given $\mathit{Re}$. In case of CO₂, the $\mathit{Sh}$ also increase with $\mathit{Re}$, but it is still lower than N_2. This is due to the effect two stage dissolution observed in CO₂, resulting in lower $\mathit{k}$.

Acknowledgments

This work was supported by JSPS KAKENHI with grant numbers 17H00790 and 20 J14975.