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.

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 CO₂saturation. 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/m³ 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).

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%.

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 .