2. Mathematical model for forced imbibition in a two-dimensional pore doublet

In order to understand the mechanism underlying forced imbibition, we first consider a simple geometry known as the pore doublet model, which is sketched in figure 1. The pore doublet consists of three parts: a feeding channel CA that supplies the wetting fluid; two capillary tubes that bifurcate from point A and reunite downstream at point B; and an exit channel BD. The branch channel at the bottom (capillary 1) has a smaller width $2r_1$ and the one at the top (capillary 2) has a greater width $2r_2$. The two branches are symmetric with the same length of $L$ along the flow direction, and the angle between the horizontal line and the centreline of each branch channel is $45^{\circ }$. Initially, the entire pore doublet is saturated with a non-wetting fluid. The wetting fluid is injected from the left inlet at a given flow rate $q$, while a constant pressure is assumed at the right outlet. The feeding and exit channels are of equal widths $h=2(r_1+r_2)$, and a constant contact angle of $\theta =30^\circ$ is considered. In the following, we present a theoretical modelling of the imbibition process based on the aforementioned pore doublet.

Figure 1.

Schematic diagram of the imbibition process in a two-dimensional pore doublet ($r_2=2r_1$).

2.1. Governing equations

Assuming that the flow through the pore doublet is a steady laminar flow and the two-phase interface advances with a constant mean curvature, the pressure difference $\Delta p$ between point A and B includes the viscous pressure drop $\Delta p_{vis}$ and the capillary pressure drop $\Delta p_{cap}$, which can be written as

(2.1)\begin{gather} \Delta p=p_A-p_B=\Delta p_{vis,1}-\Delta p_{cap,1}=\frac{3q_1}{2r_1^3}[\eta_w L_1 +\eta_n(L-L_1)]-\frac{\sigma \cos \theta}{r_1}, \end{gather}

(2.2)\begin{gather} \Delta p=p_A-p_B=\Delta p_{vis,2}-\Delta p_{cap,2}=\frac{3q_2}{2r_2^3}[\eta_w L_2 +\eta_n(L-L_2)]-\frac{\sigma \cos \theta}{r_2}, \end{gather}

where $p_A$ and $p_B$ are the pressures at points A and B, respectively; $q_1$ and $q_2$ are the volumetric flow rates (m$^2$ s$^{-1}$) in capillaries 1 and 2 with $q_i=2r_i\cdot u_i$ ($i=1,2$), and $u_1$ and $u_2$ are the corresponding average velocities; $\eta _w$ and $\eta _n$ are the dynamic viscosities of the wetting and non-wetting fluids, respectively; $\sigma$ is the interfacial tension coefficient; and $L_1$ and $L_2$ are the lengths that are occupied by the wetting fluid in the small and large capillaries. From mass conservation, one can write the total volumetric flow rate $q$ as

(2.3)\begin{equation} q=q_1+q_2. \end{equation}

2.2. Non-dimensionalisation of governing equations

In order to non-dimensionalise the governing equations, the scaling parameters for length, time and pressure are introduced:

(2.4a–c)\begin{equation} l_s= r_1, \quad t_s=\frac{2r_1 L}{q}, \quad p_s=\frac{3\eta_n qL}{2r_1^3}, \end{equation}

where the subscript $s$ refers to scaling parameters. Denoting the width ratio of both branch channels as $k=r_2/r_1$, which also represents the square root of the permeability ratio in capillary 2 to capillary 1, one easily obtains $\hat {r}_1=1$ and $\hat {r}_2=k$, where the hat over a variable means that the variable is non-dimensional. Substituting (2.4a–c) into (2.1) and (2.2) leads to

(2.5)\begin{equation} {\Delta \hat{p}_{vis,i}}=\hat{q}_i\left[ \frac{\lambda \cdot \left(\dfrac{\hat{L}_i }{\hat{L}}\right)}{\hat{r}_i^3}+\frac{\left(1-\dfrac{\hat{L}_i}{\hat{L}}\right) }{\hat{r}_i^3}\right ],\quad i=1,2, \end{equation}

where $\hat {q}_i=\hat {r}_i\hat {u}_i/\hat {L}$, $\hat {u}_i=\mathrm {d} \hat {L}_i/\mathrm {d} \hat {t}$ and $\lambda =\eta _w/\eta _n$ is the viscosity ratio of wetting to non-wetting fluid. Similarly, the capillary pressure drop can be written as

(2.6)\begin{equation} \Delta \hat{p}_{cap,i}=\frac{1}{Ca_m}\cdot \frac{1}{\hat{r}_i},\quad i=1,2, \end{equation}

where $Ca_m=3\eta _n qL/(2r_1^2\sigma \cos \theta$) is the preferential capillary number. It is clear that this capillary number takes into account the influence of pore length and size, and is different from the standard one, which is defined by the inlet mean velocity $u_{in}$ as $Ca=u_{in}\eta _n/(\sigma \cos \theta )$. Specifically, its value could be 2–4 orders of magnitude higher than that of the standard capillary number. Combining (2.5) and (2.6), one can obtain the total pressure drop as

(2.7)\begin{align} \Delta \hat{p}&=\Delta \hat{p}_{vis,i}+\Delta \hat{p}_{cap,i}=\hat{r}_i\cdot \frac{\mathrm{d}\left(\dfrac{\hat{L}_i }{\hat{L}}\right) }{\mathrm{d} \hat{t}}\left[ \frac{\lambda \left(\dfrac{\hat{L}_i }{\hat{L}}\right)}{\hat{r}_i^3}+\frac{ \left(1-\dfrac{\hat{L}_i }{\hat{L}}\right) }{\hat{r}_i^3}\right]\nonumber\\ &\quad -\frac{1}{Ca_m}\cdot \frac{1}{\hat{r}_i},\quad i=1,2. \end{align}

Mass conservation can also be written in dimensionless form as

(2.8)\begin{equation} \frac{\mathrm{d} \left(\dfrac{\hat L_1}{\hat L}\right)}{\mathrm{d}\hat t}\cdot \hat{r}_1+\frac{\mathrm{d} \left(\dfrac{\hat L_2}{\hat L}\right)}{\mathrm{d}\hat t}\cdot \hat{r}_2=1. \end{equation}

To solve the interface movement, we write (2.7) for each daughter channel. Equating the resulting two equations gives

(2.9)\begin{align} &\hat{r}_1\cdot \frac{\mathrm{d}\left(\dfrac{\hat L_1}{\hat L}\right) }{\mathrm{d} \hat{t}}\left[ \frac{\lambda \dfrac{\hat{L}_1}{\hat L} }{\hat{r}_1^3}+\frac{ \left(1-\dfrac{\hat{L}_1}{\hat L}\right) }{\hat{r}_1^3}\right ]-\frac{1}{Ca_m}\cdot \frac{1}{\hat{r}_1}\nonumber\\ &\quad =\hat{r}_2\cdot \frac{\mathrm{d}\left(\dfrac{\hat{L}_2}{\hat L}\right)}{\mathrm{d} \hat{t}}\left[ \frac{\lambda \dfrac{\hat{L}_2}{\hat L} }{\hat{r}_2^3}+\frac{ \left(1-\dfrac{\hat{L}_2}{\hat L}\right) }{\hat{r}_2^3}\right ]-\frac{1}{Ca_m}\cdot \frac{1}{\hat{r}_2}. \end{align}

Substituting (2.8) into (2.9), we obtain an ordinary differential equation for $\hat {L}_i(t)$, i.e.

(2.10)\begin{equation} \frac{\mathrm{d}\left(\dfrac{\hat{L}_1}{\hat L}\right) }{\mathrm{d} \hat{t}}=\frac{\dfrac{1}{Ca_m}\cdot \left(\dfrac{1}{\hat{r}_1}-\dfrac{1}{\hat{r}_2}\right)+\phi \left(\dfrac{\hat{L}_2}{\hat L}\right)}{\hat{r}_1 \left[\phi \left(\dfrac{\hat{L}_1}{\hat L}\right)+\phi \left(\dfrac{\hat{L}_2}{\hat L}\right)\right]}, \end{equation}

where $\phi (\hat {L}_i/\hat L)= [\lambda \hat {L}_i/\hat L+(1-\hat {L}_i/\hat L)]/\hat {r}_i^3$. The ordinary differential equation for $\hat {L}_2(t)$ can be obtained by exchanging subscripts 1 and 2:

(2.11)\begin{equation} \frac{\mathrm{d}\left(\dfrac{\hat{L}_2}{\hat L}\right) }{\mathrm{d} \hat{t}}=\frac{\dfrac{1}{Ca_m}\cdot \left(\dfrac{1}{\hat{r}_2}-\dfrac{1}{\hat{r}_1}\right)+\phi\left(\dfrac{\hat{L}_1}{\hat L}\right)}{\hat{r}_2 \left[\phi\left(\dfrac{\hat{L}_1}{\hat L}\right)+\phi\left(\dfrac{\hat{L}_2}{\hat L}\right)\right]}. \end{equation}

To prevent the flow in the branch channels from moving backwards, the following constraints must be satisfied (Al-Housseiny et al. Reference Al-Housseiny, Hernandez and Stone2014):

(2.12)\begin{equation} 0 \leq \frac{\mathrm{d} \left(\dfrac{\hat{L}_i}{\hat L}\right)}{\mathrm{d} \hat{t}} \leq \frac{1}{\hat{r}_i}, \quad i = 1, 2. \end{equation}

2.3. Semi-analytical solutions

We first consider a pore doublet geometry with $k=2$, $\hat {L}=63.11$ and $\theta =30^\circ$. To obtain the location of the meniscus in each capillary, we numerically solve (2.10) and (2.11) subject to the constraint (2.12) using the first-order forward difference scheme for different values of $Ca_m$ and $\lambda$. Solutions are found with the initial condition that $[\hat {L}_1,\hat {L}_2]=[0,0]$ at $\hat {t}=0$.

Numerical results for several typical capillary numbers at $\lambda =0.025$, 1 and 20 are shown in figure 2, where the penetration lengths $\hat {L}_1$ and $\hat {L}_2$ are plotted as a function of time $\hat {t}$, normalised by the breakthrough time $\hat {t}_B$. For each viscosity ratio, at low $Ca_m$ (figure 2a,d,g), we can see that $\hat {L}_1=\hat {L}>\hat {L}_2$ when the breakthrough occurs, although $\hat {L}_1$ lags behind $\hat {L}_2$ until $\hat {t}/\hat {t}_B=0.96$ in figure 2(a); at high $Ca_m$ (figure 2c,f,i), the meniscus in channel 2 breaks through first, i.e. $\hat {L}_1<\hat {L}_2=\hat {L}$. This suggests that there exists a critical value of $Ca_m$ between low and high $Ca_m$, known as the critical preferential capillary number ($Ca_{m,c}$), at which the breakthrough of wetting fluid occurs simultaneously in both branch channels. In the case of simultaneous breakthrough, it can be easily obtained that two menisci break through at $t_B =2L(r_1+r_2)/q$, or $\hat {t}_B=(1+k)\hat {L}$ in the dimensionless form. As shown in figure 2(b,e,h), the values of $Ca_{m,c}$ are 4.02, 2 and 0.2 for the viscosity ratios of 0.025, 1 and 20. Clearly, the critical preferential capillary number is strongly dependent on the viscosity ratio. In addition, we also interestingly find that for $\lambda =1$ in figure 2(e), the imbibition rates are constant and exactly the same in both branch channels, and $Ca_{m,c}=k=2$, consistent with the theoretical prediction as shown in Appendix A.

Figure 2.

The lengths of the wetting fluid in the branch channels as a function of time obtained by solving (2.10) and (2.11) at $\lambda =0.025$ for (a) $Ca_m=3.64$, (b) $Ca_m=4.02$ and (c) $Ca_m=5.83$; at $\lambda =1$ for (d) $Ca_m=0.7$, (e) $Ca_m=2$ and (f) $Ca_m=5.25$; at $\lambda =20$ for (g) $Ca_m=0.146$, (h) $Ca_m=0.2$ and (i) $Ca_m=0.488$.

Different imbibition behaviours at low and high values of $Ca_m$ are attributed to the competition between the capillary pressure and the viscous resistance. At low flow rates (small $Ca_m$), the viscous resistance is negligibly small while the capillary pressure is dominant, which acts a driving force for the wetting fluid to progress; since the capillary pressure is inversely proportional to the channel width, the penetration length in channel 1 is larger than that in channel 2, i.e. $\hat L_1>\hat L_2$, at breakthrough. However, at high flow rates, the viscous force is dominant; because of the lower viscous resistance in channel 2, the penetration length in channel 2 would be greater than in channel 1, i.e. $\hat L_1<\hat L_2$.

To understand the effect of the viscosity ratio on the imbibition process, a theoretical analysis is then conducted for a wide range of viscosity ratios, varying from $10^{-4}$ to $10^{3}$. Figure 3 depicts the imbibition preference of the meniscus at breakthrough in the $\lambda$–$Ca_m$ diagram. Three typical regions are identified due to the competition between capillary and viscous forces: (I) the region below the solid blue line, where the meniscus in channel 1 outpaces that in channel 2 at breakthrough, i.e. $\hat L_1=\hat L>\hat L_2$; (II) the region above the solid blue line, where the meniscus in channel 2 outpaces that in channel 1 at breakthrough, i.e. $\hat L_1<\hat L_2=\hat L$; and (III) the border of regions (I) and (II), on which the menisci in channels 1 and 2 arrive at the downstream junction at the same time, i.e. $\hat L_1=\hat L_2=\hat L$ at breakthrough. It is noted that the border corresponds to the critical curve of $Ca_m$, i.e. the $Ca_{m,c}$ curve. We can observe that for $\lambda \geq 10$, the critical capillary number $Ca_{m,c}$ obeys a scaling relation $Ca_{m,c}=3.314\lambda ^{-1}$; whereas for $\lambda \leq 0.1$, it tends to converge to a value of around 3.9. Through figure 3, we are able to predict the filling order of the wetting fluid for varying viscosity ratio and $Ca_m$ in a pore doublet. In a previous work (Sorbie et al. Reference Sorbie, Wu and McDougall1995), the existence of critical parameters for characterising the simultaneous filling of both branch channels was discussed in terms of the aspect ratio ($r_i/L$) and the channel width ratio, and the influence of aspect ratio was explained as a result of the fluid inertia; however, the aspect ratio is incorporated into the definition of the preferential capillary number in the present study. In addition to the viscosity ratio, we also vary the value of $k$ from 0.1 to 10, and the resulting $Ca_{m,c}$ surfaces viewed from two different angles are shown in figure 4. It is interestingly observed that for any given viscosity ratio, $\log Ca_{m,c}$ increases linearly with $\log k$, and in particular $Ca_{m,c}=k$ for $\lambda =1$, consistent with the theoretical derivation in Appendix A. Moreover, as in the case of $k=2$, $Ca_{m,c}$ first remains nearly constant and then decreases with increasing $\lambda$ for any fixed $k$.

Figure 3.

The $\lambda$–$Ca_m$ diagram showing the imbibition preference in a two-dimensional pore doublet. The various symbols represent the cases where simultaneous breakthrough occurs. Connecting these symbols gives the blue solid line which divides the plane into two regions, i.e. (I) and (II). In (III), the meniscus first breaks through channel 1, whereas in (II) the breakthrough first occurs in channel 2. The border on which $\hat {L}_1=\hat {L}_2=\hat {L}$ at breakthrough is denoted as (III), and it follows a scaling relation $Ca_{m,c}=3.314\lambda ^{-1}$ for $\lambda \geq 10$. The dashed line is added to show the proportional relationship between $Ca_{m,c}$ and $\lambda ^{-1}$.

Figure 4.

The $\lambda$–$k$–$Ca_m$ diagram showing the $Ca_{m,c}$ surfaces viewed from two different angles.

2.4. Comparison between LBM simulations and semi-analytical solutions

In this section, the colour-gradient LBM (see Appendix B for details) is used to simulate the imbibition behaviour in a pore doublet and its capability is assessed by comparing with the semi-analytical solutions in § 2.3. The simulations are run in a $1575\times 409$ lattice domain with $r_1=15$ lattices and $r_2=30$ lattices, which are found fine enough to produce grid-independent results. Figure 5 shows the simulation results corresponding to the same values of $Ca_m$ and $\lambda$ as in figure 2. It is clear that the simulation results at breakthrough agree well with the semi-analytical solutions qualitatively, and for each $\lambda$, two menisci in branch channels are found to arrive at the downstream junction simultaneously at $Ca_{m,c}$, consistent with the semi-analytical predictions in figure 2 as well. Although figure 5(b,e,h) appears to be the same at breakthrough, it exhibits different evolution scenarios over time, which can be seen in figure 6. For example, at $\hat {t} = 0.5\hat {t}_B$, the meniscus in capillary 1 lags behind that in capillary 2 in figure 6(a) while the result is the converse in figure 6(c), which agree with the semi-analytical predictions in figure 2(b,h).

Figure 5.

Fluid distributions at breakthrough obtained from the LBM simulations for the same parameters as those in figure 2. Specifically, the first row: $\lambda =0.025$ with (a) $Ca_m=3.64$, (b) $Ca_m=4.02$ and (c) $Ca_m=5.83$; the second row: $\lambda =1$ with (d) $Ca_m=0.7$, (e) $Ca_m=2$ and (f) $Ca_m=5.25$; the third row: $\lambda =20$ with (g) $Ca_m=0.146$, (h) $Ca_m=0.2$ and (i) $Ca_m=0.488$. The non-wetting and wetting fluids are shown in red and blue, respectively.

Figure 6.

Fluid distributions at $\hat {t}=0.5\hat {t}_B$ obtained from the LBM simulations corresponding to the parameters used in figure 5(b,e,h). Specifically, these parameters are: (a) $\lambda =0.025$, $Ca_m=4.02$, (b) $\lambda =1$, $Ca_m=2$ and (c) $\lambda =20$, $Ca_m=0.2$. The non-wetting and wetting fluids are shown in red and blue, respectively.

To assess the transient behaviour, as an example, we present snapshots of the imbibition process for $Ca_m=3.64$ and $\lambda =0.025$ in figure 7, where the upper and lower rows represent the simulation results and the semi-analytical predictions, respectively. Again, good agreement between the simulation results and semi-analytical predictions is obtained, although the LBM simulation slightly overestimates $\hat L_1$ in figure 7(e). Having verified the colour-gradient LBM, we use it to investigate the imbibition displacement in a dual-permeability pore network in the next section, where the theoretical predictions are not applicable due to the inherent complex geometry.

Figure 7.

Snapshots of the imbibition process obtained by the LBM simulations (top) and the semi-analytical solutions (bottom) for $Ca_m=3.64$ and $\lambda =0.025$ at (a) $\hat {t}/\hat {t}_B=0$, (b) $\hat {t}/\hat {t}_B=0.19$, (c) $\hat {t}/\hat {t}_B=0.5$, (d) $\hat {t}/\hat {t}_B=0.69$, (e) $\hat {t}/\hat {t}_B=0.88$ and (f) $\hat {t}/\hat {t}_B=1.0$. In the top images, the non-wetting and wetting fluids are shown in red and blue, respectively. In the bottom images, the non-wetting and wetting fluids are shown in white and black, respectively.

3. Forced imbibition in a dual-permeability pore network

In this section, we first describe the geometry set-up of the problem along with the boundary conditions. Then the simulation results of imbibition displacement in the pore network are presented and compared with those previously obtained from the pore doublet.

As shown in figure 8(a), the porous media geometry used in this study consists of an inlet and an outlet section, connected by a pore network. The pore network includes two distinct permeability zones with each occupying approximately a half-width of the domain. Each homogeneous zone contains a staggered periodic array of uniform circular grains (see figure 8b). We run the simulations in a $1428\times 1441$ lattice domain, which corresponds to a physical size of $0.714\times 0.721$ cm$^2$. The length of the pore network is $1078$ lattices. The diameter of solid grains is 64 lattices in the high-permeability zone and 32 lattices in the low-permeability zone. The diameter of pore bodies in the high-permeability (low-permeability) zone is 56 (28) lattices, and the corresponding pore throat width is 20.8 (10.4) lattices. Both permeability zones have equal porosity of $\epsilon =0.55$. Initially, the pore network is saturated with the non-wetting (red) fluid, and the wetting (blue) fluid is injected from the left-hand inlet continuously with a parabolic velocity profile of ${\boldsymbol u}=(6u_{in}({(y(H-y))}/{H^2}),0)$, which is imposed by the velocity boundary scheme of Zou & He (Reference Zou and He1997). Here, $H$ is the inlet width of the porous media geometry and $u_{in}$ is the inlet mean velocity. A constant pressure is set at the right-hand outlet via the pressure boundary condition of Zou & He (Reference Zou and He1997), and the top and bottom boundaries are no-slip walls. The densities of the two fluids are assumed to be equal since the displacement mainly occurs in the horizontal direction, where the effect of gravity can be negligible. Each simulation is run until the wetting fluid breaks through the right-hand boundary of the pore network.

Figure 8.

(a) The initial fluid distribution and set-up of the boundary conditions for the imbibition simulations in the dual-permeability porous geometry. The white circles represent the solid grains, while the blue and red regions represent the wetting and non-wetting fluids, respectively. The whole computational domain has a size of $1428\times 1441$ lattices, which consists of an inlet and an outlet section, connected by a pore network. (b) Representation for staggered array of circular grains in the pore network. A pore body is defined by the largest circle fitting locally the pore space. The size of a pore throat is defined by the narrowest width between two nearest solid grains.

We first consider a viscosity ratio of 0.1 for various values of $Ca_m$, where $Ca_m$ is defined by $Ca_m=3\eta _n qL/(2r_1^2\sigma \cos \theta$), with $r_1$ and $L$ taken as the average of pore body radius and half-throat width (see figure 8b) in the low-permeability zone and the length of the pore network. Figure 9 shows the corresponding fluid distributions in the dual-permeability pore network at breakthrough. It is found that at low (high) values of $Ca_m$, the wetting fluid prefers to invade the low-permeability (high-permeability) zone and the breakthrough first occurs in the low-permeability (high-permeability) zone, consistent with the previous observations for the pore doublet. In all cases, very few drops of the non-wetting fluid are trapped as the residual phase in the imbibition process as the wetting fluid progresses. For figure 9(a–e) ($Ca_m=1.4574\text {--}72.8723$), we notice an oblique advancing pattern of the wetting fluid in both high- and low-permeability zones, but this phenomenon disappears when $Ca_m$ is increased up to $Ca_m=145.7446$ (figure 9f). This suggests that the oblique advancing of the wetting fluid arises from the non-trivial interfacial tension.

Figure 9.

Fluid distributions in the dual-permeability pore network at breakthrough for (a) $Ca_m=1.4574$, (b) $Ca_m=2.9149$, (c) $Ca_m=4.3723$, (d) $Ca_m=14.5745$, (e) $Ca_m=72.8723$ and (f) $Ca_m=145.7446$. The viscosity ratio of wetting to non-wetting fluids is fixed at 0.1. The non-wetting and wetting fluids are shown in red and blue, respectively.

To better understand the oblique advancing pattern, as an example, we plot the evolution of fluid distributions during the early imbibition at $Ca_m=1.4574$, which is shown in figure 10. It is known that the dominant capillary pressure is larger in the smaller pores and throats according to the Young–Laplace equation, so the smallest pores and throats are filled first. For the present grain arrangement, let us take a close look at the interface between two vertically aligned solid grains A and B, as shown in figure 8(b). It is seen that a flat interface (represented by the blue solid line) with zero capillary pressure is able to touch the solid grain C, and thus the advancing meniscus of the wetting fluid always progresses towards the next column of grains through a triangle shape, as marked by the black triangles in figure 10. In addition, as shown in figure 10(f), as the wetting fluid invades the region marked by the black triangle, it cannot infiltrate in the direction highlighted by the dashed arrow due to the requirement of a positive pressure difference between the wetting and non-wetting fluids to overcome the capillary valve resistance (Xu et al. Reference Xu, Liu and Valocchi2017), but progress towards the direction highlighted by the solid arrow due to the merging with the neighbouring interface. As a result, the wetting fluid penetrates layer by layer along the direction indicated by the solid arrow, forming an oblique advancing pattern. A similar process occurs in the high-permeability zone, but in a direction perpendicular to the invading direction in the low-permeability zone. On the other hand, at the highest $Ca_m$ in figure 9(f), the aforementioned pore filling order is disrupted and no longer applicable, as here the viscous force dominates the imbibition behaviour.

Figure 10.

Fluid distributions in the porous media geometry for $Ca_m=1.4574$ and $\lambda =0.1$ at (a) $\hat t/\hat {t}_B=0.0145$, (b) $\hat t/\hat {t}_B=0.0217$, (c) $\hat t/\hat {t}_B=0.029$, (d) $\hat t/\hat {t}_B=0.0362$, (e) $\hat t/\hat {t}_B=0.0434$ and (f) $\hat t/\hat {t}_B=0.0507$. Each inset shows a close-up view of the region indicated by the black rectangular box in the lower left corner.

We then study the effect of viscosity ratio on the imbibition preference. A wide range of viscosity ratios, varying from $\lambda =0.02$ to 50.0, is considered. For each viscosity ratio, at least three different values of $Ca_m$ are simulated, covering three typical patterns observed at breakthrough. The saturation data at breakthrough for various viscosity ratios and capillary numbers are listed in table 1, where $S_1$, $S_2$ and $S_w$ are the wetting fluid saturations in the low-permeability zone, the high-permeability zone and the entire pore network. Among all the cases considered, the maximum imbibition efficiency is obtained under the conditions of $\lambda =50$ and $Ca_m=0.08745$, where the wetting fluid saturations in both permeability zones are roughly the same (the corresponding values $S_1=0.8456$ and $S_2=0.8645$). In addition, for each viscosity ratio, the highest imbibition efficiency is always achieved when $S_1$ is closest to $S_2$. This implies that the critical capillary numbers $Ca_{m,c}$ represent the optimal condition to improve the imbibition efficiency.

Table 1.

Saturations $S_1$, $S_2$ and $S_w$ at breakthrough for various values of viscosity ratio ($\lambda$) and capillary number ($Ca_m$), where $S_{1}$ and $S_{2}$ are the wetting fluid saturations in the low- and high-permeability zone and $S_w$ is the wetting fluid saturation in the whole pore network.

To locate the values of $Ca_{m,c}$ for different viscosity ratios, we extract the data regarding the imbibition preference from table 1 and plot them in the $\lambda$–$Ca_{m}$ diagram, as shown in figure 11. In this figure, the open symbols represent the cases where $S_1>S_2$, while the filled symbols represent the cases where $S_1 < S_2$. This means that for each value of $\lambda$, the critical capillary number $Ca_{m,c}$ lies between two nearest open and filled symbols. As such, the $Ca_{m,c}$ curve can be approximately obtained, which is represented by the green solid line. For the sake of comparison, figure 11 also plots the $Ca_{m,c}$ curve from the pore doublet model (represented by the blue dashed line and directly taken from figure 3). It is clear that the present $Ca_{m,c}$ curve overlaps well with the one from the pore doublet model. As shown in Appendix C and Appendix D, we also vary the permeability ratio ($k^2$) and the contact angle ($\theta$), and find that the results overall agree with the predictions from the pore doublet model. All of the results suggest that the simplified pore doublet model can provide insights into the physics of immiscible displacement in the more complex dual-permeability pore network.

Figure 11.

The $\lambda$–$Ca_m$ diagram showing preferential imbibition in a dual-permeability pore network. The open symbols represent the cases where $S_{1}>S_{2}$ at breakthrough (I), while the filled symbols represent the cases where $S_{1} < S_{2}$ at breakthrough (II). Two images of fluid distributions are shown as inserts depicting regions (I) and (II). The green solid line represents the $Ca_{m,c}$ curve, on which $S_1=S_2$ at breakthrough. The $Ca_{m,c}$ curve (represented by the blue dashed line) from the pore doublet model is also plotted for comparison.

Although the pore doublet model can predict the variation of $Ca_{m,c}$ with $\lambda$ in a dual-permeability pore network, it is not clear whether the transient imbibition behaviour in the dual-permeability pore network can be correctly captured by the pore doublet model. In order to clarify this, we plot the time evolution of $S_1$ and $S_2$ (normalised by their maximum value at breakthrough) at three typical viscosity ratios in figure 12, where the semi-analytical solutions $\hat {L}_1$ and $\hat {L}_2$ (normalised by $\hat {L}$), obtained from (2.10) and (2.11) with the dimensionless numbers $Ca_m$ and $\lambda$ identical to those in the pore network, are also shown for comparison. For each viscosity ratio, the agreement between the LBM results and the semi-analytical solutions is generally better at higher $Ca_m$ where $S_1 < S_2$, but worse when $S_1>S_2$ where the interfacial tension is dominant. The larger discrepancy when $S_1>S_2$ (see figure 12a,c,e) is attributed to the fact that in the dual-permeability pore network, the interface varies and thus the capillary pressure varies when the meniscus moves from the throat to the pore body or from the pore body to the throat, while the capillary pressure remains constant in the pore doublet. We also increase the porosity to 0.62 in a similar pore network, and again observe a good agreement between the saturation evolution and the semi-analytical predictions, which is shown in Appendix E. In addition, we interestingly notice in figure 12(c) that after $\hat t/\hat {t}_B = 0.5$, the wetting fluid infiltrates into the high- and low-permeability zones alternately. Figure 13 shows the corresponding snapshots, from which it is seen that the wetting fluid invades only into the high-permeability zone in figure 13(a–c) and only into the low-permeability zone in figure 13(d–f).

Figure 12.

The saturations in the low- and high-permeability regions (normalised by their maximum value at breakthrough) as a function of time in the dual-permeability pore network at $\lambda =0.025$ for (a) $Ca_m=2.9149$ and (b) $Ca_m=5.8298$; at $\lambda =1$ for (c) $Ca_m=0.4372$ and (d) $Ca_m=7.2872$; at $\lambda =20.0$ for (e) $Ca_m=0.1457$ and (f) $Ca_m=0.2186$. The semi-analytical solutions from the pore doublet model at the same values of $\lambda$ and $Ca_m$ are also shown for the comparison.

Figure 13.

Snapshots of the imbibition for $Ca_m=0.4372$ and $\lambda =1$ at (a) $\hat t/\hat {t}_B=0.5106$, (b) $\hat t/\hat {t}_B=0.5532$, (c) $\hat t/\hat {t}_B=0.5957$, (d) $\hat t/\hat {t}_B=0.7234$, (e) $\hat t/\hat {t}_B=0.766$ and (f) $\hat t/\hat {t}_B=0.8085$. The snapshots from (a) to (f) correspond to the solid dots marked by A to F in figure 12(c).