Numerical Simulation Study on Temporary Well Shut-In Methods in the Development of Shale Oil Reservoirs

Abstract

Field tests indicate that temporary well shut-ins may enhance oil recovery from a shale reservoir; however, there is currently no systematic research to specifically guide such detailed operations in the field, especially for the design of the shut-in scheme and multiple rounds of shut-ins. In this study, the applicability of well shut-in operations for shale oil reservoirs is studied, and a numerical model is built using the finite element method. In order to simulate the production in a shale oil reservoir, two separate modules (i.e., Darcy’s law and phase transport) were two-way coupled together. The established model was validated by comparing its results with the analytical Buckley–Leverett equation. In this paper, the geological background and parameters of a shale oil reservoir in Chang-7 Member (Chenghao, China) were used for the analyses. The simulation results show that temporary well shut-in during production can significantly affect well performance. Implementing well shut-in could decrease the initial oil rate while decreasing the oil decline rate, which is conducive to long-term production. After continuous production for 1000 days, the oil rate with 120 days shut-in was 9.85% larger than the case with no shut-in. Besides, an optimal shut-in time has been identified as 60 days under our modeling conditions. In addition, the potential of several rounds of well shut-in operations was also tested in this study; it is recommended that one or two rounds of shut-ins be performed during development. When two rounds of shut-ins are implemented, it is recommended that the second round shut-in be performed after 300 days of production. In summary, this study reveals the feasibility of temporary well shut-in operations in the development of a shale oil reservoir and provides quantitative guidance to optimize these development scenarios.

1. Introduction

Shale oil is considered a critical unconventional source to meet future energy demands [1]. For the development of a shale oil reservoir, multistage fracturing in the horizontal wells (MFHW) has been broadly utilized to produce the trapped hydrocarbon resources [2]. Through fracturing treatment, fractures with high conductivity are created in the formation rock, which is crucial for the development of unconventional reservoirs [3,4]. Field evidence shows that, after hydraulic fracturing, temporary well shut-ins can improve well performance and enhance ultimate recovery in low-permeability unconventional reservoirs. Due to recent fluctuations in the oil price due to the COVID pandemic, a considerable number of these wells around the world were temporarily shut-in over the past few years [5,6], and this validates our conclusion. Therefore, the well shut-in has become a hot topic in the development of shale and tight oil reservoirs, especially in China. In the Changqing oil field in the northwest of China, well shut-in has been utilized extensively during the development of its shale oil reservoirs. Since 2017, there have been hundreds of new horizontal wells fractured within the Chang-7 Member in the upper Triassic Yanchang formation. All those production wells underwent a shut-in before flowback and formal production were resumed. Even though temporary well shut-in has been widely used in this region, the mechanism behind this process is still unclear to the operators.

Currently, imbibition, including dynamic and spontaneous imbibition, is considered one of the major causes and can explain the mechanism behind the well shut-in [7]. Therefore, understanding how to take an advantage of the imbibition effect to enhance well performance and achieve better long-term productivity could be substantial to the development of shale oil reservoirs [8]. In terms of the studies concerning imbibition, nuclear magnetic resonance (NMR) imaging techniques [9] are widely used in the experiments. Generally speaking, in those experiments, researchers saturated core samples with oil and then immersed the core samples into water to see the imbibition. During the experiment process, NMR testing is used to image the oil–water displacement [10,11,12,13,14,15]. Karimi et al. [16] leveraged centrifugation and NMR techniques at the same time to study the oil–water displacement pattern, focusing on the effect of capillary force. They found that, due to the imbibition effect, water enters into the tiny pores much more easily than oil in a water-wet environment. Tu and Sheng [17] studied the effect of pressure on imbibition in a shale oil reservoir, utilizing experimental and numerical methods. Based on the NMR imaging results, Cheng et al. [18] claimed that submicropores contribute more to the ultimate recovery of spontaneous and dynamic imbibition, despite the nanopores having stronger capillary forces. In addition to experimental studies on imbibition, there has been some progress made in the theoretical and numerical study of this phenomenon [19]. Schmid et al. [20,21] presented a semianalytical solution to describe spontaneous imbibition behavior. In their study, they also provided a numerical solution for this problem. Their simulation efforts considered two wettability cases and three viscosity ratio cases [21]. Besides, Khan et al. [22] simulated the oil–water displacement based on a fully implicit black oil simulator, considering strongly water-wet, weakly-wet, and mixed wet cases, respectively.

In terms of the effect of shut-in on well performance in unconventional reservoirs, studies are still relatively limited. Wang [23] studied the characteristic of reservoir energy balance and energy storage after a well shut-in. It can be concluded that a temporary well shut-in after fracturing could accelerate the diffusion and pressure propagation and slow down the reservoir energy depletion. Zhang et al. [24] evaluated the potential for oil recovery enhancement when considering imbibition and the corresponding time delay in a shale oil reservoir for a shut-in. In their study, the effects of the pore geometry and clay content were focused on. Liu [25] used NMR to study the effect of well shut-in on the imbibition rate and recovery of tuffite, shale, tight sandstone, and clastic volcanic rock. Many studies looked at the effect of well shut-in on aqueous phase trapping (APT) caused by formation damage due to drilling, completion, etc. Based on the experimental results, they found that APT could be “auto-removed” after temporary well shut-in [26,27,28] and the optimal postfrac shut-in time could be determined [29]. Wang et al. [30] presented a pressure drop model for postfracturing shut-in simulation, considering multiple effects, including fluid imbibition and oil replacements. Their research revealed that the pressure drop during shut-in can be divided into eight sequential stages, and their results can be used for the interpretation of fracture and reservoir parameters. Eltahan et al. [31] studied the impact of well shut-in after hydraulic-fracture treatments on the productivity and recovery of tight oil reservoirs. In their study, the embedded discrete fracture model (EDFM) [32] was used for fracture representation, and their results indicated that the well shut-in could improve the recovery by as much as 5%. Jia et al. [33] investigated the shut-in effect on production performance, considering stress sensitivity. They found that reservoir permeability and capillary force have the most obvious effects on shut-in performance.

Regarding the gaps in the current research on well shut-in methods, although a great amount of research has been presented, there still exist some problems to be solved in this area. First, most of the experiments mentioned before were carried out in a core scale, for which the conclusion may not be convincing when repeated at the field scale. Second, it is known that well shut-in operation can be effective for the long-term production of shale oil reservoirs. However, several key points still need to be addressed for optimizing future production scenarios, including determining the optimum shut-in interval, what is the difference between optimum shut-in time in different reservoirs, and the potential for multiple rounds of shut-ins. Finally, for The Chang-7 Member, the shale oil reservoir in the Changqing oil reservoir, extensive field tests of well shut-in have been carried out while there lacks specific guidance for field production, especially for the potential of multiple rounds of well shut-ins. This issue demands we present a specific investigation [34]. In summary, the studies mentioned before still cannot fully explain the mechanism behind the well shut-in and provides guidance on the shut-in operation in the field. Further studies are still necessary in this area.

The objective of this research is to study the feasibility of the temporary shut-in method in the development of shale oil reservoirs. Besides, we also want to design reasonable production schemes and explore the potential for multiple rounds of well shut-ins. This paper is organized as follows: First, the geological background of the research area will be described. Then, a coupled mathematical model for the simulation of a well shut-in and its numerical implementation is explained in Section 3. The numerical model is verified by comparing the results with the Buckley–Leverett equation in Section 4.1. After that, the detailed results and discussion section are presented in Section 5, including oil–water displacement during shut-in, determining the optimal well shut-in interval, and the potential for multiple rounds of well shut-ins. This study provides new insight into the optimal well shut-in time and the potential benefits of multiple rounds of well shut-ins. The results and conclusions from this paper can be expected to provide quantitative guidance to optimize the operation scenario in the development of shale oil reservoirs. Some suggestions derived from this study can be directly applied in the Chenghao area and for the Changqing oil field.

2. Geological Background

This study focuses on the shale oil developments in the Chenghao area, which is located in the Ordos Basin, northwest of China. In this area, production wells are mainly drilled to exploit shale oil in the Chang-7 Member in the upper Triassic Yanchang formation [35]. The Chang-7 Member consists primarily of profundal laminated shale, occasional tuffaceous mudstone, and muddy siltstone interbeds. It is one of the most organic-rich intervals in the central and southern parts of the Ordos Basin [36]. The Chang-7 Member in the Ordos Basin has porosities ranging from 6% to 11% and an average porosity of 8.8%. The permeability of the reservoir matrix ranges from 0.08 to 0.3 mD, and the average permeability value is 0.13 mD [35]. According to the high-pressure mercury intrusion tests carried out for this field, the pore throat radius mainly distributes from 0.02 to 1 μm. Besides, the organic matter carbon content in the matrix ranges from 0.45% to 35.85%, with an average of 9.02%. Hydrocarbon generation potential ranges from 0.19 to 116.17 mg/g, with an average value of 34.55 mg/g [37]. Since 2017, there are now hundreds of new horizontal wells fractured in the Chang-7 Member in the upper Triassic Yanchang formation. It is worth mentioning that all the production wells underwent shut-ins before flowback and production. Those newly developed wells are (on average) drilled with 1709 m of horizontal length and fracked with an injection of 28612 m³ of slickwater and 3268 m³ of proppants placed into 23 clusters with a total of 111 induced fractures [38,39,40]. Geological information of the Chang-7 Member can be found in the reference [41].

3. Mathematical Modeling

In order to describe the two-phase fluid flow in a shale oil reservoir, some idealizations and assumptions are made as follows: (1) the model is simplified into two-phase isothermal water–oil flows in a shale oil reservoir. (2) It is assumed that the permeability of the matrix is isotropic. (3) The formation is fully stimulated (penetrated) by the hydraulic fractures with finite conductivity. The fracture width is assumed to be uniform. (4) The fluid flow is presumed to be into the horizontal plane, and gravity drainage is ignored; the fluid flow in the horizontal wellbore is not taken into account. (5) There is no mass transfer/exchange happening between the water phase and the oil phase.

3.1. Pressure Formulation of Oil–Water Mixture

In this study, two independent numerical modules were utilized to calculate the pressure distribution and saturation distribution in a shale oil reservoir, respectively. For the phase transport in the shale oil reservoir, a finite element method (FEM) based model, called “Phase Transport in Porous Media” (PTPM), was utilized to calculate the saturation distribution. For the pressure distribution in the shale oil reservoir, Darcy’s law was adopted to solve the pressure distribution of the water–oil mixture. Then, the calculated pressure and saturation field are coupled together to constitute a complete numerical model for further analysis. Compared to the classical two-phase Darcy’s model, this method has some advantages, including less computational demand, lower manipulation difficulty, and good convergent performance [42].

In order to calculate the pressure distribution of the oil–water mixture, the properties of the mixture need to be generated as follows [43].

$$\overline{\rho }={{\displaystyle \sum }}_{\alpha }{\rho }_{\alpha }{S}_{\alpha },\alpha =\mathrm{o},\mathrm{w},$$

(1)

$$\overline{\mu }=\frac{\overline{\rho }}{{{\displaystyle \sum }}_{\alpha }\frac{k_{r\alpha }{\rho }_{\alpha }}{{\mu }_{\alpha }}},\alpha =\mathrm{o},\mathrm{w},$$

(2)

when α = o, this represents the oil phase, and when α = w, this indicates the water phase. $\overline{\rho }$ denotes the density of the oil–water mixture. ρ_α represents the density of the corresponding phase. S_α represents the saturation of the corresponding phase, and $\overline{\mu }$ denotes the viscosity of the oil–water mixture. k_rα represents the relative permeability for the oil and water phases. μ_α represents the viscosity of the corresponding phase.

Darcy’s law for the oil–water mixture can be expressed as follows.

$$\overline{\mathit{u}}=-\frac{k}{\overline{\mu }}\nabla \overline{p},$$

(3)

where $\overline{\mathit{u}}$ is Darcy’s velocity vector for the oil–water mixture. k is the permeability of porous media. $\overline{p}$ is the pressure of the oil–water mixture.

The continuity equation for the fluid flow can be expressed as

$$\nabla \cdot \left(\overline{\rho }\overline{\mathit{u}}\right)+\frac{\mathsf{\partial }\left(\varphi \overline{\rho }\right)}{\mathsf{\partial }t}=0,$$

(4)

where ϕ is the reservoir porosity; t is the time.

Through Equation (4), the pressure of the oil–water mixture can be calculated. However, what we really need is the pressure of the oil and water phases rather than the average pressure of the oil–water mixture. Therefore, additional steps should be taken. The pressure of the water phase p_w and the oil phase p_o can be calculated through the equation of weighted fluid pressure [43] and the capillary force equation [44].

$$\overline{p}=S_{\mathrm{o}}{p}_{\mathrm{o}}+S_{\mathrm{w}}{p}_{\mathrm{w}},$$

(5)

$$p_{\mathrm{o}}-p_{\mathrm{w}}=p_{\mathrm{c}}=p_{\mathrm{entry}}{\left(\frac{S_{\mathrm{w}}-S_{\mathrm{wi}}}{1-S_{\mathrm{wi}}-S_{\mathrm{or}}}\right)}^l,$$

(6)

where p_c is the capillary pressure. p_entry is the entry capillary pressure. S_wi is the irreducible water saturation. S_or is the residual oil saturation. l is the capillary pressure exponent. Since the saturation information would be transferred from the calculation of PTPM module, the only unknown values in Equations (5) and (6) are p_w and p_o, and thus, the pressure of oil and water phases can be solved.

3.2. Phase Transport in Porous Media

Based on the oil pressure p_w and water pressure p_o calculated in the previous section, the flow velocity for the oil phase u_o and the water phase u_w in the PTPM module can be calculated as follows.

$${\mathit{u}}_{\alpha }=-\frac{k{k}_{\mathrm{r}\alpha }}{{\mu }_{\alpha }}\nabla p_{\alpha },\alpha =\mathrm{o},\mathrm{w},$$

(7)

Then, the continuity equation for each phase can be expressed as

$$\frac{\mathsf{\partial }}{\mathsf{\partial }t}\left(\varphi {\rho }_{\alpha }{S}_{\alpha }\right)+\nabla \cdot \left({\rho }_{\alpha }{\mathit{u}}_{\alpha }\right)={\rho }_{\alpha }{q}_{\alpha },\alpha =\mathrm{o},\mathrm{w},$$

(8)

where q_α is the sink or source term in the reservoir.

Besides, an additional constraint equation needs to be added to the model, i.e.,

$$S_{\mathrm{o}}+S_{\mathrm{w}}=1,$$

(9)

Here, the only unknown variable in Equation (8) is water saturation, and thus, the saturation profile can be obtained. By coupling Darcy’s flow module and PTPM module, the pressure and saturation distributions of oil and water phases can be acquired.

3.3. Fluid Flow in Hydraulic Fractures

For the development of a shale oil reservoir with multistage hydraulic fracturing, fractures play a crucial role in the production. In terms of the simulation of fluid flow inside the hydraulic fractures, here a discrete fracture model [45] is adopted.

The kinematic equation in the discrete fracture model can be written as follows in order to replace Equation (3).

$${\overline{\mathit{u}}}_{\mathrm{f}}=-\frac{k_{\mathrm{f}}}{\overline{\mu }}{\nabla }_{\mathrm{T}}{\overline{p}}_{\mathrm{f}},$$

(10)

where ${\overline{\mathit{u}}}_{\mathrm{f}}$ is the fluid velocity for the oil–water mixture inside the fractures. ${\nabla }_{\mathrm{T}}{\overline{p}}_{\mathrm{f}}$ is the pressure gradient tangent to the fracture surface [46]. k_f is the permeability of the fracture.

The continuity equation for fluid flow along the hydraulic fractures is used to replace Equation (4), which is

$${\nabla }_{\mathrm{T}}\cdot \left(d_{\mathrm{f}}\overline{\rho }{\overline{\mathit{u}}}_{\mathrm{f}}\right)+d_{\mathrm{f}}\frac{\mathsf{\partial }}{\mathsf{\partial }t}\left({\varphi }_{\mathrm{f}}\overline{\rho }\right)=0,$$

(11)

where d_f is the fracture width. ϕ_f is the porosity of the propped fracture.

3.4. Auxiliary Boundary Conditions

The outer boundaries of the model, ∂Ω₁, are assumed to form a closed Neumann boundary condition as follows [47]

$$-\mathit{n}\cdot \rho \overline{\mathit{u}}|\mathsf{\partial }{\Omega }_1=0,$$

(12)

where n is the normal vector to the outer boundary ∂Ω_1.

The simulation task can be divided into several steps, including (1) water injection, (2) well shut-in, and (3) production. For each step, different well treatments are implemented. In this study, we are assuming the horizontal well is infinitely conductive, while the fractures have a finite fracture conductivity.

(1)

For the water injection step, the boundary condition ∂Ω₂, i.e., the perforation holes, for a horizontal well is defined as an inhomogeneous Neumann condition with a flow rate.

$$-\mathit{n}\cdot \rho \overline{\mathit{u}}|\mathsf{\partial }{\Omega }_2=\overline{\rho }{\mathit{u}}_0,$$

(13)

$$S_{\mathrm{w}}|\mathsf{\partial }{\Omega }_2=1,$$

(14)

where u₀ is an input parameter to represent the injected water influx. Here we distribute the total injection rate average to each fracture.

(2)

For the shut-in step, the boundary condition ∂Ω₂ for the horizontal well is defined as the homogeneous Neumann condition or closed boundary.

$$-\mathit{n}\cdot \rho \overline{\mathit{u}}|\mathsf{\partial }{\Omega }_2=0,$$

(15)

(3)

For the production step, the boundary condition ∂Ω₂ again for a horizontal well is defined as Dirichlet boundary condition with a bottomhole pressure.

$$\overline{p}|\mathsf{\partial }{\Omega }_2=p_{\mathrm{w}},$$

(16)

where p_w is an input parameter to represent the bottomhole pressure at the horizontal well. The initial conditions for the reservoir would be

$${\overline{p}|}_{t=0}=p_{\mathrm{i}},$$

(17)

$${S_{\mathrm{w}}|}_{t=0}=S_{\mathrm{i}},$$

(18)

where p_i is the initial reservoir pressure; S_i is the initial reservoir water saturation.

In order to achieve numerical stability and avoid any numerical diffusion/oscillation, detailed numerical settings need to be designed in advance. In this study, the coupling of fluid flow and phase transfer in the reservoir should be especially highlighted. In terms of the coupling of PTPM for the saturation field and Darcy’s law for the pressure field, the Segregated algorithm is recommended in this model rather than using Fully Coupled algorithms, considering the numerical performance [48]. The Fully Coupled approach forms a single large system of equations that solve for all of the unknowns (the fields) and includes all of the couplings between the unknowns (the multiphysics effects) at once, within a single iteration. On the other hand, the Segregated approach will not solve all of the unknowns at one time. Instead, it subdivides the problem into two or more segregated steps. Each step will usually represent a single physics, but sometimes even a single physics can be subdivided into steps, and sometimes one step can contain multiple physics. These individual segregated steps are smaller than the full system of equations that are formed with the Fully Coupled approach. The segregated steps are solved sequentially within a single iteration, and thus, less memory is required. Since the Fully Coupled approach includes all coupling terms between the unknowns, it often converges more robustly and in fewer iterations as compared to the Segregated approach. However, each iteration will require relatively more memory and time to solve, so the Segregated approach can be faster overall. In the Segregated algorithm, the constraint for the saturation range, i.e., 0–1, can be easily prescribed mathematically using its value-limit function. In addition, an Anderson acceleration method is also used to improve numerical stability [49].

The calculation process is divided into three steps: the water injection process, well shut-in process, and the production process, respectively. Therefore, it is necessary to consider the delivery of the simulation results. When the calculation of the water injection study is terminated, the simulation results in the last time step should be transported into the well shut-in study as the initial value by using the “study reference” interface in COMSOL Multiphysics. The only change between these two steps is the boundary condition. Then, the same process needs to be set for the other studies. By doing so, different physical processes with different boundary conditions and numerical settings could be connected together, as shown in Figure 1.

4. Numerical Implementation

4.1. Model Validation

In order to validate the numerical model, a simulation case is built to compare its calculation results with the analytical solutions. For two-phase flow, the most widely used analytical model is the Buckley–Leverett model [50]. This model estimates the advance of a fluid displacement front in an immiscible displacement process. According to the Buckley–Leverett equation, the distance x traveled by a particular S1 contour can be expressed as the following.

$$x-x_0=\frac{\left(\frac{d{f}_1}{d{S}_1}\right)}{A\cdot \varphi }\cdot {{\displaystyle \int }}_0^{t}q\left(t\right)dt,$$

(19)

where x₀ is the radius of the wellbore. f₁ is the fractional flow of phase 1. A is the cross-section in the flow direction. q is the injection rate. ϕ is the porosity of the reservoir matrix. t represents time.

In this validation section, we assume that the relative permeability of phase 1 is expressed as k_r1 = S₁², and the relative permeability of phase 2 is expressed as k_r2 = (1 − S₁)². Then, the fraction flow of phase 1 can be derived as follows.

$$f_1=\frac{1}{1+\frac{S_1{}^2{\mu }_1}{{\left(1-S_1\right)}^2{\mu }_2}},$$

(20)

Then, the derivative of f₁ can be expressed as follows.

$$\frac{d{f}_1}{d{S}_1}=\frac{\frac{2S_1{\mu }_1}{{\left(1-S_1\right)}^2{\mu }_2}+\frac{2S_1{}^2{\mu }_1}{{\left(1-S_1\right)}^3{\mu }_2}}{{\left(1+\frac{S_1{}^2{\mu }_1}{{\left(1-S_1\right)}^3{\mu }_2}\right)}^2},$$

(21)

Therefore, the Buckley–Leverett equation utilized in this section can be derived as follows.

$$x-x_0=\frac{\left(\frac{2S_1{\mu }_1}{{\left(1-S_1\right)}^2{\mu }_2}+\frac{2S_1{}^2{\mu }_1}{{\left(1-S_1\right)}^3{\mu }_2}\right)\cdot q\cdot t}{{\left(1+\frac{S_1{}^2{\mu }_1}{{\left(1-S_1\right)}^3{\mu }_2}\right)}^2\cdot A\cdot \varphi },$$

(22)

The parameters for the numerical simulation and analytical calculation are as shown in

Table 1. Parameters for the validation of the numerical model [51].

Parameter	Value	Unit
Injection rate of phase 1	1 × 10−6	m/s
Area of cross section	1	m2
Viscosity of phase 1	1	mPa·s
Viscosity of phase 2	1	mPa·s
Porosity of reservoir matrix	15	%
Total calculation time	300	days
Permeability of reservoir matrix	100	10−3 μm2
Initial saturation of phase 2 in reservoir	100	%

.

The numerical model is validated by comparing its results with the analytical Buckley–Leverett equation. Figure 2 shows the comparison of the flood front curve between the numerical model and the analytical model. It can be seen from Figure 2 that the numerical model runs strictly according to the setting regimes. In general, the numerical model has relatively high accuracy.

4.2. Description of the Simulation Model

In the development of a shale oil reservoir, multistage hydraulic fracturing technology is widely adopted. After the hydraulic fracturing, the reservoir matrix area can be mainly divided into two separate flow regions, which are the stimulated reservoir area (SRV) and the unstimulated reservoir matrix area, as shown in Figure 3. Based on the mathematical model presented in Section 3, a numerical model for the shale oil reservoir was established for the simulations.

Due to the symmetry of the physical model, the model geometry can be simplified to reduce the computational quantity by only calculating half of the reservoir area, as shown in Figure 4a. The horizontal well is presumed to be located at the center of the reservoir, which is the bottom of the geometry model after the simplification and is represented using blue lines, as shown in Figure 4a. The main hydraulic fractures are distributed alongside the horizontal well. Unstructured grids were generated for model discretization, as shown in Figure 4b. The mesh grids near the wellbore and hydraulic fractures are refined significantly, as shown in the figure.

The parameters for the simulation of the multistage fracturing horizontal well in a shale oil reservoir are derived from a typical production well (H-1 in The Chang-7 Member, Chenghao area, northwest of China). The detailed parameter settings are as listed in

Table 2. Parameters for the simulation model of a multistage fracturing horizontal well in a shale oil reservoir.

Parameter	Value	Unit
Reservoir area after symmetrical processing	2350 × 1000	m2
Main hydraulic fracture width	0.5	cm
Permeability of main hydraulic fractures	1800	10−3 μm2
Permeability of unstimulated reservoir area	0.2	10−3 μm2
Permeability of stimulated reservoir area	50	10−3 μm2
Initial porosity of the unstimulated reservoir area	7.62	%
Initial porosity of the stimulated reservoir area	10.4	%
Initial porosity of the main hydraulic fractures	15.3	%
Initial reservoir pressure	2.5 × 107	Pa
Capillary force at Sw = 0.5	1.87 × 106	Pa
Production pressure of the horizontal well	1.5 × 107	Pa
Half-length of the main hydraulic fractures	150	m
Initial oil density	850	kg/m3
Oil viscosity	1	mPa·s
Water density	1000	kg/m3
Water viscosity	0.2	mPa·s
Comprehensive compressibility coefficient for matrix area	5 × 10−9	1/Pa
Comprehensive compressibility coefficient for fractures	1.0 × 10−8	1/Pa
Oil saturation range	20–80	%
Calculation time of well shut-in	60	days
Calculation time of production	1000	days

.

The model is discretized using finite element methods. The governing equations for the reservoir domain and the fracture domain are formulated separately and were coupled at COMSOL’s interface for the subsurface flow module and the discrete fracture model. Detailed numerical setting details can be found in

Table 3. Parameters for the numerical setting.

Parameter	Value	Unit
Time-dependent solver	MUMPS [52]	#
Time stepping method	Backward Differentiation Formulas (BDF)	#
Maximum BDF order	1	#
Minimum BDF order	5	#
Event tolerance	0.01	#
Consistence initialization	Backward Euler	#
Fraction of initial step for Backward Euler	0.01	#
Initial time step	0.001	day
Maximum time step	10	day

The symbol # means null.

. By implementing the presented numerical settings and mesh generation scenarios, the numerical simulations of well shut-in operation in a shale oil reservoir are conducted. The convergence plot is shown in Figure 5. One may see from the figure that the model solution converges fast and smoothly for the given parameters.

5. Results and Discussion

5.1. Analysis of Oil–Water Displacement during the Well Shut-In Process

After the water injection, a well shut-in interval of 60 days was considered in the numerical model. The oil saturation distribution at the beginning time of the shut-in is shown in Figure 6a. The distribution of oil saturation after a continuous shut-in for 60 days is shown in Figure 6b.

In Figure 6a, we see that there is a clear boundary between the reservoir oil and the front edge of the injected water. By contrast, it can be seen in Figure 6b that the boundary between the oil in the reservoir and the injected water becomes less distinctive or even vanishes after shut-in. In this process, the injected water moves deep into the reservoir matrix due to capillary forces and creates higher water pressure near the (horizontal) well.

A magnified comparison between Figure 6a,b is presented in Figure 7 to reveal the effect of the shut-in. Figure 7 indicates the oil saturation distribution near the wellbore and fractures for the 1st day and 60th day of the shut-in. According to the comparison, we can see that an oil–water displacement is clear. After a continuous shut-in for 60 days, the isosaturation line, S_o = 0.75, moves from y = 203 m to y = 251 m, which means the injected water is imbibed into the reservoir matrix. At the same time, the isosaturation line, S_o = 0.57, moves back slightly, from 181 m to 169 m, which means some oil in the reservoir is displaced to areas near the wellbore and fractures. These characteristic movements of the isosaturation lines indicate oil–water displacement during the shut-in due to a dynamic imbibition mechanism.

5.2. Effect of Shut-In Time on Well Performance

In the last section, we compared the distribution pattern of oil saturation during shut-ins. In this section, we are going to study the effect of the duration of shut-in on well performance. Five simulation cases with different shut-in intervals are implemented for comparison. The shut-in time in those cases varies from no shut-in to 120 days.

The daily water rate for the first 1000 days after shut-in is shown in Figure 8. It can be seen that as the shut-in time increases, the decline rate of the daily water rate decreases. At the same time, as the shut-in time increases, the change in the decline rate becomes insignificant, and it is difficult to distinguish the difference between the water rate after 90 days of shut-in versus 120 days of shut-in.

The daily oil production rate for the first 1000 days after shut-in (for different shut-in times) is plotted in Figure 9. It can be seen that, in the case of no shut-in, the oil production rate not only rises very fast at the beginning of the production process but also decreases severely after reaching the peak. In this simulation scenario, the reservoir energy is depleting quickly, which is detrimental to long-term oil production. In contrast, the oil rate curves with well shut-in vary with time at a relatively lower rate. The longer the shut-in time, the lower the rate for the curve. Besides, as the shut-in time increases, the decline rate of the oil production rate decreases. After continuous production for 1000 days, the oil production rate with 120 days shut-in was 9.85 % greater than that with no shut-in operation, which is considerable and should not be neglected when considering the long-term production of shale oil reservoirs. From Figure 9, We can also conclude that the temporary well shut-in method might lead to productivity loss at the early stage of production. However, the implementation of wells shut-in will decrease the production decline rate (or slow down the energy depletion) and increase the oil production rate after long-term production. As a result, we also need to consider the economic efficiency of the well shut-in method, given the estimated well production lifetime.

Besides,

Table 4. Peak production rate and their timings according to five simulation scenarios.

Shut-In Time (Day)	Oil Rate Peak (t/Day)	Timing (Day)
0	23.33	32
30	18.96	63
60	17.09	84
90	16.32	93
120	15.94	97

summarizes the information on oil production peak for different scenarios described in Figure 9. From the table, we can see that the peak of oil production and its timing for different scenarios vary with the shut-in time. As shut-in time increases from 0–30 days, the value of the oil rate peak decreases dramatically from 23.33 t/day to 18.96 t/day. Then, this change gradually sees a lower rate. As the shut-in time increases from 30 days to 60 days, the well productivity peak decreases from 18.96 t/day to 17.09 t/day. Then, as this time changes from 60 days to 120 days, the peak value drops slightly from 17.09 t/day to 15.94 t/day.

We can also look at the cumulative oil production rate for this period to decipher the importance of shut-in time as shown in Figure 10. As shut-in time increases, the cumulative oil production in the first 1000 days increases dramatically. Then, beyond 60 days shut-in time, cumulative oil production only sees a miniscule improvement. According to Figure 10, a positive relationship between the cumulative oil rate and well shut-in time can be identified up to a specific shut-in time.

However, it is worth noting that the positive correlation between cumulative oil production and shut-in time does not mean that a longer shut-in time leads to better development efficiency. This is because different shut-in scenarios imply different operation times. Basically, if there is no shut-in, the total operation time is 1000 days. However, if the well shut-in time is 30 days, the total time would be 1030 days. When we calculate the average daily oil rate, we need to consider the time cost and calculate cumulative oil production by dividing it by the total days of operation. This approach may draw different conclusions from the analysis. The operation time for different simulation scenarios is shown in

Table 5. Average production rates for different simulation cases.

Well Shut-In Time (Days)	0	30	60	90	120
Cumulative oil production (t)	7031	7657	7891	7925	7936
Total operation time (days)	1000	1030	1060	1090	1120
Average daily oil rate (t/day)	7.031	7.434	7.445	7.271	7.086

. The average daily oil rate according to different well shut-in times is shown in Figure 11.

It could be seen from Figure 11 that there exists an optimal well shut-in time, which is 60 days for the given parameters. In other words, for this specific production well, a well shut-in operation of 60 days could achieve the best development efficiency for the exploitation of the shale oil reservoir for the first 1000 days of production. To sum up, well shut-in operation could generate better development efficiency for the production well. A reasonable shut-in operation might be significantly beneficial to the long-term production in shale oil reservoirs.

5.3. Potential for Multiple Rounds of Shut-In

According to the field tests and calculation results from Section 5.2, well shut-in has proved to be an effective method for enhancing well performance in shale oil reservoirs. Currently, most of the field tests in China are only carried out with a single round of shut-in to enhance well performance. Multiple rounds of shut-ins have never been tested yet. As a result, in this section, we are going to explore the potential of multiple rounds of shut-ins in the development of shale oil reservoirs. Simulation scenarios with one, two, and three rounds of well shut-in operations are shown in

Table 6. Operation design for simulation scenarios with multiple rounds of shut-ins.

Rounds	Scenario Description
1	Shut-in 60 days → Produce 940 days (Total 1000 days)
2	Shut-in 60 days → Produce 300 days → Shut-in 60 days → Produce 580 days (Total 1000 days)
3	Shut-in 60 days → Produce 300 days → Shut-in 60 days → Produce 300 days → Shut-in 60 days → Produce 220 days (Total 1000 days)

. Other simulation parameters are shown in Table 2. For these three simulation scenarios, a total operation time is set as 1000 days which is consistent with the real conditions considered by the operators of unconventional wells. It should be noted that the conclusion might be slightly different if a different total operation time is incorporated.

The daily oil rate with one, two, and three rounds of well shut-ins are shown in Figure 12, Figure 13 and Figure 14.

The red lines in the above figures represent the oil production rate without shut-in, which are provided for comparative purposes. Based on the comparison between Figure 12 and Figure 13, it could be known that the second round of shut-ins could temporally enhance the oil production rate and reduce the decline rate. This change can be credited to the dynamic imbibition effect during the shut-in. By contrast, based on the comparison between Figure 13 and Figure 14, it can be seen that the effect of the three rounds of shut-ins is limited. This is because after 700 days, a great amount of fracturing water around the wellbore has flowed back and the effect of the dynamic imbibition is limited. The comparison of average daily oil rates after 1000 days, including shut-in and production, is shown in Figure 15. According to the figure, when there is no shut-in, the average oil rate after 1000 days is 7.031 t/day. When two rounds of shut-ins are implemented, the average oil rate reaches its peak, 7.434 t/day, which is 5.75 % greater than the case with no shut-in. To sum up, under the simulation conditions of this paper, the oil production of one and two rounds of well shut-in operations are close to each other. Therefore, one or two rounds of shut-ins are equally beneficial and are recommended here.

Based on the calculated results in this section, we have already shown that multiple rounds of well shut-in might significantly benefit reservoir development. However, assuming two rounds of shut-ins are implemented in the operation process, the timing for the second round of shut-in would be of vital importance for production. Since there is no related discussion in previous studies, in this section, a simulation case with different timings for the second round shut-in, i.e., the time length of the first round of production is designed for further comparison and discussion. Different production times for the first round, i.e., 200 days, 300 days, 400 days, and 500 days, are selected for the simulations.

Figure 16, Figure 17, Figure 18 and Figure 19 show the daily oil rate with different times (length) for first round of production, i.e., the timing of the second round of shut-in, respectively. It can be seen from the figures that, after the second round of shut-in, the daily oil rate would slightly rise again. However, as the time length of the first round of production increases, this production rise would diminish accordingly. When there are 500 days of production before the second round of shut-in, this rise in production rate could hardly be identified and under this condition, the second round of shut-in would have few positive effects on the oil production. Therefore, the timing of the shut-in is of vital importance for the reservoir development when multiple rounds of shut-ins are implemented.

Figure 20 shows a comparison between the average daily oil rate and different lengths of production for the first round. According to the figure, when the time length of the first round of production is 200 days, the average oil rate in 1000 days is 7.377 t/day; for 300 days, the average oil rate reaches its peak, i.e., 7.434 t/day. Then, as the time length keeps increasing, the average oil rate starts to decrease rapidly. To sum up, under the simulation conditions of this paper, it is recommended to perform 200 to 300 days of production before the second round of shut-ins.

6. Conclusions

In this study, the feasibility of the temporary well shut-in in a shale oil reservoir was studied using numerical modeling. A mathematical model for the simulation of immiscible two-phase flow in a shale oil reservoir was presented based on the coupling of the phase transport in the porous media module (PTPM) and Darcy’s law. The model was validated by comparing its results with an available analytical solution. The geological background of the Chang-7 Member (Chenghao, China) is presented, and the corresponding parameters were built into the establishment of the model simulation. Based on the simulation results, a quantitative analysis was carried out, and several conclusions can be drawn.

According to the calculated saturation distribution, an oil–water displacement can be easily identified during well shut-ins, which reveals the effect of the dynamic imbibition. Based on the calculated production rates, it can be concluded that the implementation of a shut-in can decrease the initial oil rate; however, this decreases the oil decline rate as well, which is beneficial to long-term production. After 1000 days of production, the oil production rate, with 120 days of shut-ins, was 9.85 % larger than that with no shut-in operation. By judging the average daily production rates, the optimal shut-in time was determined as 60 days for our given field conditions. Besides, the potential of several rounds of shut-ins was also explored. Single and double rounds of shut-ins are equally beneficial to long-term production under the presented reservoir condition. When two rounds of shut-ins are implemented, it is recommended to perform the second shut-in round after 300 days of production. To sum up, this study reveals a workflow for feasibly studying temporary well shut-in operations in any shale oil reservoir and provides guidance for optimizing overall development scenarios.