2.2 Digital rock model
We used a 3D digital rock model obtained from microtomographic images of
Berea sandstone (Fig. 1a). Berea sandstone is chosen due to its
relatively well-known properties (Øren & Bakke, 2003). The rock has a
relatively large mean grain size of ∼250 μm and consists of quartz,
feldspar, carbonates, and clay minerals. Its average pore size is ∼20 μm
with a relatively high size deviation (Jiang & Tsuji, 2015). The
simulation requires a rock model that is larger than its representative
elementary volume (∼1 mm for Berea sandstone), and the multi-pore nature
of the rock must be considered. Our rock model, obtained from the
database of Dong and Blunt (2009), is a stack of 320 tomographic images,
each one measuring 320 × 320 pixels. The resolution of the image is
5.345 µm; thus, the digital rock represents an actual size of 1.710 ×
1.710 × 1.710 mm.
2.3 Relative permeability calculation
The viscosity ratio M is defined by dividing the nonwetting
fluid’s viscosity (μ nw) by the wetting fluid’s
viscosity (μ w) (equation 1). In the LBM, μcan be calculated as follows:
\begin{equation}
\mu\ =\frac{1}{3}(\tau\ -\ \frac{1}{2})\text{\ \ \ \ \ \ \ \ \ \ }(5)\nonumber \\
\end{equation}where 𝜏 represents the relaxation parameter. Thus, the viscosities of
both fluids can be modified by changing the 𝜏 value. Because 𝜏 must be
larger than 0.5 for a positive viscosity and the simulation becomes
unstable as 𝜏 approaches 0.5, it is important that the simulation uses a
𝜏 value that ensures its stability and accuracy. In this study, we
altered M only by changing the 𝜏 value of the nonwetting fluid to
manipulate μ nw and keepingμ w constant throughout the simulation to preserve
its stability.
The capillary number Ca can be modified based on equation (2) by
changing the σ and μ nw parameters. Becauseμ nw is also used to determine M , we
achieved the desired Ca by adjusting σ alone. The average
velocity of the nonwetting fluid phase after the simulation has
converged to equilibrium, defined as a change less than log Ca =
±0.1, was used to calculate the Ca number. When the fluid
velocity change in the last 1,000 simulation steps was less than 2% for
all cases, the simulations were assumed to have converged at that time.
The time-averaged Ca value in the last 3,000 steps was calculated
for further investigation.
To remove the complications arising from wettability effects, we assumed
that the solid was completely nonwetting (θ = 180°) in all conditions.
Given that condition, the wetting fluid fills small pore spaces and
completely coats the surface of the rock, whereas the nonwetting fluid
is not in contact with the pore wall and can only occupy the central
parts of large pores.
We applied a steady-state simulation condition, which assumes that the
wetting fluid saturation (S w) and nonwetting
fluid saturation (S nw) are kept constant. In the
initial state, both fluids were specified as being randomly distributed
in the pore spaces. This initial condition was chosen because the
initial fluid connectivity is reduced, which can lead to increased
capillary trapping (Jiang & Tsuji, 2016). This initial condition is
desirable in a CCS site to increase the residual CO2saturation. Periodic boundary conditions were applied in the x ,y , and z directions of our 3D model. A constant body force
was applied in z direction to both fluids to mimic the pressure
gradient: g = dP /dz . For the interaction between
rock solid nodes and fluid voxels, no-slip boundary conditions were
applied using the halfway bounce back scheme. The density of both fluids
was set at 1.0 in lattice units, which corresponds to 1,000
kg/m3 in the physical unit. We did this because our
study was focused on viscosity differences, and the effect of density
contrast is minor if inertial force can be neglected.
We applied a periodic boundary condition in the flow direction;
therefore, the rock was mirrored in the z direction to ensure
that the pore spaces on the right side were connected to the left side
of the digital rock (Fig. 1b).
2.3.1 Relative permeability curve
The relative permeability curve is a plot of k versus saturationS . As S nw increases,k nw increases and k wdecreases. The shape of this curve can change as a result of viscous
drag force and capillary force; thus, an investigation of its variation
with changes in M and Ca is important to determine the
optimum conditions to achieve the desired value of k . We plottedk nw and k w as a function
of S nw to create relative permeability curves for
various conditions.
In this study, the wetting fluid viscosity was kept constant at 0.15 in
lattice units in all simulations, and M was manipulated by
altering only μ nw. Changes in M can be
interpreted as a change in either μ w orμ nw; thus, modifying either viscosity value will
produce the same results if their ratio is maintained.
It must also be noted that, because we conducted a steady-state
simulation, the S nw value used for plotting the
permeability curve might be unrealistic in actual injection settings.
Previous research has shown that the maximum S nwor irreducible S w depends on M andCa (Tsuji et al., 2016). For instance, at low Ca and lowM conditions, the maximum S nw can be less
than 50%. However, for the sake of consistency, in this study we
plotted all permeability curves for S nw values of
10% to 90%.
2.3.2 M –Ca –permeability color map
After confirming the influence of M and Ca onk nw and k w in a two-phase
flow system, we conducted simulations at various M and Caconditions with S nw held at a constant value of
20% to create a plot of M –Ca –k nwcorrelation or color diagram. The 20% S nwcondition was chosen because it can realistically be achieved under allM –Ca conditions (Tsuji et al., 2016). The color diagram
is useful to evaluate the degree of influence of M and Caon changes of relative permeability when both parameters influence a
system, and can also provide estimates of relative permeability in a
reservoir under a wide range of M and Ca .