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.

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 porescale description of the coupled imbibition-flowback behavior is crucial to understand the flowback mechanism and its impacts. In this paper, a pseudopotential 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.

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

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 Frontiers in Earth Science frontiersin.org 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: x is the spatial coordinate, e α is the lattice speed in α direction, Δt is the time step, t is the time, Ω k coll is the collision operator of the k component, and Ω k forces is the external force operator of the k component.
The collision operator Ω k coll of the k component can be expressed as: Δx/Δt is the lattice sound velocity with c 1 in this study, and u eq represents the equilibrium velocity.
The conversion matrix M is expressed as The external force operator Ω k forces of the k component can be expressed as: Let f k α f k α − Δt 2 f F,k α , then Eq. 1 can be rewritten as: Fluid density ρ k , fluid velocity u, and fluid viscosity μ k are defined as follows: In the pseudo-potential model, the total force F k of the k component mainly includes the gas-liquid force F c,k and the fluid-solid force F ads,k , which are defined as: 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;

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 × 130 × 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).

Frontiers in Earth Science
frontiersin.org   Frontiers in Earth Science frontiersin.org 04 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 80MPa, the reservoir temperature T 130°C, and the contact angle θ 10°. Based on National Institute of Standards and Technology (NIST) chemical database, some physical parameters of gas and water under p 80MPa; T 130°C are obtained as follows: the methane density ρ g 259kg/m 3 , the water density ρ w 972kg/m 3 , the methane viscosity μ g 33.9μPa · s, the water viscosity μ w 233μPa · s, and the gas-water interfacial tension σ 40mN/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 × 130 × 130 lattices. The pressure boundary condition (i.e., the non-equilibrium bounce-back boundary condition) is applied to the inlet and outlet, and the bounceback 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 gasliquid 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 6 shows the slice diagrams of gas-liquid distribution at the equilibrium state (t 6.0 × 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.

Microcosmic flowback mechanism of fracturing fluid
The flowback of fracturing fluid is simulated by LBM. the gasliquid distribution in the digital core at the initial moment of the

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

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.

Frontiers in Earth Science
frontiersin.org 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,

FIGURE 6
Slice diagrams of gas-liquid distribution at the equilibrium state (t 6.0 × 10 6 ) in the imbibition process.

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

Frontiers in Earth Science
frontiersin.org 06 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 9 shows slice diagrams of gas-liquid distribution at the equilibrium state (t 4.0 × 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 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
Slice diagrams of gas-liquid distribution at the equilibrium state (t 4.0 × 10 6 ) in the flowback process.

Frontiers in Earth Science
frontiersin.org 07 breaks through in the middle of the macro-pores, and the fracturing fluid in small pores and blind ends is difficult to flowback. 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.

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.

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

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

Frontiers in Earth Science
frontiersin.org   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 × 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 θ 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 θ 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 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.
In the simulation of the flowback process under contact angles θ 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 × 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 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 × 10 6 ), the flowback rates for contact angles θ 10°, 30°, 60°are 21.03%、 26.98%、45.31%, respectively, which illustrates that the

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.

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

Frontiers in Earth Science
frontiersin.org wettability has a great influence on the flowback of the fracturing fluid.

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 Δp =0.5 MPa、 1.5 MPa、3.0MPa, respectively. Figure 15 shows the gas-liquid distribution at the equilibrium state (t 4.0 × 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 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 × 10 6 ), the flowback rates for flowback pressure differences Δ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.

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 −4 mD, respectively. The porosity and permeability of Core 2 are 3.90% and 2.3×10 −4 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 18 shows the gas-liquid distribution at the equilibrium state (t 6.0 × 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
Gas-liquid distribution at the equilibrium state (t 6.0 × 10 6 ) in the imbibition process under different pore structures. (A) Core 1, (B) Core 2.

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

Frontiers in Earth Science
frontiersin.org

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.
Frontiers in Earth Science frontiersin.org 12 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 20 shows the gas-liquid distribution at the equilibrium state (t 4.0 × 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 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 × 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 22 shows the gas-liquid distribution at the equilibrium state (t 6.0 × 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 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 24 shows the gas-liquid distribution at the equilibrium state (t 4.0 × 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 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 × 10 6 ), the flowback rates for fracturing fluid viscosity μ 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 23
The water saturation versus time steps in the imbibition process under different fracturing fluid viscosities.

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.

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.