Lattice Boltzmann modeling of the coupled imbibition-flowback behavior in a 3D shale pore structure under reservoir condition

Imbibition and flowback of fracturing fluid usually occur in the shale matrix after hydraulic fracturing, which significantly impacts shale gas production and environmental protection. The rocks of deep shale gas reservoirs are under high-temperature and high-temperature conditions. There are rich micro-nano pores with various pore structures in deep shale. In addition, the flowback behavior is significantly affected by the imbibition behavior because the flowback begins after the end of the imbibition. Therefore, an accurate pore-scale description of the coupled imbibition-flowback behavior is crucial to understand the flowback mechanism and its impacts. In this paper, a pseudo-potential lattice Boltzmann method is employed to simulate the coupled imbibition-flowback behavior in a digital shale core, where the digital core is reconstructed by Markov Chain-Monte Carlo method based on scanning microscope images of deep shale cores. The microcosmic mechanism of the imbibition and flowback is studied under deep shale conditions. The influence of some factors, such as pore structure, fluid viscosity, wettability, and flowback pressure difference, on the flowback behavior of fracturing fluid is investigated. It is found that the fracturing fluid advances almost uniformly throughout the pore space during the imbibition process. The fracturing fluid is easy to adsorb on the pore wall, and the shale gas is located in the middle of the pore space. The viscous fingering is clearly observed during the flowback process, where shale gas flows through large pores to form a flow channel, and the fracturing fluid stays in tiny pores. The flowback rate increases gradually with the flowback time and eventually tends to be almost constant. The wettability, flowback pressure difference, and pore structure significantly influence the flowback behavior, while the fracturing fluid viscosity has a smaller effect on the flowback process.

1 Introduction

As an unconventional gas resource with massive reserves, shale gas has received wide attention since it entered the commercial development stage (Zou et al., 2015a; Zou et al., 2016). According to the US Energy Information Administration (EIA) forecast, shale gas reserves in China account for about 20% of the world, ranking first globally (Ma et al., 2017), and thus shale gas has a broad development prospect in China (Dong et al., 2015; Hao et al., 2017; Sun et al., 2020). In this decade, shale gas production in China has grown steadily every year, but it mainly comes from middle and shallow shale gas reservoirs with a depth of smaller than 3,500 m. Deep shale gas reservoirs with a depth of over 3,500 m have received particular attention and have been regarded as the new growth point for China’s shale gas production in the future (Ma et al., 2021). Deep shale gas reservoirs are difficult to exploit because of the deep shale’s low porosity and ultra-low permeability. Horizontal well drilling technology (Nie et al., 2013) and hydraulic fracturing technology (Rongved et al., 2019) are the key technologies for the commercial development of shale gas reservoirs. Fracturing technology has been continuously improved and developed in recent years (Wang et al., 2015; Li et al., 2018; Cheng et al., 2019; Wu et al., 2022; Yang et al., 2022). After hydraulic fracturing, the flowback stratagem of fracturing fluid is one of the critical factors for shale gas production (Wang et al., 2020a; Qu et al., 2021).

The flowback rate of fracturing fluid is generally low and varies significantly for different shale gas reservoirs (Lutz et al., 2013; Vengosh et al., 2014; Gao et al., 2017). Compared with conventional gas reservoirs, deep shale gas reservoirs under high-temperature and high-temperature conditions have rich micro-nano pores with various pore structures (Guo et al., 2014; Chen et al., 2019), which affects the gas-liquid flow law and the flowback characteristics of gas wells. Because the flowback characteristics of gas wells are the macroscopic manifestation of the microcosmic mechanism of the imbibition and flowback, it is crucial to clarify the microcosmic mechanism of imbibition and flowback of fracturing fluid under reservoir conditions for the efficient development of deep shale gas reservoirs (Zheng et al., 2021a; Qu et al., 2021; Shao et al., 2022; Xu et al., 2022).

Lattice Boltzmann method (LBM) is considered as a powerful tool for simulating many complex flows such as multiphase flows (Wang et al., 2023), multicomponent flows (Ren et al., 2019), suspended particle flows (Jiang et al., 2022), blood flows (Cherkaoui et al., 2022), liquid slip flows (Ren et al., 2023), turbulent flows (Kareem and Asker, 2022), nanofluid heat flows (Khan et al., 2022), and so on. Recently LBM has been widely used to simulate pore-scale flows of multiphase fluid in porous media because it is easy to deal with complex geometric boundaries and the interactions between multiphase fluid and fluid and between fluid and solid (Ansumali and Karlin, 2005; Chikatamarla and Karlin, 2009; Zhang et al., 2009). At present, LBM models used for multiphase flow simulations mainly include pseudo-potential model (Shan and Chen, 1993), color gradient model (Gunstensen et al., 1991), free energy model (Swift et al., 1995) and phase field model (Shao et al., 2013). Combined with digital core technology (Zhu et al., 2012), LBM is widely used to study the seepage mechanism of multiphase fluid in porous media (Cantisano et al., 2013). Furthermore, compared with LBM, the molecular dynamics method is also a powerful tool to study the microcosmic mechanisms (Huang et al., 2020; Huang et al., 2021a; Huang et al., 2021b). But molecular dynamics is subjected to expensive computational costs, and so it is usually applied to multiphase flows in a simple nanoscale pore structure rather than in porous media with complex pore structure. In recent years, with the wide application of hydraulic fracturing technology in the development of unconventional oil and gas reservoirs (Zou et al., 2015b; Zhang et al., 2018; Li et al., 2022), some scholars have begun to use LBM to study the microscopic imbibition and flowback mechanism of fracturing fluid. In 2015, Zhang et al., 2015 used LBM to simulate the process of gas flooding in saturated water cores based on digital shale cores with a porosity of 16.5% and confirmed that most of the water would remain in the pores and could not be discharged. In 2018, Zheng et al., 2018 used LBM to simulate the influence of wettability heterogeneity distribution on the propagation and distribution of fracturing fluid based on a digital shale core with a porosity of 26.98%. In 2020, Wang et al., 2020b used LBM to simulate the process of water imbibition and oil displacement under the action of capillary force based on the tight sandstone with porosity ranging from 10.3% to 18.6% and studied the influence of pore structure and wettability on spontaneous imbibition. In 2021, Zheng et al., 2021b used LBM to simulate the imbibition and displacement of heavy oil under reservoir conditions based on the digital core of tight sandstone with a porosity of 9.47% and analyzed the distribution law of the oil-water phase. Zhang et al., 2021 developed a pseudopotential-based lattice Boltzmann (LB) method to simulate gas/water two-phase flow at pore scale and simulate gas/water two-phase flow in dual-wettability nanoporous media.

Although some scholars use LBM to simulate the imbibition (or flowback) process of fracturing fluid, they only simulate a single imbibition (or flowback) process and do not study the coupled imbibition-flowback process. The imbibition process affects the flowback behavior because the gas-liquid saturation distribution at the beginning of flowback is determined by the gas-liquid saturation distribution after imbibition. It is necessary to simulate the coupled imbibition-flowback process to thoroughly investigate the flowback mechanism of fracturing fluid. In addition, the porosity of the digital core used to simulate the imbibition (or flowback) process of fracturing fluid is relatively large, resulting in the connectivity of the pore structure being good, which is inconsistent with the actual pore structure of deep shale.

In this work, we reconstruct a digital core of deep shale, use LBM to simulate the coupled imbibition-flowback process of fracturing fluid, study the microcosmic mechanism of the imbibition and flowback, and analyze the influence law of critical parameters on the flowback of fracturing fluid under reservoir conditions.

2 Lattice Boltzmann model

The three-dimensional nineteen-discrete-velocity (D3Q19) pseudo-potential model with multiple relaxation times (MRT) (Shan and Chen, 1993; Zheng et al., 2018) is carried out to simulate the imbibition and flowback processes of fracturing fluid. The MRT model has much better numerical stability than the Bhatnagar-Gross-Krook (BGK) model in the simulation of the flows involving complex boundaries such as porous media flow. In the pseudo-potential model, $k$ particle distribution functions are used to represent $k$ component fluids. The evolution equation of the distribution function of the $k$ component is:

$f_{\alpha }^k\left(\mathit{x}+{\mathit{e}}_{\alpha }\Delta t,t+\Delta t\right)-f_{\alpha }^k\left(\mathit{x},t\right)={\mathrm{\Omega }}_{coll}^k+{\mathrm{\Omega }}_{forces}^k(1)$

where $f_{\alpha }^k\left(\mathit{x},t\right)$ is the equilibrium distribution function of the $k$ component in ${\mathit{e}}_{\alpha }$ direction, $\mathit{x}$ is the spatial coordinate, ${\mathit{e}}_{\alpha }$ is the lattice speed in $\alpha$ direction, $\Delta t$ is the time step, $t$ is the time, ${\mathrm{\Omega }}_{coll}^k$ is the collision operator of the $k$ component, and ${\mathrm{\Omega }}_{forces}^k$ is the external force operator of the $k$ component.

The collision operator ${\mathrm{\Omega }}_{coll}^k$ of the $k$ component can be expressed as:

${\mathrm{\Omega }}_{coll}^k=M^{-1}{\mathrm{\Lambda }}^{k}M\left[f_{\alpha }^{eq,k}\left(\mathit{x},t\right)-f_{\alpha }^k\left(\mathit{x},t\right)\right](2)$

where $M$ is the conversion matrix; ${\mathrm{\Lambda }}^k$ is a diagonal matrix about relaxation times, expressed as ${\mathrm{\Lambda }}^k={\left[s_c^k,s_e^k,s_{\epsilon }^k,s_c^k,s_q^k,s_c^k,s_q^k,s_c^k,s_q^k,s_v^k,s_{\pi }^k,s_v^k,s_v^k,s_v^k,s_v^k,s_v^k,s_m^k,s_m^k,s_m^k\right]}^T$, $f_{\alpha }^{eq,k}\left(\mathit{x},t\right)$ is the equilibrium distribution function, which is expressed as:

$f_{\alpha }^{eq,k}\left(\mathit{x},t\right)=w_{\alpha }{\rho }^k\left(\mathit{x},t\right)\left[1+3\frac{{\mathit{e}}_{\alpha }\cdot {\mathbf{u}}^{eq}}{c^2}+9\frac{{\left({\mathit{e}}_{\alpha }\cdot {\mathbf{u}}^{eq}\right)}^2}{2c^2}-3\frac{{\left({\mathbf{u}}^{eq}\right)}^2}{2c^2}\right](3)$

where $w_{\alpha }=1/3\left(\alpha =0\right)$, $w_{\alpha }=1/18\left(\alpha =1-6\right)$, $w_{\alpha }=1/36\left(\alpha =7-18\right)$, ${\rho }^k$ is the density of the $k$ component, $c=\Delta x/\Delta t$ is the lattice sound velocity with $c=1$ in this study, and ${\mathbf{u}}^{eq}$ represents the equilibrium velocity.

The conversion matrix $M$ is expressed as

$M=\left[\begin{array}{ccccccccccccccccccc}1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1\\ -30& -11& -11& -11& -11& -11& -11& 8& 8& 8& 8& 8& 8& 8& 8& 8& 8& 8& 8\\ 12& -4& -4& -4& -4& -4& -4& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1\\ 0& 1& -1& 0& 0& 0& 0& 1& -1& 1& -1& 1& -1& 1& -1& 0& 0& 0& 0\\ 0& -4& 4& 0& 0& 0& 0& 1& -1& 1& -1& 1& -1& 1& -1& 0& 0& 0& 0\\ 0& 0& 0& 1& -1& 0& 0& 1& 1& -1& -1& 0& 0& 0& 0& 1& -1& 1& -1\\ 0& 0& 0& -4& 4& 0& 0& 1& 1& -1& -1& 0& 0& 0& 0& 1& -1& 1& -1\\ 0& 0& 0& 0& 0& 1& -1& 0& 0& 0& 0& 1& 1& -1& -1& 1& 1& -1& -1\\ 0& 0& 0& 0& 0& -4& 4& 0& 0& 0& 0& 1& 1& -1& -1& 1& 1& -1& -1\\ 0& 2& 2& -1& -1& -1& -1& 1& 1& 1& 1& 1& 1& 1& 1& -2& -2& -2& -2\\ 0& -4& -4& 2& 2& 2& 2& 1& 1& 1& 1& 1& 1& 1& 1& -2& -2& -2& -2\\ 0& 0& 0& 1& 1& -1& 1& 1& 1& 1& 1& -1& -1& -1& -1& 0& 0& 0& 0\\ 0& 0& 0& -2& -2& 2& 2& 1& 1& 1& 1& -1& -1& -1& -1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& -1& -1& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& -1& -1& 1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 1& -1& -1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& -1& 1& -1& -1& 1& -1& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& -1& -1& 1& 1& 0& 0& 0& 0& 1& -1& 1& -1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 1& -1& -1& -1& -1& 1& 1\end{array}\right](4)$

The external force operator ${\mathrm{\Omega }}_{forces}^k$ of the $k$ component can be expressed as:

${\mathrm{\Omega }}_{forces}^k=\frac{\Delta t}{2}\left[f_{\alpha }^{F,k}\left(\mathit{x}+{\mathit{e}}_{\alpha }\Delta t,t+\Delta t\right)-f_{\alpha }^{F,k}\left(\mathit{x},t\right)\right](5)$

where

$f_{\alpha }^{F,k}=\frac{3{\mathit{F}}_k\cdot \left({\mathit{e}}_{\alpha }-{\mathbf{u}}^{eq}\right)}{{\rho }^k{c}^2}{f}_i^{eq,k}(6)$

Let ${\overline{f}}_{\alpha }^k=f_{\alpha }^k-\frac{\Delta t}{2}{f}_{\alpha }^{F,k}$, then Eq. 1 can be rewritten as:

${\overline{f}}_{\alpha }^k\left(\mathit{x}+{\mathit{e}}_{\alpha }\Delta t,t+\Delta t\right)-{\overline{f}}_{\alpha }^k\left(\mathit{x},t\right)=M^{-1}{\mathrm{\Lambda }}^{k}M\left[f_{\alpha }^{eq,k}\left(\mathit{x},t\right)-{\overline{f}}_{\alpha }^k\left(\mathit{x},t\right)-\frac{\Delta t}{2}{f}_{\alpha }^{F,k}\right]+\Delta t{f}_{\alpha }^{F,k}(7)$

Fluid density ${\rho }^k$, fluid velocity $\mathbf{u}$, and fluid viscosity ${\mu }^k$ are defined as follows:

${\rho }^k={\displaystyle \sum }_{\alpha }{\overline{f}}_{\alpha }^k(8)$

${\rho }^k{\mathbf{u}}^k={\displaystyle \sum }_{\alpha }{\overline{f}}_{\alpha }^k{\mathit{e}}_{\alpha }+\frac{\Delta t}{2}{\mathit{F}}_k(9)$

$\mathbf{u}=\frac{{\displaystyle \sum }_k{\rho }^k{\mathbf{u}}^k}{{\displaystyle \sum }_k{\rho }^k}(10)$

${\mu }^k={\rho }^k\frac{c^2}{3}\left(\frac{1}{s_v^k}-0.5\right)\Delta t(11)$

In the pseudo-potential model, the total force ${\mathit{F}}_k$ of the $k$ component mainly includes the gas-liquid force ${\mathit{F}}_{c,k}$ and the fluid-solid force ${\mathit{F}}_{ads,k}$, which are defined as:

${\mathit{F}}_k={\mathit{F}}_{c,k}+{\mathit{F}}_{ads,k}(12)$

${\mathit{F}}_{c,k}\left(\mathit{x},t\right)=-G_c{\rho }^k\left(\mathit{x}\right){\displaystyle \sum }_i{w}_i{\rho }^{\overline{k}}\left(\mathit{x}+{\mathit{e}}_{\alpha }\Delta t,t\right){\mathit{e}}_{\alpha }(13)$

${{\mathit{F}}_{ads}}_{,k}\left(\mathit{x},t\right)=-{G_{ads}}_{,k}{\rho }^k\left(\mathit{x},t\right){\displaystyle \sum }_i{w}_{i}s\left(\mathit{x}+{\mathit{e}}_{\alpha }\Delta t\right){\mathit{e}}_{\alpha }(14)$

where $G_c$ is the gas-liquid interaction parameter, which is related to the interfacial tension; $G_{ads,k}$ is the fluid-solid interaction parameter, which is related to wettability; $s\left(\mathit{x}+{\mathit{e}}_{\alpha }\Delta t\right)$ is a marking function, if $\left(\mathit{x}+{\mathit{e}}_{\alpha }\Delta t\right)$ is a solid point, then $s\left(\mathit{x}+{\mathit{e}}_{\alpha }\Delta t\right)=1$, otherwise $s\left(\mathit{x}+{\mathit{e}}_{\alpha }\Delta t\right)=0$.

3 3D digital core reconstruction for deep shale gas reservoirs

Compared with the other reconstruction methods, Markov Chain-Monte Carlo (MCMC) method (Wu et al., 2006) has higher computational efficiency and can obtain the pore structure of porous media with good connectivity. Therefore, the MCMC method is used to reconstruct the shale digital core in this study. Deep shale cores are from Luzhou block in southern Sichuan. The scanning electron microscope (SEM) scanning image of the deep shale core is shown in Figure 1A. The SEM scanning image is denoised by the median filter algorithm (Dong and Blunt, 2009). The images before and after denoising are shown in Figure 1A and Figure 1B, respectively. It is obvious that the median filter algorithm can obtain a clearer grayscale image. Then the grayscale image is converted to a binarization image, shown in Figure 2, by the Otsu algorithm (Otsu, 1979). The conditional probabilities in three orthogonal directions are obtained by the traversing algorithm. Based on the conditional probabilities, a 3D digital core is reconstructed by the MCMC method. Figure 3 shows the reconstructed 3D digital cores of deep shale with and without isolated pores. The size of the digital core is $130\times 130\times 100$ with the resolution of 34 nm per pixel. The effective porosity of the digital core is 5.09%. Figure 4 shows the comparison between autocorrelation functions of the SEM scanning image and digital core. It is seen that the autocorrelation function of the digital core is consistent with that of the SEM scanning image, which confirms the representativeness of the reconstructed digital core (Wu et al., 2006).

FIGURE 1

FIGURE 1. SEM scanning images of deep shale core in Luzhou blockflow (A) before denoising and (B) after denoising.

FIGURE 2

FIGURE 2. Grayscale image (A) and binarization image (B) of deep shale core in Luzhou blockflow.

FIGURE 3

FIGURE 3. 3D digital cores of deep shale (A) with isolated pores and (B) without isolated pores.

FIGURE 4

FIGURE 4. Comparison between autocorrelation functions of SEM scanning image and digital core.

4 Microcosmic imbibition mechanism of fracturing fluid

According to the reservoir condition of Luzhou block in southern Sichuan, the basic parameters of LBM simulation are set as follows: the reservoir pressure $p=80\text{MPa}$, the reservoir temperature $T=130°\mathrm{C}$, and the contact angle $\theta =10°$. Based on National Institute of Standards and Technology (NIST) chemical database, some physical parameters of gas and water under $p=80\text{MPa}$; $T=130°\mathrm{C}$ are obtained as follows: the methane density ${\rho }_g=259\text{kg}/{\mathrm{m}}^3$, the water density ${\rho }_w=972\text{kg}/{\mathrm{m}}^3$, the methane viscosity ${\mu }_g=33.9\mathrm{\mu }\text{Pa}\cdot \mathrm{s}$, the water viscosity ${\mu }_w=233\mathrm{\mu }\text{Pa}\cdot \mathrm{s}$, and the gas-water interfacial tension $\sigma =40\text{mN}/\mathrm{m}$. In LBM simulation, 15 layers of pore lattices are added to the inlet and outlet of the 3D digital core to eliminate the end effect. The final size of the digital core used for simulation is $130\times 130\times 130$ lattices. The pressure boundary condition (i.e., the non-equilibrium bounce-back boundary condition) is applied to the inlet and outlet, and the bounce-back condition is applied to the solid walls.

In the simulation of the imbibition process, the pore space of the digital core is fully saturated with shale gas (i.e., methane). The fracturing fluid (i.e., water) derived by the capillary force infiltrates into the digital core and replaces the shale gas. The driving force in the imbibition process is the capillary force, which is related to the size of the pore throat. The stopping time of the simulation is 6 million time steps when the stable distribution of the two-phase fluid in the digital core can be obtained. Figure 5 shows the gas-liquid distribution in shale at different time steps and the water saturation versus time steps in the imbibition process. It is seen that the fracturing fluid advances almost uniformly throughout the pore space during the imbibition process. More liquid phase is distributed near the inlet than near the outlet. Water saturation increases with time, but the imbibition rate decreases until it finally becomes almost zero. This observation is because capillary force is the driving force in the imbibition process, and the greater the distance of fracturing fluid penetrating the core, the larger the viscous resistance of fluid migration to be overcome. Therefore, the imbibition rate gradually decreases with time, and eventually, the imbibition becomes equilibrium.

FIGURE 5

FIGURE 5. The gas-liquid distribution in shale at different time steps (A) and the water saturation versus time steps (B) in the imbibition process.

Figure 6 shows the slice diagrams of gas-liquid distribution at the equilibrium state ($t=6.0\times {10}^6$) in the imbibition process. It is seen that the fracturing fluid is easy to adsorb on the pore wall, and the shale gas is located in the middle of the pore space. Figure 7 shows the gas-liquid distribution of slice 3 in Figure 6 at different time steps. It is observed that the fracturing fluid first invaded the pore along the pore wall and then gradually occupied the central position of the pore. The larger capillary force in tiny pores results in the fracturing fluid preferentially filling smaller pores and subsequently invading larger pores.

FIGURE 6

FIGURE 6. Slice diagrams of gas-liquid distribution at the equilibrium state ($t=6.0\times {10}^6$) in the imbibition process.

FIGURE 7

FIGURE 7. Gas-liquid distribution of slice 3 in Figure 6 at different time steps.

5 Microcosmic flowback mechanism of fracturing fluid

The flowback of fracturing fluid is simulated by LBM. the gas-liquid distribution in the digital core at the initial moment of the flowback process is the same as that at the end of the imbibition process, shown in Figure 6. The shale gas is injected continuously at the right outlet, and thus the shale gas displaces the fracturing liquid from right to left. The flowback pressure difference is set as 3.0 MPa. The boundary conditions and the basic parameters are the same as those applied in the case of the imbibition process. The stopping time of the simulation is 4 million time steps when the gas-liquid distribution in the digital core almost tends to be stable.

Figure 8 shows the gas-liquid distribution in shale at different time steps and the water saturation, gas saturation, and flowback rate versus time steps in the flowback process. It is seen that the viscous fingering is clearly observed during the flowback process, where shale gas flows through large pores to form a flow channel, and the fracturing fluid stays in tiny pores, resulting in a low flowback rate of the fracturing fluid. The flowback rate increases gradually with the flowback time and eventually tends to be almost constant.

FIGURE 8

FIGURE 8. The gas-liquid distribution in shale at different time steps (A) and the water saturation, gas saturation, and flowback rate versus time steps (B) in the flowback process.

Figure 9 shows slice diagrams of gas-liquid distribution at the equilibrium state ($t=4.0\times {10}^6$) in the flowback process. The capillary force is resistance in the flowback process. Because of a smaller capillary force in a large pore, the gas phase preferentially breaks through in the middle of the macro-pores, and the fracturing fluid in small pores and blind ends is difficult to flowback.

FIGURE 9

FIGURE 9. Slice diagrams of gas-liquid distribution at the equilibrium state ($t=4.0\times {10}^6$) in the flowback process.

Figure 10 shows the gas-liquid distribution of slice 2 in Figure 9 at different time steps. It is seen that with the increase of time, the proportion of the area occupied by the gas phase in the slice gradually increases, and the increased gas phase was mainly concentrated in the large pores of the slice.

FIGURE 10

FIGURE 10. Gas-liquid distribution of slice 2 in Figure 9 at different time steps.

6 Sensitivity analysis for imbibition and flowback of fracturing fluid

In this section, based on the LBM simulations of the imbibition and flowback processes, the influence of the wettability, flowback pressure difference, pore structure, and fracturing fluid viscosity, on the imbibition and flowback behaviors are studied in detail. The simulation conditions and basic parameters are the same as those applied in imbibition and flowback cases in section 4 and 5, respectively.

6.1 Wettability of fracturing fluid

The fracturing fluid is driven by capillary force to enter the pore and to replace the gas in shale in the imbibition process, while the capillary force is the resistance in the flowback process. Because the capillary force is related to the wettability of fracturing fluid, the fluid distribution and flowback rate are usually affected by the wettability. The LBM is used to simulate the imbibition and flowback processes under different contact angles.

Figure 11 shows the gas-liquid distribution at the equilibrium state ($t=6.0\times {10}^6$) in the imbibition process under different contact angles. It is seen that the smaller the contact angle is, the stronger the wettability is, the stronger the imbibition effect of the fracturing fluid becomes, and the more the fracturing fluid intrudes into the core. When the contact angle $\theta =60°$, the imbibition effect of the fracturing fluid is significantly reduced, and the imbibition behavior mainly occurs in the pore space connected with the inlet of the digital core. Under weak wet condition, the capillary force as the driving force is small, and thus it is difficult for the fracturing fluid to further impermeable into the porous medium, and the distance of the fracturing fluid penetrating into the digital core is small. When the contact angle $\theta =10°$, the fracturing fluid has the highest imbibition efficiency and the largest penetration distance into the core, and the retention of the non-wetting phase after the imbibition front is more obvious. This is because under strong wet condition, the wetting fracturing fluid tends to fill the smaller pore structure, causing the small pores around the large pores to be filled preferentially, so that the non-wetting phase gas in the large pores is bypassed.

FIGURE 11

FIGURE 11. Gas-liquid distribution at the equilibrium state ($t=6.0\times {10}^6$) in the imbibition process under different contact angles. (A) θ = 10°, (B) θ = 30°, (C) θ = 60°.

Figure 12 shows the water saturation versus time steps in the imbibition process under different contact angles. It is observed that under different wetting conditions, the imbibition effect of the fracturing fluid gradually weakens with time, and eventually tends to the imbibition equilibrium. However, with the decrease of the contact angle, the stronger the imbibition effect of the fracturing fluid is, the higher the final water saturation of the core at the equilibrium state is.

FIGURE 12

FIGURE 12. The water saturation versus time steps in the imbibition process under different contact angles.

In the simulation of the flowback process under contact angles $\theta =10°,30°,60°$, the initial gas-liquid saturation distributions of the flowback process are the same as that at the end of the imbibition process, respectively. The flowback pressure difference is set as 1.5 MPa.

Figure 13 shows the gas-liquid distribution at the equilibrium state ($t=4.0\times {10}^6$) in the flowback process under different contact angles. As the contact angle decreases, the capillary resistance increases, and the fracturing liquid is more difficult to be displaced under the same flowback pressure difference, and thus the more fracturing fluid is retained in the pore space and the lower the flowback rate is.

FIGURE 13

FIGURE 13. Gas-liquid distribution at the equilibrium state ($t=4.0\times {10}^6$) in the flowback process under different contact angles. (A) θ = 10°, (B) θ = 30°, (C) θ = 60°.

Figure 14 shows the flowback rate versus time steps and the water saturation versus time steps in the flowback process under different contact angles. It is clear that with the increase of time, the water saturation in core gradually decreases, and the flowback rate gradually increases and finally reaches stability. With the increase of the contact angle, the flowback rate increases significantly. At the equilibrium state (i.e., $t=4.0\times {10}^6$), the flowback rates for contact angles $\theta =10°,30°,60°$ are 21.03%、26.98%、45.31%, respectively, which illustrates that the wettability has a great influence on the flowback of the fracturing fluid.

FIGURE 14

FIGURE 14. The flowback rate versus time steps (A) and the water saturation versus time steps (B) in the flowback process under different contact angles.

6.2 Flowback pressure difference

The effect of flowback pressure difference on the flowback process is investigated in this section. The gas-liquid distribution in the digital core at the initial moment of the flowback process is the same as that at the end of the imbibition process, shown in Figure 6. The flowback pressure differences are set as $\Delta p$ =0.5 MPa、1.5 MPa、3.0MPa, respectively.

Figure 15 shows the gas-liquid distribution at the equilibrium state ($t=4.0\times {10}^6$) in the flowback process under different flowback pressure differences. It is seen that the larger the flowback pressure difference is, the less the fracturing fluid is retained in the pores, and the higher the flowback rate is.

FIGURE 15

FIGURE 15. Gas-liquid distribution at the equilibrium state ($t=4.0\times {10}^6$) in the flowback process under different flowback pressure differences. (A) Δp = 0.5 MPa, (B) Δp = 1.5 MPa, (C) Δp = 3.0 MPa.

Figure 16 shows the flowback rate versus time steps and the water saturation versus time steps in the flowback process under different flowback pressure differences. It is seen that with the increase of the flowback pressure difference, the flowback rate increases significantly. At the equilibrium state (i.e., $t=4.0\times {10}^6$), the flowback rates for flowback pressure differences $\Delta p$ =0.5 MPa、1.5 MPa、3.0 MPa are 8.62%、21.03%、39.14%, respectively, which shows that the flowback pressure difference has a significant influence on the flowback of the fracturing fluid.

FIGURE 16

FIGURE 16. The flowback rate versus time steps (A) and the water saturation versus time steps (B) in the flowback process under different flowback pressure differences.

6.3 Pore structure

In order to study the effect of the pore structure on the imbibition and flowback behaviors of the fracturing fluid, another digital core with a tighter pore structure is reconstructed by the MCMC method. The digital core reconstructed in Section 3 is labeled as Core 1 and the new digital core here is labeled as Core 2 (shown in Figure 18). The permeabilities of the two cores are calculated by the LBM. The porosity and permeability of Core 1 are 5.09% and 4.0×10⁻⁴ mD, respectively. The porosity and permeability of Core 2 are 3.90% and 2.3×10⁻⁴ mD, respectively. Obviously, the porosity and permeability of Core 2 are lower than that of Core 1. The pore size distribution of the two cores is shown in Figure 17, which demonstrates that Core 2 has much tighter than Core 1.

FIGURE 17

FIGURE 17. The pore size distribution of the two digital cores.

Figure 18 shows the gas-liquid distribution at the equilibrium state ($t=6.0\times {10}^6$) in the imbibition process under different pore structures. As can be seen from Figure 18, the pore structure of Core 2 is tighter than that of Core 1, which leads to the more retention of the non-wetting gas phase after the imbibition front. This is because the wetting fracturing fluid tends to fill the smaller pore structure, causing the small pores around the large pores to be filled preferentially, so that the non-wetting phase gas in the large pores is bypassed.

FIGURE 18

FIGURE 18. Gas-liquid distribution at the equilibrium state ($t=6.0\times {10}^6$) in the imbibition process under different pore structures. (A) Core 1, (B) Core 2.

Figure 19 shows the water saturation versus time steps in the imbibition process under different pore structures. It is interesting that the tighter the pore structure is, the greater the imbibition capillary force is, the stronger the imbibition effect becomes, and the more the fracturing fluid intrudes into the core in the early stage. However, the tight pore structure leads to the increase of gas-water two-phase flow resistance, and the equilibrium time of the imbibition is shortened, and the final imbibition amount is reduced.

FIGURE 19

FIGURE 19. The water saturation versus time steps in the imbibition process under different pore structures.

Figure 20 shows the gas-liquid distribution at the equilibrium state ($t=4.0\times {10}^6$) in the flowback process under different pore structures. Since fracturing fluid is driven by gas in the flowback process, capillary force is the resistance. The tighter the pore structure is, the greater the capillary resistance is, the less the flowback amount of the fracturing fluid is, and the more the fracturing fluid is stayed far away from the inlet end.

FIGURE 20

FIGURE 20. Gas-liquid distribution at the equilibrium state ($t=4.0\times {10}^6$) in the flowback process under different pore structures. (A) Core 1, (B) Core 2.

Figure 21 shows the flowback rate versus time steps and the water saturation versus time steps in the flowback process under different pore structures. It is seen that as the pore structure becomes tighter, the flowback rate decreases. At the equilibrium state (i.e., $t=4.0\times {10}^6$), the flowback rates for Core 1 and Core 2 are 39.14% and 29.74%, respectively. Therefore, the pore structure has important effects on the flowback of the fracturing fluid.

FIGURE 21

FIGURE 21. The flowback rate versus time steps (A) and the water saturation versus time steps (B) in the flowback process under different pore structures.

6.4 Fracturing fluid viscosity

Figure 22 shows the gas-liquid distribution at the equilibrium state ($t=6.0\times {10}^6$) in the imbibition process under different fracturing fluid viscosities. In the process of imbibition, the smaller the viscosity of the fracturing fluid is, the stronger the effect of the imbibition becomes, the more amount of the fracturing fluid penetrates into the core, and the more the final imbibition amount is.

FIGURE 22

FIGURE 22. Gas-liquid distribution at the equilibrium state ($t=6.0\times {10}^6$) in the imbibition process under different fracturing fluid viscosities. (A) μw = 170 μPa.s, (B) μw = 233 μPa.s, (C) μw = 300 μPa.s.

Figure 23 shows the water saturation versus time steps in the imbibition process under different fracturing fluid viscosities. It is seen that with the decrease of fracturing fluid viscosity, the faster the imbibition speed of fracturing fluid is, and the earlier the imbibition equilibrium is reached.

FIGURE 23

FIGURE 23. The water saturation versus time steps in the imbibition process under different fracturing fluid viscosities.

Figure 24 shows the gas-liquid distribution at the equilibrium state ($t=4.0\times {10}^6$) in the flowback process under different fracturing fluid viscosities. In the flowback process, the larger the viscosity of the fracturing fluid is, the slower the flowback speed becomes, the more obvious the viscous finger advance shows, and the smaller the effect scope is.

FIGURE 24

FIGURE 24. Gas-liquid distribution at the equilibrium state ($t=4.0\times {10}^6$) in the flowback process under different fracturing fluid viscosities. (A) μw = 170 μPa.s, (B) μw = 233 μPa.s, (C) μw = 300 μPa.s.

Figure 25 shows the flowback rate versus time steps and the water saturation versus time steps in the flowback process under different fracturing fluid viscosities. It is seen that as the fracturing fluid viscosity increases, the flowback rate decreases. At the equilibrium state (i.e., $t=4.0\times {10}^6$ ), the flowback rates for fracturing fluid viscosity ${\mu }_w=$ 170 μPa s, 233 μPa s, 300 μPa s are 29.78%, 26.98%, and 22.87%, respectively. Therefore, the viscosity of the fracturing fluid has a smaller effect on the flowback of the fracturing fluid under reservoir conditions.

FIGURE 25

FIGURE 25. The flowback rate versus time steps (A) and the water saturation versus time steps (B) in the flowback process under different fracturing fluid viscosities.

7 Conclusion

In this paper, based on the SEM images of deep shale cores in Luzhou block, a 3D digital core was reconstructed by using the MCMC method. The LBM simulation of the imbibition and flowback processes in the reconstructed core was carried out, and the micro-mechanism of the imbibition and flowback of fracturing fluid was deeply studied. At the same time, the influence law of related parameters on the imbibition and flowback of the fracturing fluid was analyzed. The main results and understandings are as follows.

(1) In the imbibition process of the fracturing fluid, the distribution of the liquid phase near the inlet is more than that at the outlet, the imbibition rate gradually weakens with time, and finally the imbibition tends to be balanced. The fracturing fluid shows a uniform invasion, and the fracturing fluid is easily adsorbed on the pore wall, and the gas phase is located in the central position of the pore.

(2) In the flowback process of the fracturing fluid, the viscous fingering is obvious, the gas flows through the large pore channels, and the fracturing fluid stays in small pores, resulting in a low flowback rate of the fracturing fluid; The flowback rate gradually increases with the flowback time, and finally the flowback rate tends to be constant.

(3) The wettability, flowback pressure difference, and pore structure significantly influence the flowback behavior, while the fracturing fluid viscosity has a smaller effect on the flowback process. The pore structure becomes tighter, the flowback pressure difference becomes smaller, the contact angle becomes larger, and the viscosity of the fracturing fluid becomes larger, resulting a smaller flowback rate of the fracturing fluid.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

SW, conceptualization, methodology, software, writing—original draft preparation; JW, writing—review and editing, supervision; YL, methodology; XY, data curation; JuZ, visualization; JiZ, writing—review and editing; DZ, methodology, software; BZ, supervision; DL, writing—original draft preparation.

Conflict of interest

Authors SW, JW, YL, XY, JuZ, JiZ, DZ, BZ, and DL were employed by PetroChina Southwest Oil and Gasfield Company.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Ansumali, S., and Karlin, I. (2005). Consistent lattice Boltzmann method. Phys. Rev. Lett. 95, 260605. doi:10.1103/physrevlett.95.260605

PubMed Abstract | CrossRef Full Text | Google Scholar

Cantisano, M., Restrepo, D., Cespedes, S., Toelke, J., Grader, A., Suhrer, M., et al. (2013). “Relative permeability in a shale formation in Colombia using digital rock physics,” in SPE/AAPG/SEG Unconventional Resources Technology Conference, Denver, Colorado, USA, August 12–14, 2013.

CrossRef Full Text | Google Scholar

Chen, D., Zhang, J., Lan, B., Wang, X., and Wang, Z. (2019). Influence of sedimentary environment on the shale gas reservoir of Sahai formation in Fuxin Basin, northeast China. Energy Sources, Part A Recovery, Util. Environ. Eff. 41, 1949–1959. doi:10.1080/15567036.2018.1548506

CrossRef Full Text | Google Scholar

Cheng, W., Jiang, G., Xie, J., Wei, Z., Zhou, Z., and Li, X. (2019). A simulation study comparing the Texas two-step and the multistage consecutive fracturing method. Petroleum Sci. 16, 1121–1133. doi:10.1007/s12182-019-0323-9

CrossRef Full Text | Google Scholar

Cherkaoui, I., Bettaibi, S., Barkaoui, A., and Kuznik, F. (2022). Magnetohydrodynamic blood flow study in stenotic coronary artery using lattice Boltzmann method. Comput. Methods Programs Biomed. 221, 106850. doi:10.1016/j.cmpb.2022.106850

PubMed Abstract | CrossRef Full Text | Google Scholar

Chikatamarla, S., and Karlin, I. (2009). Lattices for the lattice Boltzmann method. Phys. Rev. E 79, 046701. doi:10.1103/physreve.79.046701

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, D., Guan, Q., Wang, S., Huang, J., Wang, Y., and Zhang, C. (2015). Shale gas in China: Reality and dream. Energy Explor. Exploitation 33, 397–418. doi:10.1260/0144-5987.33.3.397

CrossRef Full Text | Google Scholar

Dong, H., and Blunt, M. (2009). Pore-network extraction from micro-computerized-tomography images. Phys. Rev. E 80, 036307. doi:10.1103/physreve.80.036307

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, J., Wang, X., He, S., Guo, X., Zhang, B., and Chen, X. (2017). Geochemical characteristics and source correlation of natural gas in Jurassic shales in the North Fuling area, Eastern Sichuan Basin, China. J. Petroleum Sci. Eng. 158, 284–292. doi:10.1016/j.petrol.2017.08.055

CrossRef Full Text | Google Scholar

Gunstensen, A. K., Rothman, D. H., Zaleski, S., and Zanetti, G. (1991). Lattice Boltzmann model of immiscible fluids. Phys. Rev. A 43, 4320–4327. doi:10.1103/physreva.43.4320

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, X., Hu, D., Li, Y., Liu, R., and Wang, Q. (2014). Geological features and reservoiring mode of shale gas reservoirs in Longmaxi Formation of the Jiaoshiba area. Acta Geol. Sinica-English Ed. 88, 1811–1821. doi:10.1111/1755-6724.12347

CrossRef Full Text | Google Scholar

Hao, Z., Fei, H., Hao, Q., and Liu, L. (2017). China's shale gas geological investigation and prospecting have made significant progress. Acta Geol. Sin. 91, 727–728. doi:10.1111/1755-6724.13128

CrossRef Full Text | Google Scholar

Huang, L., Khoshnood, A., and Firoozabadi, A. (2020). Swelling of Kimmeridge kerogen by normal-alkanes, naphthenes and aromatics. Fuel 267, 117155. doi:10.1016/j.fuel.2020.117155

CrossRef Full Text | Google Scholar

Huang, L., Ning, Z., Lin, H., Zhou, W., Wang, L., Zou, J., et al. (2021). High-pressure sorption of methane, ethane, and their mixtures on shales from Sichuan Basin, China. Energy and Fuels 35, 3989–3999. doi:10.1021/acs.energyfuels.0c04205

CrossRef Full Text | Google Scholar

Huang, L., Zhou, W., Xu, H., Wang, L., Zou, J., and Zhou, Q. (2021). Dynamic fluid states in organic-inorganic nanocomposite: Implications for shale gas recovery and CO2 sequestration. Chem. Eng. J. 411, 128423. doi:10.1016/j.cej.2021.128423

CrossRef Full Text | Google Scholar

Jiang, M., Li, J., and Liu, Z. (2022). A simple and efficient parallel immersed boundary-lattice Boltzmann method for fully resolved simulations of incompressible settling suspensions. Comput. Fluids 237, 105322. doi:10.1016/j.compfluid.2022.105322

CrossRef Full Text | Google Scholar

Kareem, W. A., and Asker, Z. M. (2022). Simulations of isotropic turbulent flows using lattice Boltzmann method with different forcing functions. Int. J. Mod. Phys. C 33, 2250145. doi:10.1142/s0129183122501455

CrossRef Full Text | Google Scholar

Khan, N. H., Paswan, M. K., and Hassan, M. A. (2022). Natural convection of hybrid nanofluid heat transport and entropy generation in cavity by using lattice Boltzmann method. J. Indian Chem. Soc. 99, 100344. doi:10.1016/j.jics.2022.100344

CrossRef Full Text | Google Scholar

Li, N., Xie, H., Gao, Z., and Li, C. (2022). Study on the hydraulic fracturing failure behaviour of granite and its comparison with gas fracturing. Sustainability 14, 14593. doi:10.3390/su142114593

CrossRef Full Text | Google Scholar

Li, Y., Yang, S., Zhao, W., Li, W., and Zhang, J. (2018). Experimental of hydraulic fracture propagation using fixed-point multistage fracturing in a vertical well in tight sandstone reservoir. J. Petroleum Sci. Eng. 171, 704–713. doi:10.1016/j.petrol.2018.07.080

CrossRef Full Text | Google Scholar

Lutz, B., Lewis, A., and Doyle, M. (2013). Generation, transport, and disposal of wastewater associated with Marcellus shale gas development. Water Resour. Res. 49, 647–656. doi:10.1002/wrcr.20096

CrossRef Full Text | Google Scholar

Ma, X., Wang, H., Zhou, S., Shi, Z., and Zhang, L. (2021). Deep shale gas in China: Geological characteristics and development strategies. Energy Rep. 7, 1903–1914. doi:10.1016/j.egyr.2021.03.043

CrossRef Full Text | Google Scholar

Ma, Z., Pi, G., Dong, X., and Chen, C. (2017). The situation analysis of shale gas development in China-based on structural equation modeling. Renew. Sustain. Energy Rev. 67, 1300–1307. doi:10.1016/j.rser.2016.06.085

CrossRef Full Text | Google Scholar

Nie, G. M., Guo, D., Wang, Y., and Cheng, X. (2013). Research and application of the deep horizontal drilling and completion technology for Liaohe oilfield buried hill reservoir. Appl. Mech. Mater. 241, 1396–1399. doi:10.4028/www.scientific.net/amm.241-244.1396

CrossRef Full Text | Google Scholar

Otsu, N. (1979). A threshold selection method from gray-level histograms. IEEE Trans. Syst. man, Cybern. 9, 62–66. doi:10.1109/tsmc.1979.4310076

CrossRef Full Text | Google Scholar

Qu, Z., Wang, J., Guo, T., Shen, L., Liao, H., Liu, X., et al. (2021). Optimization on fracturing fluid flowback model after hydraulic fracturing in oil well. J. Petroleum Sci. Eng. 204, 108703. doi:10.1016/j.petrol.2021.108703

CrossRef Full Text | Google Scholar

Ren, J., Wang, S., Wu, Q., and Song, Y. (2023). Axisymmetric lattice Boltzmann model for liquid flows with super-hydrophobic cylindrical surfaces. Eur. J. Mech./B Fluids 98, 120–142. doi:10.1016/j.euromechflu.2022.12.001

CrossRef Full Text | Google Scholar

Ren, J., Zheng, Q., Guo, P., Peng, S., Wang, Z., and Du, J. (2019). Pore-scale lattice Boltzmann simulation of two-component shale gas flow. J. Nat. Gas Sci. Eng. 61, 46–70. doi:10.1016/j.jngse.2018.11.011

CrossRef Full Text | Google Scholar

Rongved, M., Holt, R., and Larsen, I. (2019). The effect of heterogeneity on multiple fracture interaction and on the effect of a non-uniform perforation cluster distribution. Geomechanics Geophys. Geo-Energy Geo-Resources 5, 315–332. doi:10.1007/s40948-019-00113-4

CrossRef Full Text | Google Scholar

Shan, X., and Chen, H. (1993). Lattice Boltzmann model for simulating flows with multiple phases and components. Phys. Rev. E 47, 1815–1819. doi:10.1103/physreve.47.1815

PubMed Abstract | CrossRef Full Text | Google Scholar

Shao, J., Shu, C., Wu, J., and Chew, Y. (2013). A stencil adaptive phase-field lattice Boltzmann method for two dimensional incompressible multiphase flows. Int. J. Numer. Methods Fluids 72, 671–696. doi:10.1002/fld.3759

CrossRef Full Text | Google Scholar

Shao, J., You, L., Kang, Y., Chen, M., and Tian, J. (2022). Salinity of flowback fracturing fluid in shale reservoir and its reservoir damage: Experimental and field study. J. Petroleum Sci. Eng. 211, 110217. doi:10.1016/j.petrol.2022.110217

CrossRef Full Text | Google Scholar

Sun, H., Zhou, D., Cai, X., Wang, F., Feng, D., and Lu, T. (2020). Progress and prospects in shale gas development of Sinopec. China Pet. Explor. 25, 14. doi:10.3969/j.issn.1672-7703.2020.02.002

CrossRef Full Text | Google Scholar

Swift, M., Osborn, W., and Yeomans, J. (1995). Lattice Boltzmann simulation of nonideal fluids. Phys. Rev. Lett. 75, 830–833. doi:10.1103/physrevlett.75.830

PubMed Abstract | CrossRef Full Text | Google Scholar

Vengosh, A., Jackson, R. B., Warner, N., Darrah, T. H., and Kondash, A. (2014). A critical review of the risks to water resources from unconventional shale gas development and hydraulic fracturing in the United States. Environ. Sci. Technol. 48, 8334–8348. doi:10.1021/es405118y

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, F., Chen, Q., Lyu, X., and Zhang, S. (2020). Fracturing-fluid flowback simulation with consideration of proppant transport in hydraulically fractured shale wells. ACS omega 5, 9491–9502. doi:10.1021/acsomega.0c00714

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Liao, X., and Ding, H. (2015). Monitoring and evaluating the volume fracturing effect of horizontal well. J. Nat. gas Sci. Eng. 22, 498–502. doi:10.1016/j.jngse.2015.01.005

CrossRef Full Text | Google Scholar

Wang, S., Chen, L., Feng, Q., Chen, L., Fang, C., and Cui, R. (2023). Pore-scale simulation of gas displacement after water flooding using three-phase lattice Boltzmann method. Capillarity 6, 19–30. doi:10.46690/capi.2023.02.01

CrossRef Full Text | Google Scholar

Wang, Y., Sun, Y., Liang, D., and Cai, J. (2020). Simulation study on spontaneous imbibition of tight sandstone based on digital core and lattice Boltzmann method. Petroleum Sci. Bull. 5, 458–466. (Chinese).

Google Scholar

Wu, K., Van, D., Couples, G., Jiang, Z., Ma, J., and Sorbie, K. (2006). 3D stochastic modelling of heterogeneous porous media-applications to reservoir rocks. Transp. porous media 65, 443–467. doi:10.1007/s11242-006-0006-z

CrossRef Full Text | Google Scholar

Wu, L., Zhang, J., Jia, D., Wang, S., and Yan, Y. (2022). Performance evaluation of multistage fractured horizontal wells in tight gas reservoirs at Block M, Ordos Basin. Energies 15, 613. doi:10.3390/en15020613

CrossRef Full Text | Google Scholar

Xu, R., Yang, S., Xiao, Z., and Jin, Y. (2022). Quantitatively study on imbibition of fracturing fluid in tight sandstone reservoir under high temperature and high pressure based on NMR technology. J. Petroleum Sci. Eng. 208, 109623. doi:10.1016/j.petrol.2021.109623

CrossRef Full Text | Google Scholar

Yang, Y., Xie, L., He, B., and Zhao, P. (2022). Numerical study on the use of alternating injection hydraulic fracturing technology to optimize the interaction between hydraulic fracture and natural fracture. Front. Earth Sci. 10, 422. doi:10.3389/feart.2022.873715

CrossRef Full Text | Google Scholar

Zhang, J., Yan, G., and Dong, Y. (2009). A new lattice Boltzmann model for the Laplace equation. Appl. Math. Comput. 215, 539–547. doi:10.1016/j.amc.2009.05.047

CrossRef Full Text | Google Scholar

Zhang, L., Kang, Q. J., Yao, J., Gao, Y., Sun, H., Liu, H., et al. (2015). Pore scale simulation of liquid and gas two-phase flow based on digital core technology. Sci. China Technol. Sci. 58, 1375–1384. doi:10.1007/s11431-015-5842-z

CrossRef Full Text | Google Scholar

Zhang, R., Zheng, X., Li, S., Jin, K., Zheng, H., and Wang, L. (2018). Environmentally friendly hydraulic fracturing and water-free fracturing technologies. Int. J. Oil, Gas Coal Technol. 17, 375–390. doi:10.1504/ijogct.2018.10012095

CrossRef Full Text | Google Scholar

Zhang, T., Javadpour, F., Li, J., Zhao, Y., Zhang, L., and Li, X. (2021). Pore-scale perspective of gas/water two-phase flow in shale. SPE J. 26, 828–846. doi:10.2118/205019-pa

CrossRef Full Text | Google Scholar

Zheng, H., Liao, R., Cheng, N., and Shi, S. (2021). Microscopic mechanism of fracturing fluid imbibition in stimulated tight oil reservoir. J. Petroleum Sci. Eng. 202, 108533. doi:10.1016/j.petrol.2021.108533

CrossRef Full Text | Google Scholar

Zheng, J., Ju, Y., and Wang, M. (2018). Pore-scale modeling of spontaneous imbibition behavior in a complex shale porous structure by pseudopotential lattice Boltzmann method. J. Geophys. Res. Solid Earth 123, 9586–9600. doi:10.1029/2018jb016430

CrossRef Full Text | Google Scholar

Zheng, J., Lei, W., Ju, Y., and Wang, M. (2021). Investigation of spontaneous imbibition behavior in a 3D pore space under reservoir condition by lattice Boltzmann method. J. Geophys. Res. Solid Earth 126, e2021JB021987. doi:10.1029/2021jb021987

CrossRef Full Text | Google Scholar

Zhu, B., Cheng, H., Qiao, Y., Liu, C., Shi, Y., Zhang, K., et al. (2012). Porosity and permeability evolution and evaluation in anisotropic porosity multiscale-multiphase-multicomponent structure. Chin. Sci. Bull. 57, 320–327. doi:10.1007/s11434-011-4874-4

CrossRef Full Text | Google Scholar

Zou, C., Dong, D., Wang, Y., Li, X., Huang, J., Wang, S., et al. (2015). Shale gas in China: Characteristics, challenges and prospects (I). Petroleum Explor. Dev. 42, 753–767. doi:10.1016/S1876-3804(15)30072-0

CrossRef Full Text | Google Scholar

Zou, C., Dong, D., Wang, Y., Li, X., Huang, J., Wang, S., et al. (2016). Shale gas in China: Characteristics, challenges and prospects (II). Petroleum Explor. Dev. 43, 182–196. doi:10.1016/S1876-3804(16)30022-2

CrossRef Full Text | Google Scholar

Zou, C., Yang, Z., Zhu, R., Zhang, G., Hou, L., Wu, S., et al. (2015). Progress in China's unconventional oil and gas exploration and development and theoretical technologies. Acta Geol. Sinica-English Ed. 89, 938–971. doi:10.1111/1755-6724.12491

CrossRef Full Text | Google Scholar