Lattice Boltzmann Modeling of Spontaneous Imbibition in Variable-Diameter Capillaries

Abstract

Previous micro-scale studies of the effect of pore structure on spontaneous imbibition are mainly limited to invariable-diameter capillaries. However, in real oil and gas reservoir formations, the capillary diameters are changing and interconnected. Applying the lattice Boltzmann color gradient two-phase flow model and the parallel computation of CPUs, we simulated the spontaneous imbibition in variable-diameter capillaries. We explored the reasons for the nonwetting phase snap-off and systematically studied the critical conditions for the snap-off in spontaneous imbibition. The effects of pore-throat aspect ratio, throat diameter, and the pore-throat tortuosity of the capillary on spontaneous imbibition were studied. Through analyzing the simulated results, we found that the variation in the capillary diameter produces an additional resistance, which increases with the increase in the pore-throat ratio and the pore-throat tortuosity of a capillary. Under the action of this additional resistance, the snap-off phenomenon sometimes occurs in the spontaneous imbibition, which makes the recovery efficiency of the non-wetting phase extremely low. In addition, the main factors affecting this phenomenon are the pore-throat ratio and the pore-throat tortuosity, which is different from the conventional concept of tortuosity. When the snap-off does not occur, the spontaneous imbibition velocity increases when the throat diameter increases and the pore-throat aspect ratio is fixed, and when the period increases, i.e., the diameter changing rate decreases, the spontaneous imbibition velocity also increases. In addition, when the capillary throat diameter is fixed, a bigger pore diameter and a smaller period of sine function both inhibit the speed of spontaneous imbibition.

1. Introduction

As a result of the increase in oil consumption and the depletion of conventional oil production, tight and shale oil and gas have become important supplementary resources of conventional oil and gas energy [1,2]. Tight and shale reservoirs are mainly developed through huge volumetric fracturing of horizontal wells, which creates complex high-conductance fracture networks in the matrix [3]. Spontaneous imbibition of water between the fracture and matrix is considered to be a fundamental mechanism for the oil recovery of tight and shale reservoirs [4,5,6,7]. Spontaneous imbibition is defined as the capillary displacement of a non-wetting phase by a wetting phase without the assistance of outside forces [8]. Regarding spontaneous imbibition, experimental studies are mainly based on the mass conservative method [9]. However, due to the oil-wall-hanging phenomenon [9], the threshold jump [10,11], the Malangoni effect, and other typical problems that seriously affect the measurement accuracy of traditional experimental methods [12], it is impossible to conduct in-depth research and provide a microscopic explanation for the mechanism of spontaneous imbibition. Therefore, most previous studies have focused on the elementary effects, including matrix contacting areas’ geometric properties (e.g., shape and size), rock physical properties (e.g., porosity and permeability) [13], fluid properties (e.g., density, viscosity, and interfacial tension), and rock matrix wettability [14,15,16]. Recent analytical studies focused on the effects of microscopic pore structure, for example, based on the modified Hagen–Poiseuille and Laplace–Young equations. Cai et al. derived a generalized model of porous media imbibition that takes into account the influence of curved capillary and non-circular cross-section shapes [17]; Zheng et al. simulated the spontaneous imbibition process in an irregular channel model, and analyzed the influence of tortuosity and channel shape on spontaneous imbibition behavior and the competing interfaces’ advancing behavior [18]; Liu and Cai studied the effect of micro-fractures on spontaneous imbibition at the pore scale [19]. At present, the microscopic study of pore structure is mostly limited to straight or tortuous tubes with no change in the pore diameter of a single tube. A pore structure without variation in the pore size is not consistent with the characteristics of real reservoir porous media, and, therefore, the simulated results of the invariant-diameter capillaries cannot well reflect the spontaneous imbibition characteristics of the real reservoir formations.

Additionally, when two-phase fluid leads to spontaneous imbibition in the porous media, such as the rock, the wetting fluid moves along the wetting film, and may accumulate in roughness or smaller spaces (throats) under the action of capillary force [20]. When the wetting phase loses contact with the solid, the unstable two-phase interface may lead to the rupture of the non-wetting phase bridge, and the non-wetting fluid becomes trapped in the larger space (pore), which is referred to as snap-off [21,22], as shown in Figure 1. To the best of our knowledge, limited studies have been conducted of the snap-off during the spontaneous imbibition process.

Pore-scale simulation methods make it possible to simulate the spontaneous imbibition in a detailed 3D pore space; such methods include the smooth particle fluid dynamics method (SPH), lattice Boltzmann method (LBM), pore network simulation, and molecular dynamics simulation. Specifically, the SPH boundary is not precise enough because the integral domain is truncated at the boundary [23,24]. The pore network model extracts the complex porous media into the ideal regular geometric-shaped pore networks, and this method can greatly improve the calculation efficiency. However, its disadvantage is that it limits the calculation ability and accuracy, and cannot reflect the detailed irregular pore structure [25,26,27]. Moreover, the molecular dynamics simulation is often used at the nanoscale and at a very small time interval, and is too computationally demanding to work on larger scales [28,29,30,31]. The lattice Boltzmann method (LBM) is a mesoscopic simulation method [32]. Compared with other methods, this method uses a mesoscopic model characteristic that is between a macro continuous model and a micro molecular dynamics model. It can not only ignore the details of molecular motion, like in the macro method, but also has the characteristics of the wide applicability of the micro method. Therefore, it has the advantages of a simple description of fluid interaction, easy setting of complex boundaries, easy parallel calculation, and easy implementation of programs. Commonly used multiphase flow models include the color gradient model [33,34], the Shan–Chen model [35,36,37], and the free energy model [38,39]. The color gradient model is the first proposed lattice Boltzmann model of multiphase flow [40]. In this model, two distribution functions are used to describe two kinds of fluids. In addition to the conventional streaming and collision steps, a disturbance operator is added to generate interfacial tension, and a recolor process is added to realize phase separation and maintain the phase interface. The Shan–Chen model cannot simulate the two-phase flow with high viscosity contrast, and it is difficult in the free energy model to set the wettability. By comparison, the color gradient model can simulate the high viscosity ratio, and the contact angles can easily be set by specifying the density of the two fluids in the solid nodes without changing the interfacial tension. Because of this simplicity, the color gradient model has an advantage for two-phase flow simulation in porous media.

In this work, we use the LBM color gradient model to simulate the spontaneous imbibition of variable-diameter capillary tubes and study the microscopic flow mechanism with the help of the MPI parallel implementation of CPUs. Firstly, the description of the color gradient model and the parallel computational method are presented. Then, the model is validated by comparing the simulated results with the classical LW equations for different viscosities, contact angles, and interfacial tensions [41]. This is followed by the simulation study of the variable-diameter capillary spontaneous imbibition processes, including the setting of the variable-diameter capillary model and its key parameters. Finally, the simulated results are carefully analyzed, with particular focus on the conditions underlying the snap-off phenomenon and the advancing velocity of the wetting phase front and its main influencing factors.

2. Methodologies

2.1. Lattice Boltzmann Color Gradient Two-Phase Flow Model

We apply the color gradient model of the lattice Boltzmann (RK model) [33], which can simulate the pore-level two-phase flow of immiscible fluids in complex porous media, where fluids are allowed to have different densities and viscosities. The overall motion of particles is represented by a distribution function, and the space where particles exist is represented by a lattice arrangement [42]. We use the D2Q9 lattice scheme (D = 2 dimensions, and Q = 9 velocities on the discrete lattice), as shown in Figure 2.

For two-phase flows, the RK model can be summarized in three steps: collision, recoloring, and streaming. $f_i^k$ represents the particle distribution function (PDF) for fluid k (k = r or b, where r and b indicate red and black fluids), and the total PDF at $\left(x,t\right)$ is $f_i\left(x,t\right)={\displaystyle {\sum }_k{f}_i^k\left(x,t\right)}$. $i$ represents the discrete velocity components, $i$ = 0, 1, 2, …, 8. Commonly, the lattice spacing and time steps are set as $\Delta x=\Delta t=1$.

The collision step Is given as [33]:

$$f_i^{k*}\left(x,t\right)=f_i^k\left(x,t\right)+{\left({\Omega }_i^k\right)}^1+{\left({\Omega }_i^k\right)}^2,$$

(1)

where ${\left({\Omega }_i^k\right)}^1$ and ${\left({\Omega }_i^k\right)}^2$ are collision operators. The first collision operator ${\left({\Omega }_i^k\right)}^1$ adopts the BGK collision form, as formulated by:

$${\left({\Omega }_i^k\right)}^1=-\frac{\Delta t}{\tau }\left(f_i^k\left(x,t\right)-f_i^{k,eq}\left(x,t\right)\right),$$

(2)

where $\tau$ and $f_i^{k,eq}\left(x,t\right)$ are respectively, the relaxation time and the equilibrium distribution function, and $f_i^{k,eq}\left(x,t\right)$ is defined as [43,44]:

$$f_i^{k,eq}\left(x,t\right)={\rho }_k\left(c_i+w_i\left(\frac{e_i\cdot u}{c_s^2}+\frac{{\left(e_i\cdot u\right)}^2}{c_s^4}-\frac{{\left(u\right)}^2}{2c_s^2}\right)\right)$$

(3)

where ${\rho }_k$ and $u$ are the density of each fluid and the momentum at $\left(x,t\right)$, respectively, and can be expressed in terms of distribution functions as follows:

$${\rho }_k={\displaystyle \sum _i^{}{f}_i^k},$$

(4)

$$\rho u={\displaystyle \sum _k^{}{\displaystyle \sum _i^{}{f}_i^k}}{e}_i,$$

(5)

where $\rho$ is the macroscopic density, and is defined as:

$$\rho ={\displaystyle \sum _k^{}{\rho }_k},$$

(6)

and the other parameters are:

$$c_i=\left\{\begin{array}{ll}{\alpha }_{\kappa }\hfill & i=0\hfill \\ \frac{1-{\alpha }_{\kappa }}{5}\hfill & i=1,2,3,4\hfill \\ \frac{1-{\alpha }_{\kappa }}{20}\hfill & i=5,6,8,9\hfill \end{array}\right.,$$

(7)

$c_i$, ${\alpha }_{\kappa }$ are constants.

$$w_i=\left\{\begin{array}{ll}\frac{4}{9}\hfill & i=0\hfill \\ \frac{1}{9}\hfill & i=1,2,3,4\hfill \\ \frac{1}{36}\hfill & i=5,6,7,8\hfill \end{array}\right..$$

(8)

When the viscosities (relaxation times) for the two fluids are different, $\mathit{\tau }$ can be written as [34]:

$$\tau =\left\{\begin{array}{ll}{\tau }_r\hfill & \frac{{\rho }_r-{\rho }_b}{{\rho }_r+{\rho }_b}>\delta \hfill \\ s_1+s_2\cdot \frac{{\rho }_r-{\rho }_b}{{\rho }_r+{\rho }_b}+s_3\cdot {\left(\frac{{\rho }_r-{\rho }_b}{{\rho }_r+{\rho }_b}\right)}^2\hfill & 0<\frac{{\rho }_r-{\rho }_b}{{\rho }_r+{\rho }_b}\le \delta \hfill \\ t_1+t_2\cdot \frac{{\rho }_r-{\rho }_b}{{\rho }_r+{\rho }_b}+t_3\cdot {\left(\frac{{\rho }_r-{\rho }_b}{{\rho }_r+{\rho }_b}\right)}^2\hfill & -\delta \le \frac{{\rho }_r-{\rho }_b}{{\rho }_r+{\rho }_b}\le 0\hfill \\ {\tau }_b\hfill & \frac{{\rho }_r-{\rho }_b}{{\rho }_r+{\rho }_b}<-\delta \hfill \end{array},\right.$$

(9)

$\delta$, $s$ and $t$ are constants.

The second collision operator ${\left({\Omega }_i^k\right)}^2$ of Equation (1) is given by [44]:

$${\left({\Omega }_i^k\right)}^2=A\left|F\right|\left(w_i{\left(\cos \left({\lambda }_i\right)\left|c_i\right|\right)}^2-B_i\right),$$

(10)

where $F\left(x,t\right)$ is the color gradient, ${\lambda }_i$ is the angle between $F\left(x,t\right)$ and $c_i$, A is a parameter that controls the interfacial tension, and B is a specified constant. Recently, an optimized color gradient $F\left(x,t\right)$ was proposed by Mora et al. (2020) [45], as given by:

$$F\left(x,t\right)={\displaystyle \sum _i^{}{b}_i{c}_i}\left({\rho }_r\left(x+c_i\Delta t,t\right)-{\rho }_b\left(x+c_i\Delta t,t\right)\right),$$

(11)

where $b_i$ are scalar coefficients of the second-order finite difference approximation of the gradient given by:

$$b_i=\left\{\begin{array}{cc}\frac{1}{2+4w}& i=1,2,3,4\\ \frac{w}{2+4w}& i=5,6,7,8\end{array}\right..$$

(12)

After the collision step, the recoloring step should be calculated to achieve the separation of the two fluids, and is written as [33]:

$$f_i^{r+}=\frac{{\rho }_r}{\rho }{f}_i^{*}+\beta \frac{{\rho }_r{\rho }_b}{{\rho }^2}{f}_i^{eq}\left(\rho ,0\right)\cos \left({\lambda }_i\right),$$

(13)

$$f_i^{b+}=\frac{{\rho }_b}{\rho }{f}_i^{*}-\beta \frac{{\rho }_r{\rho }_b}{{\rho }^2}{f}_i^{eq}\left(\rho ,0\right)\cos \left({\lambda }_i\right).$$

(14)

Next, through the streaming step, the particle distribution function of the next time step is obtained [44]:

$$f_i^k\left(x+e_i\Delta t,t+\Delta t\right)=f_i^{k+}\left(x,t\right),$$

(15)

The three steps (collision, recoloring, and streaming) are repeated to realize the simulation of two-phase flow. In addition, the relationship between kinematic viscosity ${\nu }_k$ and relaxation time ${\tau }_k$ can be obtained as:

$${\nu }_k=c_s^2\left({\tau }_k-0.5\right)\Delta t,$$

(16)

where $c_s^{}=\Delta x/\left(\sqrt{3}\Delta t\right)$ is the speed of sound in the lattice, and the contact angle ${\theta }_w$ on the solid boundary can be analytically determined by [33]:

$${\theta }_w={\cos }^{-1}\left(\frac{{\rho }_{wr}-{\rho }_{wb}}{{\rho }_i}\right).$$

(17)

2.2. The Parallel Implementation of the Lattice Boltzmann Method

The study problem was computed using the Python programming. Although the LBM is simple to program, it requires a large amount of computation and takes a long time to run on a single personal computer. Parallel processing is required to increase the running speed. There are two approaches to parallel computing: CPUs in parallel; and GPUs in parallel. Here, we used the Message Passing Interface (MPI) [46] and cloud computing servers, which adopt CPU parallelism. There are also two ways to use MPI to run the LBM in parallel. One is to decompose the time domain and run this in parallel. However, the D2Q9 model has nine speed components, and the maximum number of parallel processes is nine, which is unfavorable to improving the running speed of the LBM. Second, spatial domain decomposition is adopted to decompose the calculation area into N blocks, so that N MPI processes can be used to calculate LBM, as shown in Figure 3a (n = 8), greatly reducing the calculation time.

The calculation of the color gradient model is divided into three steps: collision, recoloring, and streaming. The recoloring step is a spatially local operation and does not require information transfer. In the collision and streaming steps, since the data are passed from the lattice to the surrounding lattices, which moves the data past the region boundary, it needs to communicate with the adjacent domains. We achieve this by adding a virtual area in the analog domain, as shown in Figure 3b, which illustrates the concept of decomposition in the 4 × 2 domain. The red area is the real lattice (simulation area), and the blue area is the virtual lattice. We take the collision step as an example to illustrate the specific process of information transmission. Before the collision step, we pass the data on the real grid at the edge of the decomposed region to the virtual grid adjacent to the decomposed region, as shown in Figure 4. In the collision step, we flow ${\rho }_k$ from the virtual lattice to the real lattice. This requires a total of four communication steps: up, down, right, and left (for each decomposed domain), as shown in Figure 4. In this work, we take the spatial domain decomposition scheme and use 120 cores (120 parallel MPI processes) to carry out parallel computing for the established model, with good computational efficiency.

2.3. Model Verification

We verified the simulated the contact angles of three phases (solid, wetting fluid, and non-wetting fluid) and the classical LW equations of spontaneous imbibition for the proposed color gradient model [41].

In the color gradient model, we can set different three-phase contact angles by adjusting the density of the two-phase fluid on the solid boundary lattice, as shown in Equation (17). At the solid boundary, we set up five different sets of two-phase fluid densities to verify the accuracy of the color gradient model. In the simulation, the setup of the solid boundary and the initial distribution of the two-phase fluid is shown in Figure 5. The simulation results of different contact angles are shown in Figure 6, from which we can see that the simulation results correspond well with the contact angle obtained by Equation (17). The coordinates of Figure 5 and Figure 6 are position coordinates.

Then, considering different wetting angles and different two-phase viscosity ratios ($\frac{{\mu }_r}{{\mu }_b}$), the LW equation (Equation (18)) was used to verify the color gradient model, and the LW equation can be written as [18]:

$$\frac{dl}{dt}=\frac{H\sigma \cos \theta }{6\left({\mu }_b\left(L-l\right)+{\mu }_{r}l\right)},$$

(18)

where, respectively, $H$ and $L$ are the width and length of the model, $l$ is the position of the imbibing leading interface, ${\mu }_b$ and ${\mu }_r$ are the viscosity of the wetting phase fluid and the non-wetting phase fluid (Equation (16)), $\sigma$ and $\theta$ are the interfacial tension between the two-phase fluid (determined by A in Equation (10)) and the three-phase antennae (Equation (17)).

The simulated results are shown in Figure 7 and Figure 8. It can be seen from Figure 7 and Figure 8 that the matching is not perfect at the beginning of the simulation, but the simulation results in other stages fit well with the analytical solution of the LW equation. In the beginning, it takes some time for the bending interface of the spontaneous imbibition leading interface to form and stabilize, but the LW equation does not consider this transition, which leads to a tiny matching error.

3. Results with Discussion

3.1. The Setting of Variable-Diameter Capillaries

We constructed the variable-diameter capillary through the sine function defined as $A\sin \left(\omega x\right)$, as shown in Figure 9, and the diameter formula of the variable-diameter capillary is expressed as:

$$d=2A\sin (\frac{2\pi x}{T})+2A+d_t.$$

(19)

The wetting phase and the non-wetting phase are water and oil, respectively, and they both have a density of 1.

Non-equilibrium extrapolation boundary conditions were set at both ends of the model to ensure the same pressure at both ends. At the initial moment, the fluid is stationary, and its distribution is shown in Figure 9.

In Figure 10 and Equation (19), $T=\frac{2\pi }{\omega }=40\pi$, $d_p=4A+d_t$ and $d_t$ are the period of the variable-diameter capillary, the pore diameter, and the throat diameter, respectively. The larger the value of $T$, the less the numbers of the pore and throat, as shown in Figure 10. We can obtain capillaries of different shapes by adjusting the size of parameters ($A$, $T$ and $d_t$) in Equation (19).

Because the wetting phase flows along the solid wall under the action of capillary force, in the pore-throat model, the actual distance of fluid movement is approximately equal to the length of the inner wall of the capillary ($L_S$). The ratio of the inner wall length to the horizontal length ($L_s/L_o$) is defined as pore-throat tortuosity (${\tau }_{pt}$), which is different from the usual tortuosity ($\mathsf{\tau }=L_t/L_o$), as shown in Figure 11. The larger the pore-throat tortuosity, the greater the fluctuation in the capillary diameter, and its minimum value is equal to 1. The pore-throat tortuosity is controlled by $A$ and $T$ in Equation (19).

3.2. The Simulated Results and Analysis

3.2.1. Results and Analysis of Spontaneous Imbibition with Snap-Off

In simulations, we found that the snap-off of the wetting phase fluid sometimes occurred during the spontaneous imbibition process. Most of the non-wetting phase fluid is trapped in the capillary after the snap-off, resulting in very low recovery of the non-wetting phase. The snap-off refers to the phenomenon in which the non-wetting phase is wrapped by the wetting phase in the macropore during spontaneous imbibition, as shown in Figure 12a and Figure 13. In Figure 12a and Figure 13, the position of the red box shows the snap-off; and in Figure 12b,c, the wetting phase fluid did not cause the snap-off, but the red circle in Figure 12b shows a tendency to snap-off.

To study the effects of capillary structures on the snap-off, we plotted the distribution of snap-off and non-snap-off phenomena to characterize the effects of the pore-throat aspect ratio and the pore-throat tortuosity of the capillary on the snap-off, as shown in Figure 12. The pore-throat aspect ratio and pore-throat tortuosity comprehensively reflect the capillary structure characteristics.

We can see that there is a clear black boundary in Figure 14. Above the black dashed line, the snap-off occurs during spontaneous imbibition, whereas, below the dashed line, the snap-off does not occur. When the pore-throat tortuosity is constant, the larger pore-throat aspect ratio, the more likely that snap-off occurs. In addition, when the pore-throat aspect ratio of the variable-diameter capillary is constant, the pore-throat tortuosity is larger, and the snap-off is more likely to happen. Overall, the pore-throat aspect ratio and pore-throat tortuosity are highly relevant to the occurrence of snap-off.

In addition, to further study the influence of diameter changes on the snap-off, we compared the spontaneous imbibition of three types of capillaries, as shown in Figure 15. The mean diameters of the three capillaries are the same, and the pore-throat tortuosity of the two curved capillaries are the same. It can be seen from Figure 15 and Figure 16 that spontaneous imbibition does not occur in the invariable-diameter capillary, but the phenomenon is completely different in the variable-diameter capillary.

3.2.2. Results and Analysis of Spontaneous Imbibition without Snap-Off

Without snap-off, when we fix the pore-throat aspect ratio, the simulated spontaneous imbibition velocities vs. throat diameters for different periods are shown in Figure 17. According to the figure, spontaneous imbibition velocity increases when the throat diameter increases with the pore-throat aspect ratio fixed. In addition, when the period increases, i.e., the diameter changing rate decreases, the spontaneous imbibition velocity also increases.

Furthermore, the fluid interface vs. time step is shown in Figure 18 for different pore diameters, with a fixed throat diameter $d_t=$ 20 and period $T=$ 40π. According to the figure, the interface moving velocity is changing during the spontaneous imbibition process, which demonstrates the fluid passing through pores and throats. The interface position in the black rectangle is around the maximum diameter, i.e., the pore region. More importantly, the figure suggests that the spontaneous imbibition velocity increases when the pore diameter decreases (the pore-throat aspect ratio also decreases) with a fixed throat diameter.

Regarding different periods, when the period ($T$) is 40 π, 50 π, and 70 π, and the throat ($d_t$) is fixed, the average spontaneous imbibition velocities for different pore diameters ($d_p$) are shown in Figure 19. It can be seen from Figure 19 that when the period and the throat are fixed, the average spontaneous imbibition velocity decreases with the increase in pore diameter, and the decreasing rate is slower and slower. Different periods lead to different exponents of the fitting equation. Generally, when the period increases, i.e., the diameter changing rate decreases, the exponent of the trendline decreases. Overall, when the throat diameter is fixed, a bigger pore diameter inhibits the speed of spontaneous imbibition.

To further investigate the effects of the diameter changing rate on the spontaneous imbibition process, we set up two contrasting scenarios. In Figure 20a, two capillaries are noted by ① and ②, where the pore diameter, throat diameter, and the average diameter of the capillary are the same, with values of 60, 20, and 40, respectively. Moreover, the filled fluid properties are completely the same. The two simulated curves of ① and ② in Figure 20b correspond to the spontaneous imbibition process of ① and ② capillaries in Figure 20a. It can be seen from the curve that the flat section in the rectangles of capillary ① is much longer than that of capillary ②. Although the pore, throat size, and average diameter of the two models are the same, the changing rate of their straight diameter is different (the model ① diameter changes suddenly, resulting in large additional resistance; the model ② diameter changes smoothly). In summary, when the throat diameter is fixed, the larger the pore diameter (the larger the pore-throat aspect ratio), the smaller the spontaneous imbibition force (the smaller the capillary force), and the larger the additional resistance, which results in a smaller spontaneous imbibition velocity.

4. Conclusions

We applied the color gradient model and the parallel computation technique using CPUs to study the spontaneous imbibition of variable-diameter capillaries. The model was first validated by comparing the simulated results with the classical LW equations for different viscosities, contact angles, and interfacial tensions. The following conclusions are drawn:

• The existence of the variable diameter makes the snap-off occur, which reduces the recovery of the non-wetting phase.
• The pore-throat aspect ratio and the pore-throat tortuosity are the main factors affecting the snap-off. The larger the pore-throat aspect ratio, the more easily snap-off occurs. When the pore-throat aspect is less than 3.7, regardless of how the pore-throat changes, it will not cause the snap-off.
• When the snap-off does not occur, the interface moving velocity changes during the spontaneous imbibition processes and the minimum velocity corresponds to the largest pore diameter.
• The spontaneous imbibition velocity increases when the throat diameter increases and the pore-throat aspect ratio is fixed; when the period increases, i.e., the diameter changing rate decreases, the spontaneous imbibition velocity also increases.
• When the capillary throat diameter is fixed, bigger pore diameter and smaller period of the sine function both inhibit the speed of spontaneous imbibition.

The limitation of the variable-diameter capillary model is that the expansibility of solids and chemical reactions between media are ignored.