Three-Dimensional Numerical Simulation of Hydrocarbon Production and Reservoir Deformation of Oil Shale In Situ Conversion Processing Using a Downhole Burner.

A three-dimensional numerical simulation of oil shale in situ conversion processing by applying the downhole burner heating technology was conducted. The evolution of the fluid vector and temperature field and the characteristic of kerogen decomposition and oil and gas production were analyzed. The effects of different burning temperatures and gas injection velocities on the thermal evolution processing of oil shale in situ conversion were investigated. The stress-strain and deformation of the oil shale stratum during in situ processing were studied. The results show that kerogen decomposition is a thermo-kinetically controlled mechanism. Both the gas injection velocity and burning temperature can enhance the kerogen decomposition and oil production, especially for the latter one. In addition, the stratum-deformation of oil shale should be considered for oil shale in situ conversion processing, especially for the long-term operational lifetime.

Introduction

With continuous and increasing consumption of conventional petrochemical energy resources, the world is facing severe challenges to find alternative or unconventional energy resources to replace oil or ease the tightness of the oil supply. Oil shale, an unconventional fossil energy resource, contains kerogen and can produce petroleum-like hydrocarbons[1−3] and has become a promising alternative energy source due to its vast in-place reserve.[4−6] According to the U.S. Geological Survey estimates, the world’s oil shale reserves are about 2.8–3.0 × 1012 barrels of shale oil,[7] which is 2.5 times the world’s conventional crude oil reserves. Generally, there are two ways to commercialize this resource, namely, generating electrical energy via combustion and recovering petrochemical products via pyrolysis. The latter, recovering petrochemical products from oil shale by surface retorting or by in situ conversion, is undoubtedly a promising way to maximize the value of this rock. As the primary and widely used method for recovering petrochemical products, the surface retorting processes, such as the ATP, Petrosix, Kiviter, and Fushun-type retorts, include mining, crushing, retorting, and upgrading and thus make it more costly than the production of conventional crude oil financially and in terms of its environmental impact.[8] However, the in situ processes, such as the Shell in situ conversion process (Shell ICP), Electrofrac, and Chevron, have directly processed the heating of the shale underground mainly by thermal conduction and convection, and therefore, their methods have gained tremendous attention and have been considered to be the most attractive practical methods for recovering hydrocarbon from oil shale due to their better applicability to deep layers and related minor environmental issues.[5,9−16] However, a major problem in performing oil shale in situ processing is how to provide sufficient heat efficiently to the underground for kerogen decomposition, which takes place at a high temperature (at least 300 °C).

Since the 1960s, numerous studies have been conducted on the experiments and simulations of oil shale in situ conversion processing, and diverse technologies have been put forward using various heat sources and heat delivery systems, attempting to find economically and environmentally acceptable methods of commercializing such resources. Among them, the Occidental modified in situ (MIS) process, inducing combustion air to permeate and burn oil shale underground after mining 20% of the oil shale and fracturing the rest to create a void space of 20–25%, had been conducted in eight in situ field pilots of different scales in Colorado by 1983, with the improved shale oil and high oil yield exceeding 60% of the Fisher assay, successfully demonstrating the technology. The MIS process was ensured of a viable technology at the time; however, Occidental company did not carry out further research due to the difficulty of combustion control and pollution to the underground aquifer due to the damage of the barrier by mining as well as extensive drilling engineering and problematic excavation. Shell ICP, currently the most mature oil shale in situ technology, heats oil shale by thermal conduction from a closely spaced array of electric resistance heaters. The ICP has been conducted by eight Colorado and one Jordan field pilots in the oil shale resource with increasing scope and complexity, which costs vast amounts of money.[11,12] The disadvantages, such as complicated process, high cost, and low oil and gas migration power, make the economics of ICP technology not warrant commercial development.[9] The conduction, convection, and reflux processes from American Shale Oil company, characterized by a faster heating rate than Shell ICP, use refluxing oil to transport heat provided from a horizontal heater in an “L”-shaped well using a downhole burner to the retort boundary.[15] Also, this technology was initially planned to take 10 years, from 2008 to 2018, to complete the pilot test and switch to commercial operations. Still, the project was announced to be discontinued in 2016 due to a technical problem with the downhole burner. Recently, a downhole burner heating technology (DBHT), hot gas convection, and reflux method using burners placed in the borehole to heat refluxing gases were devised by the Jilin Zhongcheng Oil Shale Group Co,. Ltd (.China) to recover hydrocarbon products from the oil shale layer.[17] This technology uses hydraulic fracturing to create channels for oil and gas transportation and re-injects the combustible gas produced during the oil shale in situ processing. The field pilot of the DBHT process was conducted in the Songliao Basin, China, and the first barrel of shale oil recovered from underground in China was produced in July 2014, proving the feasibility and adaptability of the DBHT process.

Although these technologies have achieved varying degrees of success, none has been put into commercial operations. A range of studies and experiments are needed to understand the impact of various factors on oil shale in situ processing. Compared to the field tests and pilots, numerical simulation is a low-cost and helpful tool for designing and controlling operations and identifying and interpreting observations that differ from model predictions and has been used widely to explore the oil shale in situ conversion processing.[6,18−34] For instance, Fan et al.[27] investigated the effects of temperature and well layout on the oil shale in situ pyrolysis by numerically simulating the electrical heating process and showed that the oil production rate was highly dependent on electrical heating temperature. Youtsos et al.[28] built a one-dimensional model to analyze the thermal front development in porous reservoirs with reaction flows. They found that hot gas injection was more effective in heating the stratum. Guo et al.[32] numerically simulated and evaluated the oil shale in situ conversion and found that downhole heating is a promising technology for high oil productivity. However, most previous studies were based on the one- or two-dimensional model and focused on the effect of technology, well layout plan, and experimental parameters on the temperature field and oil and gas yields of oil shale in situ processing.[6,33−37] With the decomposition of kerogen and the increase of the pores in oil shale, the stratum strength decreases, which may cause the deformation of the stratum and consequently have an adverse effect on the well pipe and the regular operation of underground equipment. However, to the best of our knowledge, no literature is available to date regarding the work on the stress–strain and deformation of the oil shale stratum after in situ conversion.

In this study, a three-dimensional (3D) model based on the field pilot of DBHT in the Songliao Basin was established, and the numerical simulation of oil shale in situ pyrolysis by applying the idea of DBHT was conducted to examine its temperature field, thermal processes, stress–strain, and deformation. The effects of different burning temperatures (800, 850, 900, 950, and 1000 K) and different gas injection velocities (2, 5, 10, 15, and 20 m/s) on the thermal evolution processing were investigated. The stress–strain and deformation of the oil shale stratum during in situ processing of various cases were analyzed for the first time. We anticipate that this study can provide comprehensive insights and suggestions for applying the DBHT process for oil shale in situ conversion, which will be essential and helpful for designing and optimizing the next field scale oil shale in situ project.

Model Description and Parameters

The model set up in the paper is based on the field pilot of DBHT in the Songliao Basin. The Songliao Basin hosts a remarkably high number of oil shale layers,[38−40] but their oil shale is lean (low oil yield and high ash) and deeply buried, making it a vital testing site for oil shale in situ exploitation. The oil shale explored in the field pilot of DBHT is the Qingshankou Formation, with a burial depth above 300 m, an average oil yield of 5.3%, and an oil shale layer thickness of 4–4.5 m. The upper and lower oil shale layers are mudstone layers, an excellent natural water barrier. The oil shale layer of 4 m was considered in this study, and the thickness of the upper and lower mudstones in the geometric model is set to 2 m.

Geometric Model and Meshing

The geometric model of the simulation calculation is shown in Figure . As the entire model is mirror-symmetrical, only half of it is taken as the simulation calculation area. The geometric model of the simulation calculation area is a rectangular parallelepiped with a length of 20 m in the x-direction, a height of 8 m in the y-direction, and a width of 10 m in the z-direction. Moreover, the size of the fracture area is 10 × 5 × 5 mm. The gas injection well and the production well with a distance of 5 m are connected through the fracture area. A fixed heat source is set in the gas injection well with 50 cm above the fracture area to simulate the flame produced by the burner.

Figure 1

Geometric and meshing diagram of the simulation calculation model.

The model is meshed by a pure hexahedral (Figure ), for there is a fracture area with an aspect ratio of 1:2000. The entire model consists of fluid calculation areas, including gas injection well, production well and fracture, and solid calculation area, including the oil shale layer and the top and bottom mudstones. A boundary layer is set at the fluid calculation areas in the fluid–solid coupling interface. Also, the fluid calculation areas are separately refined in quantification during meshing to meet the calculation requirements of the solver FLUENT.

Selected Parameters and Boundary Conditions

The physical parameters of injected gas, oil shale, and mudstone based on experimental data from the field pilot of DBHT in the Songliao Basin used in this paper are shown in Table . In the simulation process, the oil shale was heated by hot gas heated by the heat exchange with a burner. The far boundary temperature was 300 K.

Table 1

Physical Parameters of the Injected Gas, Oil Shale, and Mudstone

property	value
injected gas	density (kg·m–3)	–1 × 10–8T3 + 3 × 10–5T2– 0.0326T + 13.192 (R2 = 0.9909)
	specific heat capacity (J·kg–1·K–1)	0.1278T + 976.1 (R2 = 0.9958)
	thermal conductivity (W·m–1·K–1)	–2 × 10–8T2+ 8 × 10–5T + 0.0034 (R2 = 0.9989)
oil shale	initial porosity	2%[41]
	initial density (kg·m–3)	2050
	specific heat capacity (J·kg–1·K–1)	0.0019T2 – 4.7631T + 3548.8 (R2 = 0.8727)
	thermal conductivity parallel to the shale bedding (W·m–1·K–1)	–0.0009T + 0.9316 (R2 = 0.8067)
	thermal conductivity of the shale perpendicular to the bedding (W·m–1·K–1)	–0.0006T + 0.7227 (R2 = 0.7766)
mudstone	density (kg·m–3)	1900
	specific heat capacity (J·kg–1·K–1)	1000
	thermal conductivity (W·m–1·K–1)	1.010

Several cases are studied to examine the effects of different burning temperatures (800, 850, 900, 950, and 1000 K) and different gas injection velocities (2, 5, 10, 15, and 20 m/s) on the thermal processes of oil shale in situ pyrolysis. Both the steady values and the transient values of 10 years are calculated. For the convenience of discussion, each case is denoted as M–x–y, where x is the gas injection velocity and y is the burner’s temperature. For example, M-10-900 represents the case in which the gas injection velocity is 10 m/s, and the burning temperature of the burner is 900 K.

Mathematical Model

Energy Conservation Equation

The temperature field of oil shale in situ conversion follows the first law of thermodynamics. The expression of energy conservation is in ref (42)where Q denotes the quantity of energy supplied to the closed system as heat; W denotes the amount of thermodynamic work done by the closed system on its surroundings; ΔU is the change of the internal energy of the closed system; ΔEk is the change of the kinetic energy of the closed system; ΔEp is the change of the potential energy of the closed system.

As the heat transfer characteristics of the oil shale in situ conversion follow , Formula can be simplified as

The heat transfer phenomenon here is heat conduction, which complies with Fourier’s law of heat conduction. The heat transfer Q in time t iswhere t is the time; Kc is the thermal conductivity; T is the temperature; A is the heat transfer area; and d is the distance between the heat transfer surfaces.

According to the energy conservation equation and Fourier’s law, the expression of the thermal conductivity differential equation is as follows[43]where ρ is the fluid density; c is the specific heat capacity; λ is the thermal conductivity coefficient; and Φ is the heat generated by the internal heat source in the micro-element body per unit time.

For the transient heat conduction process studied in this paper, the definite solution includes the initial condition, the temperature distribution of the stratum temperature field at the initial moment, and the boundary condition—the heat transfer on the stratum boundary, using the Dirichlet boundary condition, which is

Fluid Motion Control Equation

The fluid motion in this model follows the conservation of mass, conservation of momentum, and conservation of energy[44,45] as follows

Continuity equation (conservation of mass)where u, v, and w are the fluid velocities in the x-direction, y-direction, and z-direction, respectively.

Equation of motion (conservation of momentum)where P is the static pressure on the fluid, τ is the viscous stress on the fluid surface, and V is the velocity vector of the fluid.

Conservation of energywhere e is the internal energy of the fluid and q is the heat energy received by the fluid.

Grid Independence Verification and Model Validation

To ensure the reliability of the simulation results, a grid independent verification is conducted for the present numerical model with the gas injection velocity of 10 m/s and the burning temperature of 900 K. Three grid systems of Grid 1, Grid 2, and Grid 3 are selected for simulations with the mesh numbers of 822733, 1810013, and 3258024, respectively. The temperature evolution results of the point with coordinates (25, 5, and 20) in the grid are shown in Figure S1. The maximum error between Grid 1 and Grid 2 is about 1.62%, which is reduced to 0.47% from Grid 2 to Grid 3. Thus, Gird 3 with the mesh number of 3258024 was chosen for simulation as a tradeoff between calculation time and model accuracy.

The model used in the present study was established based on the real physical parameters of injected gas, oil shale, and mudstone as well as the experimental data from the field pilot of DBHT in the Songliao Basin and a small lab experimental system.

Evaluation Parameters

Function Fitting of the Oil Production Rate of Oil Shale

Figure shows the kerogen convertible and oil production curves obtained by Shell ICP pilot data as well as the fitting function of these two curves, which are used for fitting and analyzing the kerogen conversion and oil production rate of various cases in this study. Moreover, several critical temperatures are defined according to the oil production rate function of oil shale to compare the heating efficiency of different models.

Figure 2

Kerogen convertible and oil production. Reprinted in part with permission from ref (12). Copyright 2010 American Chemical Society.

Excessively high temperature: the temperature over which no more oil will be produced. The lowest excessively high temperature is 688.5 K, which means that there is no oil produced when the temperature of oil shale is over 688.5 K.

Invalid temperature: the temperature below which no oil is yielded. The highest invalid temperature is 573.5 K, which means the oil shale starts to produce oil only when the oil shale is heated at more than 573.5 K.

Ideal oil production temperature: the temperature at which the oil production rate of oil shale is more than 50% and the kerogen content is less than 30%. According to Figure , the ideal oil production temperature is between 643.5 and 688.5 K.

Evaluation Parameters of In Situ Pyrolysis of Oil Shale

Based on the kerogen convertible, oil production, and the critical temperatures defined above, the following evaluation parameters are defined and used to analyze the oil production efficiency of various cases.

(1) The mass-weighted average temperature (T̅, K), the temperature of each discrete unit of oil shale times its weight divided by the total weight of oil shalewhere ρos is the density of oil shale (kg/m3), V is the volume of the ith discrete unit of oil shale (m3), and T is the temperature of the ith discrete unit of oil shale (K).

(2) The ideal oil production temperature volume (Vp, m3), the total volume of all discrete units of oil shale whose temperature exceeds 643.5 kwhere VT is the volume of the discrete unit at temperature T (m3).

(3) The excessively high-temperature volume (Vh, m3), the total volume of all discrete units of oil shale whose temperature exceeds 688.5 K

(4) The invalid temperature volume (Vn, m3), the total volume of all discrete units of oil shale whose temperature is lower than 573.5 K

(5) The Kerogen content of oil shale (Ke, %), the ratio of kerogen content of oil shale at a certain temperature to the total kerogen before pyrolysis.where K(T) is the Kerogen content of the ith discrete unit of the oil shale at the temperature of T (%).

(6) The actual oil production rate of oil shale (W, %), the mass ratio of the shale oil produced at a certain temperature to the total shale oil after complete pyrolysiswhere λos is the oil yield of oil shale (%) and Poil(T) is the oil production rate of the ith discrete unit of oil shale at the temperature of T (%).

Results and Discussion

Simulation Results of the Fluid Vector

Figure shows the fluid vector results of M-10-900. As can be seen in Figure , when the gas injection velocity is 10 m/s, the maximum fluid velocity in the entire flow field occurs at the intersection of the gas injection well and the fracture, reaching 117.7 m/s, which is much higher than the fluid velocity at the inlet of the gas injection well. This is due to the sharp decrease in space from the gas injection well of 150 mm to the fracture of 5 mm, consequently increasing the fluid velocity instantaneously. However, after entering the fracture, the fluid was subject to tremendous resistance, and the flow rate dropped obviously. After the gas flowed through the fracture and near the production well, the velocity started to increase and reached 79 m/s due to the low temperature and pressure in the production well. However, the flow velocity decreased again when the gas entered the production well because of the increased space. In addition, the flow velocity of a portion of gas that did not flow to the production well after entering the fracture became smaller and smaller. The entire injected gas flow process conforms to the principles of fluid mechanics. The equivalent volume diagram of fluid velocity in Figure S2 shows that the fluid velocity of the gas injection and production wells and the intersection of the wells and the fracture is 10 m/s or more. The equivalent volume area with a low fluid velocity of 4 m/s connected the entire area between the injection and production wells.

Figure 3

Fluid vector diagram of M-10-900.

Temperature Field Simulation Results

Figure shows the evolution of the temperature field of M-10-900, including both the transient values of 1 to 10 years and the steady values. The mass-weighted average temperature and the annual temperature increment of each of the transient and steady results are shown in Table . It can be found that a large part of oil shale is at a low temperature in the first 3 years and the T̅ of the first 3 years are all lower than 600 K. From the sixth year, the T̅ of oil shale has exceeded the ideal oil production temperature (643.5 K), indicating that the energy utilization efficiency of oil shale pyrolysis is relatively suitable at this time. However, the T̅ of oil shale has exceeded the excessively high temperature of 688.5 K from the ninth year, indicating that the energy utilization efficiency of oil shale pyrolysis began declining. As shown in Table , the annual temperature increment decreases year by year. The largest annual temperature increment is 225 K in the first year due to the high-temperature gradient between hot gas and oil shale in the early stage. As the oil shale was gradually heated to high temperatures, the temperature gradient decreased; meanwhile, the heat spread to the farther periphery of the injection well and production well. The annual temperature increment became less than 20 K since the sixth year. The 10th annual temperature increment is only 9.42 K. In the early stage of the reaction, the temperature difference between the hot gas and the oil shale layer is relatively large, which gradually decreased as the oil shale layer became warm, leading to the decrease in the heat absorbed by the layer, and thus the heating efficiency of oil shale, which needs to be fully considered in commercial operation, decreases as the time of gas injection increases.

Figure 4

Evolution of the temperature field of the M-10-900 case.

Table 2

Mass-Weighted Average Temperature and Annual Temperature Increment of the Oil Shale of M-10-900

time (year)	T̅ (K)	annual temperature increment (K)
1	525.01	225
2	563.59	38.59
3	593.23	29.64
4	617.54	24.31
5	637.71	20.17
6	655.16	17.45
7	669.92	14.76
8	682.83	12.91
9	693.99	11.16
10	703.40	9.42
Steady	740.67

The corresponding Vp, Vh, and Vn based on eqs –12 are presented in Table . It can be seen that the Vp increases year by year, and its growth rate first drops in the first 4 years and then rises. The Vp/Vpsteady exceeds 30% by the end of the 5th year, but it is only 51.47% in the 10th year, far from the ideal oil production temperature volume of the steady calculation result. The Vh also increased with increasing gas injection time, and its growth rate first drops in the first 7 years and then rises. The value of Vh/Vp remains between 73 and 80%, showing an upward trend in the first 4 years and then starting to decline. The smaller the value of Vh/Vp, the higher the energy efficiency of the oil shale pyrolysis process. The Vn decreases with the increase of gas injection time, and the Vn/Voil shale exceeds 60% in the first 5 years, indicating that less than 40% of oil shale is heated to a proper temperature for producing oil and gas. It takes more than 7 years to get 50% of oil shale to be heated to a proper temperature, and there is still more than 20% of oil shale unheated to a proper temperature by the end of the 10th year. This indicates that it takes a long time to have oil shale layer heated, and thus it is difficult to realize economic value with a short heating time for an in situ pyrolysis project.

Table 3

Annual Values of the Vp, Vh, Vn, and Vh/Vp

time (year)	Vp (m3)	<>Vp/Vpsteady<>a (%)	Vh (m3)	Vh/Vp (%)	Vn (m3)	Vn/Voil shaleb (%)
1	100.70	15.94	74.65	74.13	651.11	81.31
2	141.97	22.47	107.63	75.81	607.05	75.81
3	167.30	26.48	131.69	78.72	575.29	71.85
4	187.91	29.74	149.23	79.41	544.10	67.95
5	207.57	32.85	164.50	79.25	509.70	63.65
6	227.73	36.04	179.03	78.61	469.47	58.63
7	249.13	39.43	193.44	77.64	420.29	52.49
8	272.25	43.09	208.08	76.43	356.73	44.55
9	297.49	47.08	223.08	74.99	270.11	33.73
10	325.21	51.47	238.52	73.34	189.68	23.69
Steady	631.89		374.16	59.21	0.27	0.03

Vpsteady is the ideal oil production temperature volume for the steady calculation.

Voil shale is the total volume of oil shale.

Kerogen Content and Oil Yield

The transient kerogen content and oil production rate of M-10-900 were calculated based on eqs and 14, and the gas yield was converted into oil yield according to the heat value; the results are shown in Figure . It can be seen that, as the content of kerogen declines year by year, the oil and gas yield increases year by year, indicating that it is the kerogen that decomposes to oil and gas. The decrease rate of the kerogen content was rapid in the first few years, which reduced to less than 50% by the end of the 4th year, while there was still 21.82% of kerogen left in oil shale after 10 years. The total oil production rate exceeded 50% in the first 6 years, and the actual oil production rate for the 10 years was 67.86%, which is only 75.4% of the oil production rate of the steady value. This highlights that it still needs a long time to get the steady value as the decomposition rate of kerogen gradually decreases. Symington and Spiecker[46] and Han et al.[47] found that it would take almost 7–8 years to recover more than 90% of oil from oil shale by multiple heating wells.

Figure 5

Kerogen content and production of oil and gas.

Sensitivity Analysis

Two parameters of the numerical simulation model, that is, gas injection velocity (2, 5, 10, 15, and 20 m/s) and burning temperature (800, 850, 900, 950, and 1000 K), were selected to assess different sensitivities.

Effect of the Gas Injection Velocity

Figure shows the evolution of T̅, Vn, Vp, and Vh of cases with different gas injection velocities when the burning temperature is 900 K. The trend of T̅ over time is approximately the same in each case. However, the T̅ increases with increasing gas injection velocity. Basically, the higher the gas injection velocity, the higher the T̅, especially when the gas injection velocity increased from 2 to 5 m/s, whereas the T̅ of the 10th year increased by 112.2 K. However, the increased rate of the T̅ slowed down when the gas injection velocity increased over 10 m/s. The T̅ only increased by around 30 K when the gas injection velocity doubled from 10 to 20 m/s. Similarly, the effect of gas injection velocity on the Vp, Vh, and Vn increased as the gas injection time increased. Also, the higher the gas injection velocity, the higher the Vp and Vh and the lower the Vn, especially when the gas injection velocity increased from 2 to 5 m/s, whereas the increments of Vp, Vh and Vn in the 10th year were 121, 95.7, and −227.8 m3, respectively. However, the change of each indicator became smaller when the gas injection velocity was over 10 m/s. When the gas injection velocity increased from 15 to 20 m/s, the increment of Vp, Vh and Vn in the 10th year were only 20.6, 13.7, and −23.6 m3, respectively. As mentioned above, the DBHT injects room-temperature gas and uses a burner placed in the injection well to heat refluxing gas; when the gas velocity is low, increasing the gas injection velocity can increase the carried heat and accelerate the increase of the strata’s temperature. However, as the total heat is determined by the burning temperature of the burner, when the flow velocity increases to a certain extent that the heat carried by the gas reaches the maximum, the effect of the further increase in the flow velocity on heating oil shale layer is not significant. Therefore, it is meaningless to increase the injection rate further after it exceeds 10 m/s in the model discussed here.

Figure 6

Evolution of (a) mass-weighted average temperature T̅, (b) ideal oil production temperature volume Vp, (c) excessively high-temperature volume Vh, and (d) invalid temperature volume Vn with different gas injection velocities.

Figure shows the trends of kerogen content and oil and gas production over time in each case. It can be seen that the produced oil and gas increased as the content of kerogen decreased. As the gas injection time increased, the decomposition of kerogen and the production of oil and gas increased. Also, increasing the gas injection velocity can enhance the decomposition of kerogen. The kerogen content left in the oil shale in the 10th year decreased from 55 to 30% when the gas injection velocity increased from 2 to 5 m/s. Still, the increment in kerogen decomposition started reducing when the gas injection velocity continued increasing. At the same time, the oil and gas production increased as the gas injection velocity increased. The oil and gas production of the 10th year increased by 20.3% with the increase of the gas injection velocity from 2 to 5 m/s. Still, it only increased by 1.96% when the gas injection velocity increased from 15 to 20 m/s, consistent with the results reported by Pei et al.[48] that the cumulative oil equivalent was significantly improved by increasing the flow rate when the flow rate was low. Moreover, the effect of the gas injection velocity on oil and gas production is similar to that on the Vp (Figure b), that is, the higher the Vp, the more the oil and gas produced, indicating that it is not the mass-weighted average temperature but the volume of oil shale to reach the ideal oil production temperature that determines the size of the reaction zone where kerogen decomposes to generate oil and gas. This is because the kerogen decomposition occurs within a certain temperature range. The high mass-weighted average temperature may be caused by a local high temperature of the stratum, which has less impact on the ideal oil production temperature volume.

Figure 7

Evolution of (a) kerogen content and (b) oil and gas production with different gas injection velocities.

Effect of Burning Temperature

Five burning temperatures of 800, 850, 900, 950, and 1000 K are considered to assess its effect. In addition, the mass-weighted average temperature, ideal oil production temperature volume, excessively high-temperature volume, and invalid temperature volume of cases with different burning temperatures are shown in Figure .

Figure 8

Evolution of (a) mass-weighted average temperature T̅, (b) ideal oil production temperature volume Vp, (c) excessively high-temperature volume Vh, and (d) invalid temperature volume Vn with the different burning temperatures.

As shown in Figure , with the increase of the burning temperature from 800 to 1000 K, the T̅, Vn, and Vp increase and the Vh decreases. The increments/decrements are almost the same for every 50 K increase in burning temperature at any time but increases with the increase in the gas injection time. This indicates that the increase in the burning temperature is of benefit to the heat transfer. A higher heat transfer rate can shorten the effective production time and increase the process efficiency.[48]

Figure shows the evolution of kerogen content and oil and gas production of each case with the different burning temperatures. Also, it can be seen that the burning temperature has a positive effect on kerogen decomposition and oil and gas production. The kerogen content decreased, and the oil and gas production increased with increasing burning temperature. The increment and decrement of each volume were similar for each 50 K increase in burning temperature. This is because kerogen decomposition is a thermo-kinetically controlled mechanism; the higher the burning temperature, the more heat transfers to the stratum and the more rapidly the stratum approaches the pyrolysis temperature, resulting in faster recovery of hydrocarbons.[48,49]

Figure 9

Evolution of (a) kerogen content and (b) oil and gas production with the different burning temperatures.

Deformation and Stress–Strain Simulation Results

Considering that there is some equipment placed at the bottom of two wells, such as the burner in the injection well and the pumping unit in the production well, the deformation of strata or wells will have an adverse effect on the operation of the bottom equipment and recovery of the oil.

To analyze the stress–strain and deformation of M-10-900, we set the bottom surface of the model as a fixed surface and the front surface as a frictionless symmetrical surface. At the same time, the other four sides were constrained by their static pressure values calculated by the static pressure of rock mechanics. The changes in the physical and mechanical properties of oil shale with temperature and anisotropy as well as the thermal stress of the temperature field and the injection pressure were also considered. The intersection of two wells and the fracture as well as the interior of the fracture are chosen as the critical points to analyze its deformation and stress strain, and the results are shown in Table S1. It can be found in Table S1 that both the maximum and minimum principal stresses at Point A are negative, indicating that point A is under compressive stress, and so is Point C. The maximum principal stress at Point B is positive while its minimum principal stress is negative, and the absolute value of the minimum principal stress is much larger than the maximum principal stress, which indicates that Point B is mainly in the state of compressive stress. Therefore, the minimum principal stress (compressive stress) is critical to the oil shale in situ processing. The evolutions of minimum principal stress, minimum principal strain, and total deformation over time for Point A, Point B, and Point C are shown in Figure .

Figure 10

Evolution of (a) minimum principal stress, (b) minimum principal strain, and (c) total deformation over gas injection time.

With the increase of gas injection time, the minimum principal stress and strain at Point A decreased in the first 7 years and then slightly increased, and its minimum principal stress and strain of the steady results are similar to those after the first year. A similar trend can also be found at Point B, whose minimum principal stress and strain decreased in the first 9 years and then started to increase, and the minimum principal stress and strain of the steady results of Point B are almost the same as those of the third year. According to the steady results, although it kept decreasing in the first 10 years, the minimum principal stress and strain of Point C should increase later. The total deformations of Point A, Point B, and Point C showed a decrease first and then increased. The lowest total deformations of Point A, Point B, and Point C were 0.03, 0.009, and 0.012 m in the fourth, fifth, and sixth years, respectively. In the first 4 years, the total deformation at Point A is the smallest one, followed by that of Point B, and the total deformation at Point C is the largest. After the sixth year, the order of the total deformation of these three points is Point A > Point B > Point C. Besides, the total deformations of the 10th year of Point A, Point B, and Point C were 0.121, 0.090, and 0.079 m, respectively, lower than those of the steady values of 0.149, 0.115, and 0.107 m, respectively.

The evolution of the total deformation of M-10-900 is shown in Figure . The graph in the upper right corner of the annual deformation graph is the corresponding isosurface graph. It can be seen that the total deformation in the first year was concentrated around the injection well, and the maximum deformation was 0.13 m. The maximum deformation of the second year increased up to 0.15 m. The deformation zone further increased along the fracture from the third to the sixth year with the increasing maximum deformation by 0.04 m. The maximum deformation increased faster since the seventh year, and it increased to 0.24 m by the end of the 10th year, which was slightly smaller than that of the steady value. However, there is still a significant difference in the cloud map of the total deformation between the 10th year and the steady results, indicating that the deformation of the stratum could be a critical problem during the oil shale in situ conversion process, especially in the long-term operational lifetime. It is necessary to consider the deformation at the bottom of the gas injection well and production well, which may have an adverse effect on the operation of the equipment there.

Figure 11

Cloud maps of the total deformation of M-10-900 over time.

Conclusions

In this paper, the numerical simulation with a 3D model of oil shale in situ processing by applying the idea of DBHT was conducted to examine its temperature field, thermal processes, and stress–strain and deformation. The main conclusions and recommendations are given below:

The maximum fluid velocity in the entire flow field occurred at the intersection of the gas injection well and the fracture.

The mass-weighted average temperature of the oil shale layer increased with increasing gas injection time, but the annual temperature increment decreased year by year. The volume of oil shale to reach the ideal oil production temperature determines the size of the reaction zone where kerogen decomposes to generate oil and gas.

The oil and gas yields rose as the kerogen content decreased. Increasing the gas injection velocity can accelerate the kerogen decomposition due to the enhanced heat transfer efficiency, but this acceleration becomes insignificant when the gas injection velocity exceeds 10 m/s. Also, increasing the burning temperature can significantly enhance the recovery rate of hydrocarbon, which is an effective way to accelerate kerogen decomposition and oil production.

For the long-term operational oil shale in situ conversion field pilot, the deformation of the oil shale formation could be a critical issue and should be considered at the design stage.