An Effective Numerical Simulation Method for Steam Injection Assisted In Situ Recovery of Oil Shale

Abstract

This paper presents an effective numerical simulation method for production prediction of in situ recovery of oil shale reservoirs with steam injection. In this method, finite volume-based discretization schemes of heat and mass transfer equations of the thermal compositional model are derived and used. The embedded discrete fracture model is used to accurately handle the fractured vertical well. A smooth non-linear solver is proposed to solve the global equations, then cell pressure, temperature, saturation, component mole fractions, and well production rates can be obtained. Compared with the existing commercial software, this new method can have a smoother non-linear solution and handle the complex fracture geometry theoretically. A numerical example is used to test this presented method and can realize accurate calculation results compared with CMG. Another numerical case with a hydraulic fracture and an open thermal boundary condition is implemented to validate the presented method and can effectively handle the actual situation of steam injection-assisted in situ recovery of oil shale, which was difficult to handle using previous methods.

1. Introduction

The world is rich in oil shale resources, and against the background of the increasing shortage of traditional oil resources, the efficient exploitation of oil shale becomes increasingly important. The kinetic reaction of oil shale is the key factor for oil shale production, which had been widely modeled [1,2,3]. Up to now, the relatively mature oil-shale in situ production technologies have been microwave heating technology [4], electric heating technology, and fluid heating technology.

Since the Shell company successfully developed electric heating technology and conducted field tests, the technology has been widely recognized by the public, and the numerical simulation analysis of the technology has been widely studied. Harold and Scott [5] used the commercial software Computer Modeling Group (CMG) STARS module to simulate the electric-heating in situ production process in an oil field in northwest Colorado, simulated and calculated the oil and gas production, and compared it with the actual production data. Fan et al. [6] simulated the process of oil shale pyrolysis and pyrolysis product seepage by electric heater heating oil shale based on the Stanford General Purpose Research simulator (GPRS). Shen et al. [7] also used the stars module of CMG, considered the wellbore reflux heat loss and skin effect, simulated the field process of shell oil shale electric heating, and achieved higher calculation accuracy of production data.

Compared with the electric-heating in situ recovery technology, the fluid-assisted in situ recovery technology of oil shale is more complex [8]. More factors need to be considered in its numerical simulation, including heat transfer, temperature field, pressure field, fluid field, and the influence of fractures. Pei et al. [9] studied the mechanism of nitrogen injection-assisted in situ recovery of oil shale. The steam injection can provide more heat to the formation, and compared with the only heat conduction effect in electric heating, the thermal convection effect of steam injection in the formation can achieve more efficient heat transfer, so as to effectively improve the oil shale recovery in theory. Wang et al. [10] made a numerical study of the actual field case of in situ recovery of oil shale by injecting superheated steam. Due to the low permeability of the oil shale matrix, steam injection is often carried out by fractured vertical wells, so handling the hydraulic fracture is a key factor in steam injection-assisted oil shale recovery. Taheri et al. [11] studied the effect of hydraulic fracture characteristics on the thermal EOR process. The embedded discrete fracture model (EDFM) [12,13,14] or distinct element method [15] is an effective method to efficiently handle the large-scale discrete fracture.

A compositional model is generally used to model the in situ recovery of oil shale. In the compositional model by Fan et al. [6], the chemical reaction of kerogen decomposition, component phase equilibrium, multiphase seepage and thermal convection diffusion are considered, and the equations are solved by a fully implicit Newton iterative method. Li et al. [16,17] proposed a multi-scale numerical simulation method for oil shale in situ mining in view of the slow calculation efficiency of oil shale numerical simulation. Zheng et al. [18] used a simplified compositional model and combined the experimental data for a kinetic reaction of oil shale to study the production process of in situ combustion of oil shale. Multi-physics coupled simulation studies for oil shale recovery were also implemented to make the production prediction more accurate [19,20].

At present, there is no relevant research on the numerical model of steam injection via a fractured vertical well for oil shale production based on a thermal compositional model, so the numerical simulation of oil shale in situ production is still undergoing further research. Based on the thermal compositional model, the finite volume method, embedded discrete fracture model, and a smooth non-linear solver of the thermal compositional model without judging phase change, this paper intends to form an effective numerical simulation method of steam injection assisted in situ recovery of oil shale, to provide an efficient and accurate production prediction.

2. Methodology

2.1. Basic Thermal Compositional Model for Oil Shale Recovery

The numerical simulation of in situ recovery of oil shale involves the coupling calculation of heat and mass transfer equation based on a kerogen pyrolysis chemical reaction, so it is a more complex thermal compositional model. It mainly includes the mass conservation equation, energy conservation equation of each component, and the auxiliary equation composed of a phase equilibrium equation and well equation. It should be pointed out that in order to reduce the difficulty of calculation and improve the calculation efficiency, constant values that do not change with the component concentration are adopted for the physical properties of each phase in this paper, and it can be understood as the average of these physical properties of each phase in the whole calculation process. The relevant equations are introduced as follows:

(1)

Mass conservation equations

For fluid components:

$$-\nabla \cdot {\displaystyle \sum _j\left({\rho }_j{X}_{ij}{u}_j\right)+{\displaystyle \sum _l\left({{\nu }^{\prime }}_{i,l}-{\nu }_{i,l}\right)r_l}+{\displaystyle \sum _j\left({\rho }_j{X}_{ij}{q}_j{}^W\right)}}=\frac{\partial }{\partial t}\left(\varphi {\displaystyle \sum _j\left(S_j{\rho }_j{X}_{ij}\right)}\right)$$

(1)

where fluid components include six types: CO₂, N₂, HO₂, C₂, C₁₃, C₃₇, and $u_j=-k\frac{k_{rj}}{{\mu }_j}\left(\nabla p_j-{\gamma }_j\nabla D\right)$.

For the components in the solid phase, there are two components, kerogen, and coke. For kerogen, the mass conservation equation is:

$${\displaystyle \sum _l\left({{\nu }^{\prime }}_{i,l}-{\nu }_{i,l}\right)r_l}=\frac{\partial C_i}{\partial t}$$

(2)

where $C_i=\varphi S_j{\rho }_j{X}_{ij}$.

For coke, due to coke being generally not the product of the chemical reaction, the equation is:

$${\displaystyle \sum _i{\displaystyle \sum _l\left({\nu }_{i,l}-{{\nu }^{\prime }}_{i,l}\right)r_l}}=\frac{\partial C_C}{\partial t}l$$

(3)

where $C_C$ is the coke concentration in the solid phase.

(2)

Energy conservation equation:

$$\nabla \cdot \left(\lambda \nabla T\right)-\nabla \cdot \left[{\displaystyle \sum _j{\rho }_j{H}_j{u}_j}\right]+{\displaystyle \sum _l{H}_l{r}_l}+{\displaystyle \sum _j\left(H_j{\rho }_j{q}_j\right)}+q_H=\frac{\partial }{\partial t}\left(\varphi {\displaystyle \sum _j{S}_j{\rho }_j{U}_j}+\left(1-\varphi \right)U_r\right)$$

(4)

(3)

Auxiliary equation

Chemical reaction model:

$$r_l=K_l\exp \left(-\frac{E_l}{RT}\right)C_{kero}$$

(5)

where $C_{kero}$ is the kerogen concentration in the solid phase, because the kerogen is generally the only reactant.

In the in situ recovery of oil shale, the phase equilibrium of each fluid component in oil phase and gas phase is considered.

Phase equilibrium model:

$$X_{i,g}=K_i{X}_{i,o}$$

(6)

Well model:

$$q_j{}^W=WI\left(p_j-p_{wf}\right)$$

(7)

where $WI=\frac{2\pi k\Delta z}{\mu \ln \left(r_0/r_w\right)}$.

$$q_H{}^W=HI\left(T_h-T_{wb}\right)$$

(8)

where $HI=\frac{2\pi k\Delta z}{\ln \left(r_0/r_w\right)}$.

2.2. Finite Volume-Based Discretization of Governing Equations

In this paper, the above equations are discretized based on the block-center finite volume method. Taking the mass conservation equation of fluid components as an example, the finite volume discrete scheme of the governing equation is derived in detail as follows:

For the control volume V_h of the h-th grid, assuming that there are n_m grids adjacent to the grid, then take the time and volume integration of both sides of Equation (1) to obtain:

$$\begin{array}{l} {\displaystyle {\int }_t^{t+\Delta t}{\displaystyle {\int }_{V_h}{\displaystyle \sum _l\left({{\nu }^{\prime }}_{i,l}-{\nu }_{i,l}\right)r_l}d\Omega dt}}+{\displaystyle {\int }_t^{t+\Delta t}{\displaystyle {\int }_{V_h}\nabla \cdot {\displaystyle \sum _j\left({\rho }_j{X}_{ij}{u}_j\right)d\Omega }dt}}\\ +{\displaystyle {\int }_t^{t+\Delta t}{\displaystyle {\int }_{V_h}{\displaystyle \sum _j\left({\rho }_j{X}_{ij}{q}_j{}^W\right)d\Omega }}dt}={\displaystyle {\int }_{V_h}{\displaystyle {\int }_t^{t+\Delta t}\frac{\partial }{\partial t}\left(\varphi {\displaystyle \sum _j\left(S_j{\rho }_j{X}_{ij}\right)}\right)dt}d\Omega }\end{array}$$

(9)

The right side of the above equation is rewritten as:

$$\begin{array}{l}{\displaystyle {\int }_{V_h}{\displaystyle {\int }_t^{t+\Delta t}\frac{\partial }{\partial t}\left(\varphi {\displaystyle \sum _j\left(S_j{\rho }_j{X}_{ij}\right)}\right)dt}d\Omega }={\displaystyle {\int }_{V_h}{\left[\varphi {\displaystyle \sum _j\left(S_j{\rho }_j{X}_{ij}\right)}\right]}^{t+\Delta t}-{\left[\varphi {\displaystyle \sum _j\left(S_j{\rho }_j{X}_{ij}\right)}\right]}^t}d\Omega \\ =V_h{\left[\varphi {\displaystyle \sum _j\left(S_j{\rho }_j{X}_{ij}\right)}\right]}^{t+\Delta t}-{\left[\varphi {\displaystyle \sum _j\left(S_j{\rho }_j{X}_{ij}\right)}\right]}^t\end{array}$$

(10)

Then the integration is furtherly approximated as:

$${\displaystyle {\int }_t^{t+\Delta t}{\displaystyle {\int }_{V_h}{\displaystyle \sum _l\left({{\nu }^{\prime }}_{i,l}-{\nu }_{i,l}\right)r_l}d\Omega dt}}={\displaystyle {\int }_t^{t+\Delta t}{V}_h{\displaystyle \sum _l\left({{\nu }^{\prime }}_{i,l}-{\nu }_{i,l}\right)r_l}dt}$$

(11)

$${\displaystyle {\int }_t^{t+\Delta t}{\displaystyle {\int }_{V_h}{\displaystyle \sum _j\left({\rho }_j{X}_{ij}{q}_j{}^W\right)}d\Omega }dt}={\displaystyle {\int }_t^{t+\Delta t}{\displaystyle \sum _j\left({\rho }_j{X}_{ij}{Q}_j{}^W\right)}dt}$$

(12)

where $Q_j{}^W={\displaystyle {\int }_{V_h}{\displaystyle \sum _j{q}_j{}^W}d\Omega }$.

The second term at the left side of Equation (9) is estimated by a Gaussian formula, and the following is obtained:

$${\displaystyle {\int }_t^{t+\Delta t}{\displaystyle {\int }_{V_h}\nabla \cdot {\displaystyle \sum _j\left({\rho }_j{X}_{ij}{u}_j\right)}d\Omega }}dt={\displaystyle {\int }_t^{t+\Delta t}{\displaystyle {\int }_{\partial V_h}{\displaystyle \sum _j\left({\rho }_j{X}_{ij}{u}_j\right)}\cdot \stackrel{\rightharpoonup }{n}d\eta }}dt$$

(13)

The two-point linear approximation scheme of interface flux is used to further estimate the right side of Equation (13) as:

$${\displaystyle {\int }_{\partial V_h}{\displaystyle \sum _j\left({\rho }_j{X}_{ij}{u}_j\right)}\cdot \stackrel{\rightharpoonup }{n}d\eta }=-{\displaystyle \sum _k^{n_h}{\displaystyle \sum _j\left({\rho }_{j,hk}{X}_{ij,hk}{T^{\prime }}_{j,hk}\left(p_{j,h}-p_{j,k}\right)\right)}}$$

(14)

where ${T^{\prime }}_{j,hk}$ is the transmissibility of phase j between grid m and grid k, which is equal to the product of j-th phase mobility and geometric factor, that is:

$${T^{\prime }}_{j,hk}={\lambda }_{j,hk}{G}_{hk}$$

(15)

For two adjacent rectangular grids, the geometric factor is calculated as:

$$G_{hk}=\frac{1}{2}{\left(G_h{}^{-1}+G_k{}^{-1}\right)}^{-1} G_h=\frac{k_h{A}_{hk}}{d_h} G_k=\frac{k_k{A}_{hk}}{d_k}$$

(16)

where $k_h$ and $k_k$ are the permeability of grid h and grid k respectively, and $A_{mk}$ are the area of the interface of grid h and grid k.

The mobility ${\lambda }_{j,hk}$ is calculated as:

$${\lambda }_{j,hk}=\frac{k_{rj,hk}}{{\mu }_{j,hk}{B}_{j,hk}}$$

(17)

where the terms related to saturation (relative permeability) adopt an upstream weight scheme, and the terms related to pressure adopts an arithmetic average scheme, which are:

$$k_{rj,hk}=\left\{\begin{array}{c}{k}_{rj,h}\left(p_{j,m}-{\rho }_{j,sc}g\frac{D_h}{B_{j,hk}}\right)\ge \left(p_{w,j}-{\rho }_{w,sc}g\frac{D_k}{B_{j,hk}}\right)\\ k_{rj,h}\left(p_{j,m}-{\rho }_{j,sc}g\frac{D_h}{B_{j,hk}}\right)<\left(p_{w,j}-{\rho }_{w,sc}g\frac{D_k}{B_{j,hk}}\right)\end{array}\right.$$

(18)

$${\mu }_{j,hk}=\frac{{\mu }_{j,h}+{\mu }_{j,k}}{2}$$

(19)

$$B_{j,hk}=\frac{B_{j,h}+B_{j,k}}{2}$$

(20)

Then, Equation (14) is rewritten as:

$${\displaystyle {\int }_t^{t+\Delta t}{\displaystyle {\int }_{V_h}\nabla \cdot {\displaystyle \sum _j\left({\rho }_j{X}_{ij}{u}_j\right)}d\Omega }}dt={\displaystyle {\int }_t^{t+\Delta t}-{\displaystyle \sum _k^{n_h}{\displaystyle \sum _j\left({\rho }_{j,hk}{X}_{ij,hk}{T^{\prime }}_{j,hk}\left(p_{j,h}-p_{j,k}\right)\right)}}}dt$$

(21)

The implicit scheme is used and the following are obtained:

$${\displaystyle {\int }_t^{t+\Delta t}{\displaystyle {\int }_{V_h}\nabla \cdot {\displaystyle \sum _j\left({\rho }_j{X}_{ij}{u}_j\right)}d\Omega }}dt=\Delta t{\left[-{\displaystyle \sum _k^{n_h}{\displaystyle \sum _j\left({\rho }_{j,hk}{X}_{ij,hk}{T^{\prime }}_{j,hk}\left(p_{j,h}-p_{j,k}\right)\right)}}\right]}^{t+\Delta t}$$

(22)

$${\displaystyle {\int }_t^{t+\Delta t}{\displaystyle {\int }_{V_h}{\displaystyle \sum _l\left({{\nu }^{\prime }}_{i,l}-{\nu }_{i,l}\right)r_l}d\Omega dt}}={\displaystyle {\int }_t^{t+\Delta t}{V}_h{\displaystyle \sum _l\left({{\nu }^{\prime }}_{i,l}-{\nu }_{i,l}\right)r_l}dt}=\Delta t{V}_h{\left[{\displaystyle \sum _l\left({{\nu }^{\prime }}_{i,l}-{\nu }_{i,l}\right)r_l}\right]}^{t+\Delta t}$$

(23)

$${\displaystyle {\int }_t^{t+\Delta t}{\displaystyle {\int }_{V_h}{\displaystyle \sum _j\left({\rho }_j{X}_{ij}{q}_j{}^W\right)}d\Omega }dt}={\displaystyle {\int }_t^{t+\Delta t}{\displaystyle \sum _j\left({\rho }_j{X}_{ij}{Q}_j{}^W\right)}dt}=\Delta t{\left[{\displaystyle \sum _j\left({\rho }_j{X}_{ij}{Q}_j{}^W\right)}\right]}^{t+\Delta t}$$

(24)

Then combining Equations (22)–(24), it is obtained that:

$$\begin{array}{l}{\left[-{\displaystyle \sum _k^{n_h}{\displaystyle \sum _j\left({\rho }_{j,hk}{X}_{ij,hk}{T^{\prime }}_{j,hk}\left(p_{j,h}-p_{j,k}\right)\right)}}\right]}^{t+\Delta t}+V_h{\left[{\displaystyle \sum _l\left({{\nu }^{\prime }}_{i,l}-{\nu }_{i,l}\right)r_l}\right]}^{t+\Delta t}\\ +{\left[{\displaystyle \sum _j\left({\rho }_j{X}_{ij}{Q}_j{}^W\right)}\right]}^{t+\Delta t}=\frac{V_h}{\Delta t}{\left[\varphi {\displaystyle \sum _j\left(S_j{\rho }_j{X}_{ij}\right)}\right]}^{t+\Delta t}-{\left[\varphi {\displaystyle \sum _j\left(S_j{\rho }_j{X}_{ij}\right)}\right]}^t\end{array}$$

(25)

Similarly, for kerogen in the solid phase, the discrete scheme of Equation (2) can be obtained as:

$$V_h{\left[{\displaystyle \sum _k\left({{\nu }^{\prime }}_{i,k}-{\nu }_{i,k}\right)r_k}\right]}^{t+\Delta t}=\frac{V_h}{\Delta t}\left(C_{kero}{}^{t+\Delta t}-C_{kero}{}^t\right)$$

(26)

For coke in the solid phase, the discrete scheme of Equation (3) is:

$$V_h{\left[{\displaystyle \sum _i{\displaystyle \sum _k\left({\nu }_{i,k}-{{\nu }^{\prime }}_{i,k}\right)r_k}}\right]}^{t+\Delta t}=\frac{V_h}{\Delta t}\left(C_C{}^{t+\Delta t}-C_C{}^t\right)$$

(27)

The discrete scheme for energy conservation equation can be obtained as:

$$\begin{array}{l}{\left[-{\displaystyle \sum _k^{n_h}{\lambda }_{hk}\left(T_h-T_k\right)}\right]}^{t+\Delta t}+{\left[-{\displaystyle \sum _k^{n_h}{\displaystyle \sum _j\left({\rho }_{j,hk}{H}_{j,hk}{T^{\prime }}_{j,hk}\left(p_{j,h}-p_{j,k}\right)\right)}}\right]}^{t+\Delta t}+V_h{\left({\displaystyle \sum _l{H}_l{r}_l}\right)}^{t+\Delta t}\\ +V_h{\left({\displaystyle \sum _j\left(H_j{\rho }_j{q}_j\right)}\right)}^{t+\Delta t}+Q_H{}^{W,t+\Delta t}\\ =\frac{V_h}{\Delta t}\left[{\left(\varphi {\displaystyle \sum _j{S}_j{\rho }_j{U}_j}+\left(1-\varphi \right)U_r\right)}^{t+\Delta t}-{\left(\varphi {\displaystyle \sum _j{S}_j{\rho }_j{U}_j}+\left(1-\varphi \right)U_r\right)}^t\right]\end{array}$$

(28)

where, $H_{j,hk}$ is the enthalpy of the j-th phase between the h-th and k-th grids, which adopts the upstream weight scheme, i.e.,

$$H_{j,hk}=\left\{\begin{array}{c}{H}_{j,h}\left(p_{j,h}-{\rho }_{j,sc}g\frac{D_h}{B_{j,hk}}\right)\ge \left(p_{w,j}-{\rho }_{w,sc}g\frac{D_k}{B_{j,hk}}\right)\\ H_{j,k}\left(p_{j,h}-{\rho }_{j,sc}g\frac{D_h}{B_{j,hk}}\right)<\left(p_{w,j}-{\rho }_{w,sc}g\frac{D_k}{B_{j,hk}}\right)\end{array}\right.$$

(29)

2.3. Treatment of Hydraulic Fracture

Because the permeability of oil shale is very low, it is difficult to directly inject hot steam into oil shale reservoir, so the method of fracturing a vertical well is often used to inject hot steam. When fracturing vertical wells are used, the injected hot steam will flow along the high-conductivity fractures and then into the reservoir. Therefore, the numerical modeling of the flow in fractures is very important for the numerical simulation of steam injection assisted in situ production of oil shale. There have been lots of treatment methods for various fractures including hydraulic fractures, natural fractures, and induced fractures in the reservoir [21,22,23,24]. In this paper, due to the embedded discrete fracture model (EDFM) being easy to couple in the calculation of the thermal compositional model in the case of Cartesian mesh, the commonly used EDFM is used to explicitly handle the large-scale discrete fractures of fractured vertical wells [12,13,14]. The main idea is to embed discrete fractures into the Cartesian matrix grid, in which the transmissibility of phase j between matrix grid m and fracture grid f is:

$${T^{\prime }}_{j,mf}={\lambda }_{j,mf}{G}_{mf}$$

(30)

where λ_j,mf is the mobility of the jth phase, and the upwind scheme is used. G_mf is:

$$G_{mf}=k_{mf}\frac{A_f}{〈d〉}, 〈d〉=\frac{{\displaystyle {\int }_{{\Omega }_m}\left|\stackrel{\rightharpoonup }{x}\cdot \stackrel{\rightharpoonup }{n}\right|d\Omega }}{V_m}$$

(31)

where $k_{mf}$ is the harmonic average value of fracture grid permeability and matrix grid permeability, and $〈d〉$ is the average distance from points in the matrix grid to the fracture grid.

2.4. Non-Linear Solver

The automatic differentiation-based Newton iteration method is generally used to solve the above non-linear equations. However, in the simulation of oil-shale in situ production, the initial condition of oil saturation is:

$$S_o\left(t=0\right)=0$$

(32)

Therefore, it is necessary to judge whether there will be oil phase in each grid during the calculation process. If there is no oil phase, there is only a gas phase for each light and heavy component generated by oil shale pyrolysis, and the above phase equilibrium equation does not hold. If phase change occurs in a grid and oil phase is produced, the phase equilibrium equation between oil phase and gas phase just holds. Therefore, for different grids, generally, due to the different phase state, the equations within each grid may be different, which poses challenges to the smoothness of the non-linear calculation.

In order to enhance the smoothness of the Newton iteration-based solution process and improve the calculation efficiency without significantly reducing the calculation accuracy, the initial condition of oil saturation is set as:

$$S_o\left(t=0\right)=0.0001$$

(33)

To make each grid initially contain oil, the equations at the initial time of each grid are the same without judging in which grid the phase change occurs. Because there will be a residual saturation term in the relative permeability curve, that is, when the oil saturation in the grid is lower than that, the oil phase in the grid will not flow out if residual oil saturation meets:

$$S_{or}>S_o\left(t=0\right)=0.0001$$

(34)

It can be ensured that each grid will have at least oil saturation $S_o\left(t=0\right)=0.0001$ in the calculation process, so that there is no need to judge the phase change in the whole calculation process, and in each grid the control equations are the same, which improves the smoothness of the non-linear solution by Newton iteration.

In order to make the obtained solution by Newton iteration physical, some limiting equations need to be added to limit the concentration of each component in the range of 0–1, that is:

$$0\le X_{ij}\le 1$$

(35)

In this paper, the inequality constraint is transformed into an equality constraint and added to the global equations by introducing relaxation variables ${\nu }_{ij}$, that is,

$$X_{ij}-{\sin }^2\left({\nu }_{ij}\right)=0$$

(36)

Therefore, based on the block-center finite volume method that satisfies local mass conservation and using the two-point linear approximation scheme of the orthogonal grid, the 17 discrete equations are obtained, and there are also 17 independent variables, including: pressure, water saturation, gas saturation, solid-phase saturation, kerogen concentration in solid phase, temperature, concentration of five fluid components in gas phase except water components, and relaxation variables constrained by inequality.

2.5. Treatment of Boundary Condition

The reservoir boundary condition is closed, but it is still thermally conductive. In the commercial numerical simulator CMG, the boundary is closed, that is, there is neither mass transfer nor heat conduction. As shown in Figure 1a, the simulator achieves the purpose of little mass transfer but normal heat conduction by adding a circle of grid with extremely low permeability around the calculation domain, so as to form a more real reservoir boundary condition. In the pressure profile shown in Figure 1b, since the permeability of the grid added outside the calculation domain is very low, and the pressure of the grid added during the simulation calculation is always the initial reservoir pressure.

3. Numerical Examples

3.1. Validation of the Smooth Non-Linear Solver: A 2D Case with a Heater Device

As shown in Figure 2, this example is a two-dimensional reservoir model, and the grid size is 5 m × 5 m, the grid division is 11 × 11. The reservoir thickness is 5 m and the initial pressure is 10 MPa. The reservoir boundary is closed both in terms of the mass and heat transfer. The heat injection grids are (4,4), (4,8), (8,4), (8,4), the constant temperature of the heating device is 370 °C, the grid where the production well is located is (6,6), and the production pressure is 5 MPa. The chemical reaction in this example is shown in Equation (37), and the corresponding chemical reaction coefficient is shown in

Table 1. Chemical reaction coefficient of kerogen pyrolysis.

Chemical Reaction Coefficient						Frequency Factor	The Activation Energy (KJ/mol)
	kerogen	IC2	IC13	IC37	CO2
reaction	1	0	0	0	0	3.74 × 1012	161.600
production	0	0.04475	0.0178	0.0096	0.00541

, The enthalpy of each component is calculated according to the fifth degree polynomial function about temperature [6]. Other basic fluids, thermal, and reservoir parameters are shown in

Table 2. Value physical properties.

Property	Value	Property	Value
Initial water saturation, fraction	0	Rock thermal conductivity, J/(m·s·K)	3
Initial gas saturation, fraction	0.1499	Solid-phase thermal conductivity, J/(m·s·K)	3
Initial oil saturation, fraction	0.0001	Oil phase thermal conductivity, J/(m·s·K)	0.6
Initial solid saturation, fraction	0.8500	Water phase thermal conductivity, J/(m·s·K)	0.2
Kerogen concentration in a solid phase, fraction	1	Gas-phase thermal conductivity, J/(m·s·K)	0.1
Initial reservoir temperature, °C	20 °C	Initial oil viscosity, mPa·s	50
Initial porosity, fraction	0.3	Water viscosity, mPa·s	0.6
Permeability, mD	1	Gas viscosity, mPa·s	0.01
Rock density, kg/m³	2700	Rock compressibility, MPa−1	1.07 × 10−4
Solid density, kg/m³	2000	Solid-phase compressibility, MPa−1	1.07 × 10−4
Oil density, kg/m³	877	Oil phase compressibility, MPa−1	3.02 × 10−3
Gas density, kg/m³	26	Water phase compressibility, MPa−1	5 × 10−4
Water density, kg/m³	1000	Gas phase compressibility, MPa−1	5 × 10−2

. The coefficients of the polynomial function can be seen in Fan et al. [6]. Figure 3 and Figure 4, respectively, show pressure profiles at 0.001, 10, 30, and 200 days calculated by this method and CMG. Figure 5 and Figure 6, respectively, show the oil saturation profiles at 0.001, 10, 30, and 200 days calculated by this method and CMG. Figure 7 and Figure 8, respectively, compare the oil saturation profiles at different times calculated by this method. Figure 9 compares the production rates of components for 500 days calculated by this method and CMG, respectively. Comparisons of calculation results from this method and CMG illustrate that this method can obtain high-accuracy results.

In addition, it can also be seen that, at the initial stage of production, the pyrolysis reaction rate of solid kerogen-producing liquid and gaseous organic matter is very low, resulting in the pressure drop of depleted exploitation of the saturated inorganic CO₂ and N₂ in the reservoir pores, shown in Figure 3a at the beginning of production; and then Figure 3b shows the depleted pressure wave diffuses to the boundary. The pressure of the whole reservoir calculation domain is very close to the bottom hole flow pressure (BHP). Then, with the increase of the surrounding temperature of the heating device shown in Figure 5, the pyrolysis reaction rate of solid kerogen increases. The IC₂, IC₁₃, IC_37, and CO₂ produced by pyrolysis increase the pressure near the heating device, that is, the pressure distribution shown in Figure 3c. Finally, as shown in Figure 3d, the organic and inorganic substances produced by pyrolysis reaction are recovered, and the pressure of the whole reservoir calculation domain is reduced to close to the BHP again. Besides, with the increase of pyrolysis reaction rate of oil shale, the rate of carbon dioxide produced by pyrolysis is also faster, resulting in the output of CO₂ in the later stage of production being higher than that of N₂.

$$\mathrm{KER}\to 0.0096{\mathrm{IC}}_{37}+0.0178{\mathrm{IC}}_{13}+0.04475{\mathrm{IC}}_2+0.00541{\mathrm{CO}}_2$$

(37)

3.2. A 3D Example with a Steam Injection Well

As shown in Figure 10, this example is an oil shale reservoir with a size of 30 m × 30 m * 18 m. It is divided into two layers longitudinally, and a circle of heat transfer layer is added outside the boundary. The initial conditions of this reservoir model are listed in

Table 3. Initial conditions.

Initial Conditions	Values
Pressure, MPa	10
Water saturation, fraction	0
Gas saturation, fraction	0.15
Temperature, fraction	20
Solid-phase saturation, fraction	0.85
N2 concentration in gas phase, fraction	0.2
CO2 concentration in the gas phase, fraction	0.8
Kerogen concentration in the solid phase, fraction	1

, and related physical and thermal parameters are the same as those in Table 1. The reservoir boundary is closed but thermally open, and the pyrolysis reaction of oil shale is the same as that in Equation (37). There is a fractured vertical well for steam injection, the injection rate is 25 tons/day, the steam temperature is 480 °C, and a production well is controlled by a fixed BHP of 5 MPa. The pressure distribution in Figure 10c can also reflect the well locations, and the simulation time lasts for 800 days, Figure 11 shows the calculated temperature and water saturation at different times. It can be seen that the water saturation distribution is positively correlated with the temperature distribution, and the injected steam is effectively heating the formation. Figure 12 and Figure 13 also show the calculated oil saturation and IC2 light component in the gas phase, respectively. Figure 12 shows that the oil saturation increases with the continuous heat injection, and the corresponding concentration of IC2 component in the gas phase in Figure 13 also increases, indicating that a large number of IC2 light components are produced by oil shale pyrolysis, and the profiles of the concentration of IC2 component in the gas phase is basically consistent with the profiles of oil saturation, which represents the distribution of the intensity of the oil shale pyrolysis reaction. Because the fracture only penetrates the lower layer, the lower layer is heated more, the reaction is stronger, and the difference between layers is obvious.

Figure 14 shows the calculated production rates of the components for 800 days. The production curves of inorganic and organic components of the well can also be obtained. It can be seen that, during early production, the production rate of saturated N₂ and CO₂ declines sharply. In less than a day, almost all the inorganic gases in the formation have been extracted. Moreover, since the chemical reaction coefficient of the light component IC2 in the pyrolysis reaction is the largest, and the proportion of the light component IC2 in the gas phase at phase equilibrium is much higher than that of the heavy components IC13 and IC37, the light component IC2 will be recovered first with the increase of the chemical reaction speed caused by the continuous temperature increase with steam injection. The output of the organic components is low in this small-scale toy model, so we can continue to simulate and calculate the production performance for a longer time, or improve the steam injection rate to produce organic components faster.

4. Conclusions and Future Work

This paper presents an effective numerical simulation method for oil-shale in situ production. Firstly, the finite volume discrete scheme of mass and heat transfer equation for oil shale in situ production is derived, and the hydraulic fracture of the steam injection well is treated by using the embedded discrete fracture model. Then, a non-linear solver which does not need to judge phase transformation and can give unified constraints on phase saturation and component concentration is given, Thus, the mass and heat transfer equation of oil shale can be solved more smoothly. A numerical example compared with CMG shows that the method based on the smooth non-linear solver can achieve accurate results. Through a numerical example with a fractured vertical well and open thermal boundary conditions, it is verified that this method can effectively deal with the actual situation of steam injection-assisted oil shale production which is difficult to deal with by the previous methods.

There are several limitations in this paper: (1) although the presented smooth nonlinear solver ensures that the equation form of each grid is always the same in the process of simulation calculation, additional constraint equations shown in Equation (20) are added, and the more components, the more additional equations are added, significantly increasing the computational cost. Therefore, it will be valuable in future work to strictly analyze and compare the calculation efficiency and accuracy between the presented non-linear solver and the widely used non-linear solver of the thermal composition model for oil-shale in situ recovery based on natural formulations. (2) In the processing method of the thermally open condition given in this paper, if the grid size added outside the boundary of the computing domain is too large, the approximation accuracy of the thermally open condition will be reduced. If the grid size is small, more grids need to be added, which will reduce the computing efficiency. How to determine the optimal strategy of adding grids or develop a better thermally open condition processing method will be valuable to ascertain in future work. (3) In this paper, only some small-scale conceptual examples are calculated and tested. The application of this method in the actual large-scale oil-shale development can be further studied. (4) The application of the presented smooth non-linear solver in the steam-flooding reservoir and the thermochemical-flooding reservoir can be further explored.