Dynamic Pore-Scale Network Modeling of Spontaneous Water Imbibition in Shale and Tight Reservoirs

Abstract

Spontaneous water imbibition plays an imperative role in the development of shale or tight oil reservoirs. Spontaneous water imbibition is helpful in the extraction of crude oil from the matrix, although it decreases the relative permeability of the hydrocarbon phase dramatically. The dynamic pore-scale network modeling of water imbibition in shale and tight reservoirs is presented in this work; pore network generation, local capillary pressure function, conductance calculation and boundary conditions for imbibition are all presented in detail in this paper. The pore network is generated based on the characteristics of Barnett shale formations, and the corresponding laboratory imbibition experiments are matched using this established dynamic pore network model. The effects of the wettability, throat aspect ratio, viscosity, shape factor, micro-fractures, etc. are all investigated in this work. Attempts are made to investigate the water imbibition mechanisms from a micro-scale perspective. According to the simulated results, wettability dominates the imbibition characteristics. Besides this, the viscous effects including viscosity, initial capillary pressure and micro-fractures increase the imbibition rate, while the final recovery factor is more controlled by the capillarity effect including the cross-area shape factor, contact angle and the average pore-throat aspect ratio.

1. Introduction

The spontaneous water imbibition phenomenon is commonly encountered in the development of shale or tight reservoirs. Water imbibition is helpful for the replacement of oil from the matrix and increases the oil recovery, while it decreases the relative permeability of the hydrocarbon phase dramatically, which is detrimental to the production rate and recovery factor [1,2]. Generally, the capillarity is believed to be the controlling factor for the imbibition of water in shale and tight formations. Analogous to waterflooded carbonate reservoirs [3,4,5], taking full advantage of water imbibition to yield the highest oil recovery is a major concern for the development of shale and tight reservoirs. Moreover, surfactant additives are proposed to alter the formation wettability and improve the performance of water imbibition [6,7,8]. Large numbers of experimental studies of spontaneous imbibition have been reported in the literature for both conventional [9,10,11,12,13] and unconventional reservoirs [14,15,16]. Large numbers of analytical models based on a bundle of capillary tubes have been proposed following the early works by Lucas [17]. Viscosity ratios, the tortuosity of capillary tubes, the shapes of cross-sectional areas and the assumption of fractal porous media are involved in these models [18,19]. An attempt has been made to find a universal upscaling equation since the imbibition rate and its recovery factor is sensitive to the sample size; i.e., the characteristic length. Unfortunately, these upscaling techniques are still not applicable for some other factors including wettability, complex rock structures and a wide range of fluid viscosities, etc.

Recently, scholars have become more focused on the micro-scale studies of water imbibition [20], and pore-scale network models have therefore been proposed to simulate this process. There are two types of pore-scale network fluid flow simulation models: quasi-static and dynamic models. For quasi-static models, the fluid flow is controlled by capillarity effects and the viscous forces are neglected. In contrast, inspired by the continuum reservoir simulation method, the dynamic pore network models consider both capillarity and viscous effects. For every increase or decrease of capillary pressures, the two-phase displacement mechanisms—including the piston-like entry of the non-wetting phase, the piston-like throat filling of the wetting phase, the snap-off of the non-wetting phase and the cooperative pore filling of the wetting phase [21,22]—are considered by quasi-static models. Based on the calculated capillary pressure and relative permeability from quasi-static models, the spontaneous imbibition process can then be simulated using conventional reservoir simulators or some simple analytical solutions [23,24,25,26]. However, we cannot simulate the process of imbibition directly from the quasi-static model, as not all calculated macro-properties have a relationship with time. Meanwhile, the quasi-static models assume capillary-dominated flow, and the capillary number is assumed to be less than 10⁻⁵ [27,28].

For the dynamic pore network models, Sheng and Thompson [29] reviewed and divided the dynamic models into three categories: semi-dynamic models (perturbative models), Washburn equation-based dynamic models and coupled dynamic models. The semi-dynamic models separately calculate each phase pressure and update the saturation profiles by ranking the pore elements using a defined global potential which is a combination of the local capillary pressure, viscous pressure drop and gravity. This idea was first derived by Blunt and Scher [30] and then further developed by Hughes and Blunt [31], Idowu and Blunt [32] and Aghaei and Piri [33]. The second type of dynamic model uses the Washburn equation to combine the viscous and capillary forces. Aker et al. [34] did pioneering work and established a framework in this area. Accounting for the wetting film flow, this model is further improved [35,36]. The third type of dynamic model is based on the mass conservation of multiphase flow in every pore element, which is similar to the conventional macroscale reservoir simulation. The Implicit Pressure and Explicit Saturation (IMPES) algorithm is applied to solve the dynamic model. As noted by Khayrat [37], there are two key difficulties resulting in the instability of numerical solutions: (1) local rules in dynamic network models cause a strong non-linearity in the non-wetting conductance because of the occurrence of non-continuity; (2) the capillary pressure is fully coupled within the displacement process, which is a sensitive function of phase saturation. Therefore, an appropriate algorithm is needed to ensure numerical stability, and several related studies have been conducted [38,39].

In this work, we propose a dynamic pore-scale network fluid flow simulation model that simulates the process of spontaneous water imbibition in shale or tight formations. The newly derived mathematical model and corresponding solution algorithm are presented. The proposed model is first applied and validated for Barnett shale formations for the single-phase flow, and then the spontaneous imbibition process is simulated and corresponding systematical sensitivity studies are conducted. Attempts are made to propose an approach to investigate the spontaneous water imbibition mechanisms from a micro-scale perspective.

2. Dynamic Pore Network Modelling of Water Imbibition

2.1. Pore Network Generation

Capillary tube models are suitable for single-phase flow; however, they fail to simulate the multiphase flow scenarios in a porous medium. Pore-scale network models bridge the gap between the single capillary tube model and the porous media system. Firstly, we need to generate the pore network for the studied porous media before simulating two-phase flow processes within them. A structured network consisting of nodes and connecting bonds is used, where nodes are assigned as pore bodies and bonds as pore-throats. A truncated lognormal density distribution function (c.f. Equation (1)) commonly yields a good fit for the profiles of the pore-throat distribution of reservoir rocks [38]. Since the larger pores tend to be connected with larger throats, we cannot randomly assign the pore bodies connected to the pore-throats; here, we sort the pore bodies by the mean values of the connected throats’ radii, and then the corresponding sorted pore body radii are assigned.

$$f\left(r,\sigma \right)=\frac{\sqrt{2}\exp \left[-\frac{1}{2}{\left(\ln \frac{r}{r_m}\right)}^2\right]}{\sqrt{\pi {\sigma }^2}{r}_i\left[\mathrm{erf}\left(\frac{\mathit{ln}\frac{r_{max}}{r_m}}{\sqrt{2{\sigma }^2}}\right)-\mathrm{erf}\left(\frac{\mathit{ln}\frac{r_{min}}{r_m}}{\sqrt{2{\sigma }^2}}\right)\right]}$$

(1)

where, r refers to the radius of pores or throats, and $r_m,r_{max},r_{min}$ refer to the mean, maximum and minimum radii of pores or throats, respectively; $\sigma$ is the standard derivation of the radius distribution.

Following the work by Mason and Morrow [40], in which they studied the capillarity behavior of drainage and imbibition in irregular perfect wetting triangular micro-tubes, the definition of the shape factor is used as shown in Equation (2).

$$G=\frac{A}{P^2}$$

(2)

where $A$ and $P$ are the cross-sectional area and the perimeter of a pore or throat.

For a cross-sectional shape, if we define the corner half angles of the triangle as ${\beta }_1<{\beta }_2<{\beta }_3$, two constraints need to be met, given a shape factor value:

$${\beta }_1+{\beta }_2+{\beta }_3=\frac{\pi }{2}$$

(3)

$$G=\frac{A}{P^2}=\frac{1}{4{{\displaystyle \sum }}_{i=1}^3\cot {\beta }_i}=\frac{1}{4}\tan {\beta }_1\tan {\beta }_2\cot \left({\beta }_1+{\beta }_2\right)$$

(4)

where ${\beta }_2$ is randomly chosen to range from the lower and upper boundaries in Equations (5) and (6).

$${\beta }_2^{max}=\mathrm{atan}\left\{\frac{2}{\sqrt{3}}\cos \left[\frac{\mathrm{acos}\left(-12\sqrt{3}G\right)}{3}+\frac{4\pi }{3}\right]\right\}$$

(5)

$${\beta }_2^{min}=\mathrm{atan}\left\{\frac{2}{\sqrt{3}}\cos \left[\frac{\mathrm{acos}\left(-12\sqrt{3}G\right)}{3}\right]\right\}$$

(6)

The mean values of Equations (5) and (6) are used; thus, the smallest half corner angle can be calculated by

$${\beta }_1=-\frac{1}{2}{\beta }_2+\frac{1}{2}\mathrm{asin}(\frac{\tan {\beta }_2+4G}{\tan {\beta }_2-4G}\sin {\beta }_2)$$

(7)

Finally, ${\beta }_3$ is calculated by applying Equation (3). In this way, we can use a single shape factor for the whole pore network. Therefore, one single shape factor will help us to use a single dimensionless local capillary function for every pore or throat with variable wetting phase saturation, which will be presented later. The schematic of the pore network and the triangular cross-section of every pore element are presented in Figure 1.

2.2. Control Equations and Conductance Calculation

The proposed dynamic pore-scale network model is based on the following assumptions:

(1) The volume of throats is negligible compared with that of pore bodies, while the throats control the conductances between pore bodies; (2) the pressure field is continuous, and an upwind scheme is used for the numerical treatment, while the saturation could be discontinuous, especially for the connecting throats; (3) fluids are assumed to be incompressible, and the pore network is rigid.

Applying the mass conservation for every pore body, the governing equations are obtained as

$$V_i\frac{\partial S_w}{\partial t}={\displaystyle \sum }_j\frac{K_{ij}{K}_{ij}^{rw}}{{\mu }_w}\left(p_j^w-p_i^w\right)$$

(8a)

$$V_i\frac{\partial S_n}{\partial t}={\displaystyle \sum }_j\frac{K_{ij}{K}_{ij}^{rn}}{{\mu }_n}\left(p_j^n-p_i^n\right)$$

(8a)

where $V, S, \mathrm{and} p$ refer to the volume, phase saturation and pressure for every pore body, respectively, and $K_{ij}$ and $K_{ij}^r$ are the absolute and relative conductance between the $i^{th}$ and $j^{th}$ pores.

The geometry of pores or throats controls the conductance of multiphase flow. For single-phase flow in an irregular triangle, the absolute conductance is obtained by [28]

$$K=0.6\frac{G{A}^2}{L}$$

(9)

where $L$ is the length of a pore or throat.

For the two-phase flow scenario of the triangular cross-section, the non-wetting phase occupies the center, while the wetting phase resides in corners, which leads to extra resistances for wetting phase flow [41]. Following the study by Sheng and Thompson [29], the relative conductances for two-phase flow in a throat during the imbibition process are further derived in Equation (10a,b):

$$K^{rn}=\{\begin{array}{c}{S}_n^2 , S_w<{\left[F\left(G,\theta \right)\frac{\sigma }{r}\frac{1}{p_c}\right]}^2\\ 0, S_w\ge {\left[F\left(G,\theta \right)\frac{\sigma }{r}\frac{1}{p_c}\right]}^2\end{array}$$

(10a)

$$K^{rw}=\{\begin{array}{c}C\left(G,\theta \right)S_w^2 , S_w<{\left[F\left(G,\theta \right)\frac{\sigma }{r}\frac{1}{p_c}\right]}^2\\ 1, S_w\ge {\left[F\left(G,\theta \right)\frac{\sigma }{r}\frac{1}{p_c}\right]}^2\end{array}$$

(10b)

where $C\left(G,\theta \right)$ is the function of the shape factor and contact angle, accounting for the extra resistance in corners of the wetting phase, and $S_w$ and $S_n$ refer to the saturations of wetting and non-wetting phases, respectively; $F\left(G,\theta \right)$ will be introduced later in Equation (13).

Following the analytical derivation of Øren et al. [21] and Valvatne and Blunt [22], the expression of $C\left(G,\theta \right)$ is rearranged in Equation (11a–c). Note that $C\left(G,\theta \right)$ does not depend on the local capillary pressure and saturation, which is only related to the cross-sectional shape factor and contact angle. Given a shape factor and a contact angle, a constant value of $C\left(G,\theta \right)$ is calculated and will be used in Equation (10b).

$$C\left(G,\theta \right)=\frac{{{\displaystyle \sum }}_{i=1}^3\frac{\left(0.364G_{ci}+0.28G_i^{*}\right)}{0.6G}{S}_{wi}^2}{S_w{}^2}=\frac{{{\displaystyle \sum }}_{i=1}^3\frac{\left(0.364G_{ci}+0.28G_i^{*}\right)}{0.6G}{\left[\frac{cos\theta \cos \left(\theta +{\beta }_i\right)}{sin{\beta }_i}+\theta +{\beta }_i-\frac{\pi }{2}\right]}^2}{{\left\{{{\displaystyle \sum }}_{i=1}^3\left[\frac{cos\theta \cos \left(\theta +{\beta }_i\right)}{sin{\beta }_i}+\theta +{\beta }_i-\frac{\pi }{2}\right]\right\}}^2}$$

(11a)

$$G_i^{*}=\frac{\sin {\beta }_i\cos {\beta }_i}{4{\left(1+\sin {\beta }_i\right)}^2}$$

(11b)

$$G_{ci}=\frac{\left[\frac{cos\theta \cos \left(\theta +{\beta }_i\right)}{sin{\beta }_i}+\theta +{\beta }_i-\frac{\pi }{2}\right]}{4{\left[\frac{\cos \left(\theta +{\beta }_i\right)}{sin{\beta }_i}+\frac{\pi }{2}-\theta -{\beta }_i\right]}^2}$$

(11c)

2.3. Local Capillary Pressure Function

For every pore element (pore body or throat), the fluid distribution is controlled by the local capillary pressure function, which is dependent on the element radius, shape of the cross area and contact angle. Based on the previous studies [27,42,43], the local capillary pressure function is reformulated as shown in Equation (12):

$$p_c=B\left(G,\theta \right)\frac{\sigma }{r}\frac{1}{\sqrt{S_w}}$$

(12)

where $\mathrm{B}\left(G,\theta \right)=2\sqrt{\mathrm{G}{\displaystyle \sum }_{i=1}^3\left[\frac{cos\theta \cos \left(\theta +{\beta }_i\right)}{sin{\beta }_i}+\theta +{\beta }_i-\frac{\pi }{2}\right]}$.

However, Equation (12) is only partially applicable, as either snap-off or piston-like-filling happens before the wetting phase saturation reaches unity. The snap-off occurs when the upstream pore element is not filled by the wetting phase, while piston-like-filling occurs when the upstream element is filled by the wetting phase. The corresponding criteria [21,22] are shown in Equations (13) and (14) for snap-off and piston-like-filling, respectively.

$$F\left(G,\theta \right)=\cos \theta -\frac{2\sin \theta }{\cot {\beta }_1+\cot {\beta }_2}$$

(13)

$$\mathrm{E}\left(G,\theta \right)=\cos \theta +\sqrt{{\cos }^2\theta +4G\left(\pi -\frac{2}{3}\theta +3\sin \theta \cos \theta -\frac{{\cos }^2\theta }{4G}\right)}$$

(14)

If using the dimensionless capillary pressure as $\frac{p_c}{\sigma /r}$, the local capillary function is shown in Figure 2. Before snap-off or piston-like-filling occurs, the local capillary pressure function follows Equation (12). After the pore filling occurs, the local capillary function becomes a horizontal line (Equation (13) or (14)) until the wetting phase saturation reaches unity. A smooth technique is used in this work to make the local capillary function monotonic. The general capillary fitting correlation by Andersen et al. [44] is used. Thus, given either saturation or capillary pressure, the other situation can be reached.

2.4. Solution Scheme and Boundary Conditions

The solution scheme follows the IMPES algorithm. Firstly, Equation (8a,b) are combined and discretized to yield the implicit pressure equation as

$${\displaystyle \sum }_j\left\{{\left[\frac{K_{ij}{K}_{ij}^{rw}}{{\mu }_w}\right]}^l\left(p_j^{w,l+1}-p_i^{w,l+1}\right)+{\left[\frac{K_{ij}{K}_{ij}^{rn}}{{\mu }_n}\right]}^l\left(p_j^{n,l+1}-p_i^{n,l+1}\right)\right\}=0$$

(15)

where the local capillary pressure is treated in an explicit form (using the value from current time step $l$); i.e., $p^{n,l+1}=p^{w,l+1}-p_c\left(S_w^l\right)$.

After solving the pressure equations, the next step is to update the saturation explicitly using Equation (8a) a in discretized form:

$$S_w^{l+1}=S_w^l+\frac{\Delta t}{V_i}{\displaystyle \sum }_j\left\{{\left[\frac{K_{ij}{K}_{ij}^{rw}}{{\mu }_w}\right]}^l\left(p_j^{w,l+1}-p_i^{w,l+1}\right)\right\}$$

(16)

Setting the maximum time step size for every time step interval is necessary to ensure that the solution scheme is stable. Here, we only allow a single pore body to be filled by the wetting phase for every step. The maximum time step size estimation from the current step $l$cannot ensure the stability of the calculation of the next step $l+1$due to the severe nonlinearity of this problem. This instability issue can be prevented by rechecking the maximum water saturation changing at the $l+1$ step. If this is less than unity, the calculation proceeds to the next step $l+1$; otherwise, the water saturation at the $l+1$ step is recalculated using a smaller time step size. The details can be found in the work by Aziz and Settari [45]. After obtaining the saturation field, the capillary pressure is updated; then, the calculation proceeds until the preset maximum time. Moreover, the throats are treated explicitly by applying the Equations (8) and (15), which significantly reduces the nonlinearity of the conductivity of non-wetting phase.

The boundary conditions for spontaneous imbibition are set as follows:

$$p^w\left(@\mathrm{Inlet}\right)=0$$

(17a)

$$p_c\left(@\mathrm{Inlet}\right)=0$$

(17b)

$$p_c\left(@\mathrm{Outlet}\right)=p_c^{max}$$

(17c)

$$p^w\left(@\mathrm{Outlet}\right)=-p_c^{max}$$

(17d)

where the initially capillary pressure of every pore element is set as the maximum pressure $p_c^{max}$.

3. Dynamic Pore-Scale Network Modeling of Water Imbibition in Shale and Tight Formations

Firstly, the pore network of Barnett shale is generated using the truncated lognormal distribution of Equation (1), and the single-phase flow simulation is validated. Moghaddam and Jamiolahmady [46] presented the experimental study of the apparent gas permeability of Barnett shale samples; the corresponding pore-throat distributions were obtained from high-pressure mercury injection capillary pressure (HPMICP) experiments. We use Equation (1) to fit the experimental throat radius distribution data of HPMICP, and the mean, minimum and maximum of the throat radii are 6.2, 1 and 50 nanometers, respectively. Since the pore throat aspect ratio is non-trivial to be obtained for shale, the value is set as 1.7 following the works of Valvatne and Blunt [22]. For gas flow in shale and tight formations, the smaller the pore size, the more important the effect of a gas–solid collision; therefore, the corresponding non-Darcy effect needs to be considered. After considering the gas non-Darcy flow (the details are provided in our previous work [47]) for every pore element, the calculated apparent gas permeability and experimentally measured data show a good match, as shown in Figure 3, which demonstrates that the generated network is representative for Barnett Shale, which forms the foundation for the following two-phase flow simulation.

Using the obtained pore network for Barnett shale, the initialization is conducted by setting every pore and throat element to have the local capillary pressure $p_c^{max}=10 \mathrm{MPa}$. Then, boundary conditions (Equation (17a) through (17d)) are applied in our dynamic pore network model, and the upper and lower boundaries are set as periodic. The water imbibition process can then be simulated directly within this proposed dynamic pore network model representing Barnett shale. Figure 4 shows the effect of the network size on the simulated spontaneous imbibition recovery with respect to time. We can see that a size of 20 × 60 is sufficient to capture the accuracy and efficiency for the following studies. The water saturation profiles of the water imbibition process are shown in Figure 5, where the average water saturation equals 0.2239. Note that the contact angle is set as 60° as the advancing contact angle is larger than the intrinsic contact angle due to the hysteresis effect [48]. All the mentioned conditions are set as the base case.

Generally, a mixed wettability is assumed for shale and tight formations. Here, we considered mixed wet conditions for Barnett shale. We assigned three types of mixed wet conditions—30% oil-wet, 50% oil-wet and 70% oil-wet—as shown in Figure 6. The corresponding spontaneous water imbibition-induced oil recoveries are shown in Figure 7 for dimensionless time $t_D=Ct\sqrt{\frac{k}{\varphi }}\frac{\sigma }{{\mu }_w}\frac{1}{L^2}$, where $C$ is the unit conversion factor, which is equal to 0.018849 if $t$ is in minutes, $k$ in md, $\varphi$ in decimal, $\sigma$ in $mN/m$, ${\mu }_w$ in cp, and L is a characteristic length in cm [49]; we follow the same units here. According to the figure, with the increasing proportion of oil-wet pores, the imbibition rate and recovery factor decrease dramatically. The results indicate that wettability dominates the water imbibition characteristics.

The existence of micro-fractures is considered in this work to investigate their effect on water imbibition, as natural and hydraulic fractures are presented within shale and tight formations. Four lines of micro-fractures are assigned within the Barnett shale pore network, and the pore radius is set as 20 nanometers (c.f. Figure 8), which is 3.23 times the mean radius of the Barnett shale matrix; therefore, the micro-fracture permeability is approximately 10 times that of the matrix permeability. The corresponding imbibition recovery with respect to dimensionless time is shown in Figure 9. With the existence of micro-fractures, the imbibition rate increases significantly while the final recovery remains similar. The micro-factures contribute to the conductance of the pore network instead of the pore volume, and the increase of conductance increases the imbibition rate but not the final recovery, as the final recovery is controlled by the capillarity of the matrix pores and throats.

Based on the above analysis, the simulated imbibition recovery factor with respect to dimensionless time is compared with the laboratory experimental work of Barnett shale samples by Morsy and Sheng [16]. The corresponding results are shown in Figure 10. Both the effects of wettability and micro-fractures are considered to fit the experimental data. Mixed wettability for Barnett shale is found within the fitting process. Since micro-fractures are observed in the imbibition experiments, the characteristic length is reduced significantly. In the fitting results, the effective characteristic length of Barnett shale samples is approximately several millimeters, and the fraction of oil-wet pores ranges from 70% to 50%.

Moreover, the effects of the initial maximum capillary pressure and the viscosity of the non-wetting phase on water imbibition are studied in this work, and the corresponding results are shown in Figure 11 and Figure 12. According to the figures, when the initial capillary pressure increases and the non-wetting phase viscosity decreases, the corresponding imbibition rates increase. However, the final imbibition recovery factor tends to be similar, which is because the water imbibition process is a capillary-dominated flow and the recovery factor is controlled by capillarity, which is affected by pore shape, aspect ratio, contact angle, etc. Therefore, the viscous forces in the reservoir conditions increase the imbibition rates but not the final recovery factor.

To further investigate the capillarity effect on the water imbibition of Barnett shale, the following sensitivity studies are conducted: the determination of the pore throat aspect ratio, contact angle and shape factor of the cross-area. The corresponding results for the water imbibition recovery factor with respect to dimensionless time are shown in Figure 13, Figure 14 and Figure 15, respectively. A higher aspect ratio tends to increase the percentage of the non-wetting phase trapped by the snap-off effect, which increases the residual non-wetting phase saturation and hence the final recovery factor (c.f. Figure 13). When the contact angle and cross-area shape factor decrease, leading to an increase of capillary pressure, the imbibition rate increases at the beginning. However, the final recovery factor decreases slightly (c.f. Figure 13 and Figure 14) because snap-off tends to occur more frequently with a lower contact angle and shape factor, as shown in Equation (13).

4. Summary and Conclusions

In this paper, a dynamic pore network model is proposed, including the pore network generation, local capillary pressure function, conductance calculation and boundary conditions for water imbibition. The proposed model is applied to Barnett shale formations, and the corresponding laboratory imbibition experiments are matched using the dynamic established pore network model. Then, systematical sensitivity studies are conducted. According to the simulated results using our dynamic pore network model, wettability dominates the water imbibition characteristics. Besides this, the viscous effects including viscosity, initial capillary pressure and micro-fractures increase the imbibition rate, while the final recovery factor is controlled to a greater extent by the capillarity effect including the cross-area shape factor, contact angle and aspect ratio.