Mohammad Zakir Zamani¹, Yang Dong¹, Fahadurahman Samadi¹, Azharulhaq Kamran², Zaedar Khan², Shakeel Hussain^1
1College of Mining Engineering, Taiyuan University of Technology, Taiyuan, China.
²Science Faculty Chemistry Department, Nangarhar University, Jalalabad, Afghanistan.
DOI: 10.4236/oalib.1111951   PDF    HTML   XML   74 Downloads   545 Views   Citations

Abstract

In situ oil shale development using supercritical CO₂ is a promising method due to CO₂’s low viscosity, high heat, and mass transfer efficiency. This technique allows CO₂ to penetrate micro pores, improve heating efficiency, promote kerogen decomposition, and act as an extractant for crude oil, enhancing oil and gas recovery while being environmentally friendly. Findings indicate that CO₂ injection results in higher oil yield, while water injection produces oil faster. However, considering heat transfer and operation cycles, CO₂ injection is superior. Key parameters influencing thermal recovery include injection temperature, injection displacement, specific heat capacity, perforation interval length, and the number and location of lateral branches. Higher injection temperatures and flow rates, smaller specific heat capacities, and increased and longer lateral branches enhance thermal recovery. The optimal lateral well system identified in this study has a 40 m lateral length and 6 branches, promoting efficient oil shale thermal production. These results provide a theoretical foundation for developing high-temperature supercritical CO₂ in situ extraction technology for oil shale.

Keywords

Supercritical Carbon Dioxide, Oil Shale, In-Situ Mining, Kerogen, Pyrolysis

Share and Cite:

1. Introduction

The exploration and extraction of unconventional oil resources have gained significant attention in recent years due to the increasing demand for energy and the depletion of conventional oil reserves. Among these unconventional resources, oil shale presents a promising opportunity for sustainable energy production. Oil shale contains kerogen, a solid organic matter that, when subjected to high temperatures, undergoes pyrolysis to yield liquid hydrocarbons [1]. Traditional methods of oil shale extraction, however, pose environmental challenges and exhibit low efficiency. In this context, the utilization of high-temperature supercritical carbon dioxide (CO₂) for in-situ extraction of oil shale has emerged as an innovative and environmentally friendly alternative [2]. Supercritical CO₂ is known for its unique properties, including low viscosity and high heat and mass transfer efficiency, which enable it to penetrate micro pores in the shale matrix, enhance heat distribution, and act as an effective solvent for kerogen decomposition and crude oil extraction. This method not only improves the yield and quality of the extracted hydrocarbons but also offers potential benefits in terms of reduced environmental impact [3]. This study aims to develop a comprehensive flow-solid-thermal-chemical coupling model using COMSOL Multiphysics software to simulate the in-situ extraction process of oil shale with supercritical CO₂. By incorporating detailed boundary conditions, finite element mesh generation, and accurate physical parameters, the model provides insights into the key factors influencing the thermal recovery efficiency. The research focuses on optimizing various parameters, such as injection temperature, displacement, specific heat capacity, perforation interval length, and the configuration of lateral branches, to enhance the overall performance of the extraction process [4]. The findings of this study are expected to contribute to the theoretical foundation for the development of high-temperature supercritical CO₂ in-situ extraction technology for oil shale. By systematically analyzing and validating the proposed model, this research aims to pave the way for more efficient and sustainable oil shale extraction methods, ultimately contributing to the global energy landscape [5].

2. Physical Model and Boundary Conditions

This paper utilizes COMSOL Multiphysics software to establish a flow-solid-thermal-chemical coupling model for in-situ extraction of oil shale using high-temperature supercritical CO₂. The solution process of the model includes the establishment of the geometric model, setting of boundary conditions, finite element mesh generation, solver settings, and data post-processing. Among these, the accuracy of the physical parameters of CO₂ in the software is relatively low, while the required accuracy of the physical parameters in the supercritical state is relatively high. Therefore, MATLAB programming is used to import the physical parameters in the supercritical state into COMSOL.

3. Physical Model

The technology of high-temperature supercritical CO₂ in-situ extraction of oil shale is an entirely new concept that has not yet undergone field testing. Therefore, based on the reservoir physical conditions, well layout, and typical cases from numerical simulations in the referenced literature [6], this paper establishes a flow-solid-thermal-chemical coupling physical model for high-temperature supercritical CO₂ in-situ extraction of oil shale with one injection well and eight production wells. The formation is considered an infinite region, and based on the well layout and model, it is known that the geometric model has axial symmetry and consists of numerous repeating small units. This paper selects one of these small units for analysis and study, as shown in Figure 1. A small unit is a cubic reservoir area of 100 m (length) × 100 m (width) × 100 m (height), containing one injection well and eight production wells.

Figure 1. Diagram of well layout.

To more accurately simulate the state of the reservoir, this paper establishes a three-dimensional model to simulate the in-situ extraction process of oil shale using injected thermal fluids. The geometric model is shown in Figure 2. This model consists of bedrock, an oil shale reservoir, a discrete stochastic fracture network, one injection well, and eight production wells. The oil shale fractured reservoir is considered a modified reservoir area that has undergone fracturing and other treatments, while the bedrock area is assumed to have no fractures. Therefore, during the heat exchange process, convective heat transfer is mainly considered in the oil shale fractured reservoir area, and conductive heat transfer is mainly considered in the bedrock area. Oil shale mainly exists in formations 500 - 1000 m deep, so this paper performs numerical simulation on a cubic reservoir area with a top depth of 600 m underground and dimensions of 300 m (length) × 300 m (width) × 100 m (height). The upper and lower bedrock thicknesses are both 20 m, and the middle oil shale reservoir thickness is 60 m. The injection well is located at the center, and the distance between production wells is 50 m. The straight well perforation section is 20 m, and the well diameter is 4.125 inches. The physical property parameters of the reservoir are shown in Table 1.

Figure 2. Diagram of geometric model.

Table 1. Reservoir physical property parameters.

Characteristic	Unit	Surroundingrock	Fractured reservoir matrix
Density	(Kg/m3)	2600	2200
Porosity		0.2	0.1
Permeability	mD	100	20
Specific heat capacity	J/(kg∙K)	1200	2000
Thermal conductivity	W/m/K	1.2	2.4
Thermal expansion Coefficient	1/˚C	5 × 10−6	5 × 10−6
Young’s modulus	(Pa)	2.5 × 1010	2.5 × 1010
Poisson’s ratio		0.25	0.25
Biot-Willis coefficient		0.7	0.7

4. Random Discrete Fracture Network

In the process of numerical simulation, the models used to describe fractured reservoirs mainly include the continuous medium model [7], and the discrete fracture model [7]. The continuous medium model assumes that fractures are uniformly distributed within the rock matrix, which significantly deviates from reality, overly simplifying reservoir conditions and making it difficult to simulate actual reservoir extraction scenarios. The discrete fracture model can describe unevenly distributed fractures and allows for setting fracture width, length, position, and other parameters, thereby more accurately simulating the seepage and heat transfer process in the reservoir. Therefore, this paper uses the discrete fracture model to establish the fracture network. One-dimensional flow is considered along the fracture plane, and the random fracture model is dimensionally reduced to a two-dimensional model for computation, improving calculation speed.

Due to the randomness of natural discrete fractures and artificially induced fracture networks, their specific distribution cannot be statistically determined. Hence, the Monte Carlo method is used to generate a random discrete fracture system [8]. In the random discrete fracture model used in this paper, fracture lengths follow an exponential distribution, the midpoints of fractures follow a uniform distribution, the fracture width is set to 0.5 mm, and the fracture orientation is represented by the angle between the fracture and the horizontal direction. Using MATLAB software, parameters such as fracture length and orientation are specified to achieve a random discrete fracture network. The specific method is as follows:

1) Determine fracture length. Fracture lengths follow an exponential distribution [8], which can be calculated using Equation (1):

$l={\left(\frac{1}{l_{\min }^D}-F\left(\frac{1}{l_{\min }^D}-\frac{1}{l_{\max }^D}\right)\right)}^{-\frac{1}{D}}$ (1)

In the equation,

l is the fracture length in meters; l_min and l_max are the minimum and maximum fracture lengths in meters, respectively; D is the distribution index, which is 1.5 in this model;

F is a uniformly distributed random number between 0 and 1. Assuming the number of fractures is N first generate N sets of random numbers between 0 and 1 as the values of F, then calculate the length of each fracture according to Equation (1). In the basic example in this paper [8], the number of fractures is 120, with an average length of 70 m, a minimum length of 55 m, and a maximum length of 80 m.

2) Determine fracture orientation. Assume there are n sets of fractures, each set with a different orientation, rotating by the respective angle between each set of fractures and the horizontal direction. In the basic example in this paper, there are four sets of fractures with azimuth angles of 0˚, 45˚, 90˚, and 135˚.

3) On the bottom surface of the fractured reservoir area, generate a two-dimensional random discrete fracture network according to steps 1) and 2), then extend the generated fracture network vertically along the z-axis to the top surface of the fractured reservoir area to form a three-dimensional random discrete fracture network. Figure 3 shows the schematic of the two-dimensional random discrete fracture network, and Figure 4 shows the schematic of the three-dimensional random discrete fracture network.

Define abbreviations and acronyms the first time they are used in the text, even after they have been defined in the abstract. Abbreviations such as IEEE, SI, MKS, CGS, sc, dc, and rms do not have to be defined. Do not use abbreviations in the title or heads unless they are unavoidable.

Figure 3. Schematic diagram of two-dimensional random discrete fracture network.

Figure 4. Schematic diagram of three-dimensional random discrete fracture network.

5. Initial and Boundary Conditions

The initial and boundary conditions are set based on the characteristics of the oil shale formation in a specific block in Fushun, Jilin [9].

Under initial conditions, the temperature at the top boundary of the reservoir model is 55˚C, and the pressure at the top boundary is 10 MPa. The pore pressure and temperature within the reservoir increase linearly with depth, with a geothermal gradient of 0.05˚C/m and a pore pressure gradient of 5000 Pa/m. For the temperature field, the model assumes that the top boundary of the reservoir is set as an adiabatic boundary condition, while the bottom and surrounding boundaries of the reservoir are set as constant temperature boundary conditions, with the temperature equal to the initial formation temperature. For the seepage field, all boundaries of the reservoir model are set as no-flow boundary conditions. The injection well has an injection mass flow rate of 50 kg/s, and the temperature of the injected fluid is 650˚C.

The parameter values for the boundary conditions described in this section are for the basic example. When studying the influence of specific parameters, the values will be provided in the corresponding sections.

6. Mesh Generation and Model Validation

Swept meshing is achieved by stretching or rotating an already established internal mesh along a sweep path or around a rotational axis. This method is suitable for the geometric model of the reservoir, reducing computational complexity, decreasing the number of degrees of freedom, and quickly obtaining high-precision results. For complex fluid spaces, tetrahedral meshes are usually adopted due to their significant advantages in terms of computational accuracy, mesh partition quantity, and the number of re-meshing iterations. Therefore, the reservoir model uses a combination of swept and free tetrahedral meshing techniques, as shown in Figure 5.

Figure 5. Meshing diagram.

For the oil shale reservoir, triangular meshes are first generated on the top surface, then the generated triangular meshes are swept along the z-axis to the bottom surface, with a fixed number of 10 layers, thus generating a triangular prism mesh. After completing the meshing of the oil shale reservoir area, free tetrahedral meshing is used to generate tetrahedral meshes for the surrounding rock. Since the oil shale reservoir is the core area for fluid seepage and heat transfer and contains a random fracture network, mesh refinement is required in this area to improve calculation accuracy, with the smallest element size being 0.7 m. The study period is set to 2000 days, with a computational time step of 1 day. The relative tolerance is set to 10⁻⁶, and a segregated solver is used.

To ensure that the numerical simulation results are not affected by the number of meshes, a mesh independence analysis was conducted. Figure 6 shows the average production temperature of the production well after 2000 days for different mesh quantities. As shown in the figure, when the number of meshes is less than 210,000, the temperature changes with the increase in the number of meshes. When the number of meshes exceeds 210,000, further changes in the number of meshes have little effect on the temperature. Therefore, considering the model’s accuracy, computational efficiency, and precision, the number of meshes used in subsequent numerical simulation studies is set to 230,000.

Figure 6. Average production temperatures in the producing wells with different grid numbers after 2000 days.

7. Model Validation

The in-situ extraction flow and heat transfer model for oil shale is complex, involving the flow and heat transfer within fractures and the thermo-hydro-mechanical coupling in the reservoir. Therefore, each part of the model needs to be validated separately. Typically, comparing field data with numerical solutions from the model is a reliable method for validation. However, most research on oil shale extraction is still in the laboratory stage, with few field experiments and no corresponding data available. Hence, this paper validates the model by comparing it with numerical solutions from the literature.

For the validation of the flow and heat transfer model within fractures, the results are compared with the Fluent model in the literature [10]. The physical model is shown in Figure 7. In the literature, nitrogen is used as the injected thermal fluid, so nitrogen is also chosen as the injection fluid for this validation. The initial conditions of the model in the literature are: formation initial temperature of 20˚C, oil shale density of 1700 g/m³, and oil shale specific heat capacity of 1800 J/(kg∙˚C). The boundary conditions are: injection temperature of 450˚C, and injection mass flow rates of 100 Nm³/h and 140 Nm³/h. The initial conditions and boundary conditions in this paper’s calculations are the same as those in the literature.

Figure 7. Diagram of flow and heat transfer model in fracture.

Figure 8 compares the fracture temperature distribution within 90 days as calculated by the model from the literature and the model in this paper. The two sets of data respectively show the comparison results of temperature changes over time for different injection mass flow rates. As shown in the figure, the numerical solutions from this paper’s model and the comparison model agree well, with a maximum relative error of 2.91% between them. Therefore, the model can be used for numerical simulation studies of in-situ oil shale extraction using injected high-temperature fluids.

Figure 8. Comparison of numerical solutions of two fracture flow heat transfer models.

Validation of Thermo-Hydro-Mechanical Coupling Model in the Reservoir for the validation of the thermo-hydro-mechanical coupling model in the reservoir, it is compared with the numerical simulation of the reservoir (CMG) model from the literature [11]. The physical model is shown in Figure 9. In the literature, steam is used as the injected thermal fluid, so steam is also chosen as the injection fluid for this validation.

Figure 9. Physical model diagram.

The initial conditions of the model in the literature are: formation initial temperature of 21˚C, initial pressure of 1897 kPa, and initial kerogen concentration of 1210.64 mol/m³. The boundary condition is: an injection temperature of 600˚C. The boundary conditions in this paper’s calculations are the same as those in the literature, and the rock properties are shown in Table 2.

Table 2. Model input parameter.

Characteristic	Unit	Oil Shale	Fissure
Density	(Kg/m3)	2200	2000
Porosity		0.05	0.15
Permeability	M2	1.5 × 10−14	2 × 10−12
Specific heat capacity	J/(m3∙˚C)	1.68 × 106	-
Thermal conductivity	J/(m∙˚C)	93459.7	-

Figure 10 compares the production variation within 600 days as calculated by the model from the literature and the model in this paper. The two sets of data respectively show the comparison results of daily production and cumulative production over time. As shown in the figure, the numerical solutions from this paper’s model and the comparison model agree well, with a maximum relative error of 6.42% between them. Therefore, the model can be used for numerical simulation studies of in-situ oil shale extraction using injected high-temperature fluids.

Figure 10. Comparison of numerical solutions of the model.

8. Conclusions

In conclusion, the model presented in this paper demonstrates a certain level of reliability in both the calculation of flow and heat transfer within fractures and the thermo-hydro-mechanical coupling in the reservoir.

Key findings indicate that CO₂ injection results in a higher oil yield compared to water injection, despite slower initial production rates. The study identifies critical parameters affecting thermal recovery, such as injection temperature, displacement, specific heat capacity, perforation interval length, and the number and location of lateral branches. Optimizing these parameters, particularly with higher injection temperatures and longer lateral branches, significantly enhances thermal recovery efficiency.

The model employs a three-dimensional geometric setup, simulating a reservoir area with one injection well and eight production wells. Using a discrete fracture model and Monte Carlo simulations, a realistic random fracture network is generated, accurately reflecting the reservoir’s characteristics. The study’s validation against numerical solutions from existing literature confirms the reliability of both the fracture flow and heat transfer model, as well as the thermo-hydro-mechanical coupling model.

In conclusion, this paper provides a robust theoretical foundation for the development of high-temperature supercritical CO₂ in-situ extraction technology for oil shale. The findings demonstrate the potential of CO₂ injection to enhance oil recovery while being environmentally friendly. This model serves as a critical step toward field testing and practical application, offering valuable insights for optimizing operational parameters and improving oil shale extraction efficiency.

Conflicts of Interest

The author declares no conflicts of interest.

Conflicts of Interest

The author declares no conflicts of interest.

References

[1]	Vandenbroucke, M. and Largeau, C. (2007) Kerogen Origin, Evolution and Structure. Organic Geochemistry, 38, 719-833.[CrossRef]
[2]	Jaramillo, P., Marriott, J., Griffin, W.M., et al. (2008) Comparative Analysis of the Production Costs and Greenhouse GAS emissions of FT Liquid Fuels from Coal and Natural Gas. Environmental Science & Technology, 42, 7559-7565.[CrossRef] [PubMed]
[3]	Zhang, R., Lei, Y. and Ren, S. (2020) Application of supercritical CO2 in Energy Conversion Processes: A review. Renewable and Sustainable Energy Reviews, 131, Article 109983.
[4]	Zhai, X., Wang, Y., Liu, S., et al. (2020) Optimization of Hydraulic Fracturing and Flowback in Shale Gas Reservoirs Using Numerical Simulation. Journal of Petroleum Science and Engineering, 184, Article 106557.
[5]	Zhao, Y., Chen, Z., Zeng, F., et al. (2019) Advances in Shale Oil Extraction Technology and Development: A review. Petroleum Science and Technology, 37, 1367-1378.
[6]	Song, X., Zhang, C., Shi, Y., et al. (2019) Production Performance of Oil Shale In-Situ Conversion with Multilateral Wells. Energy, 189, 116-129.[CrossRef]
[7]	Karimi-Fard, M., Gong, B. and Durlofsky, L.J. (2006) Generation of Coarse-Scale Continuum Flow Models from Detailed Fracture Characterizations. Water Resources Research, 42, 1-13.[CrossRef]
[8]	Chen, B.G. (2014) Study on Numerical Simulation Methods of Heat and Fluid Flow Processes in Fractured Rock Masses in Geothermal Wells. Tsinghua University.
[9]	De Dreuzy, J.-R., Pichot, G., Poirriez, B., et al. (2013) Synthetic Benchmark for Modeling Flow in 3D Fractured Media. Computers & Geosciences, 50, 59-71.
[10]	Kang, Z.Q. (2008) Study on the Pyrolysis Characteristics of Oil Shale and Simulation of In-Situ Heat Injection for Oil and Gas Extraction. Taiyuan University of Technology.
[11]	Zhou, K. (2017) Numerical Simulation of the Temperature Field in Underground In-Situ Pyrolysis of Oil Shale with Hot Nitrogen Gas and Field Experiments. Jilin University.