Cyclic Subcritical Water Injection into Bazhenov Oil Shale: Geochemical and Petrophysical Properties Evolution Due to Hydrothermal Exposure

Cyclic Subcritical Water Injection into Shale: Geochemical Petrophysical Properties Evolution Abstract: In situ shale or kerogen oil production is a promising approach to developing vast oil shale resources and increasing world energy demand. In this study, cyclic subcritical water injection in oil shale was investigated in laboratory conditions as a method for in situ oil shale retorting. Fifteen non-extracted oil shale samples from Bazhenov Formation in Russia (98 ◦ C and 23.5 MPa reservoir conditions) were hydrothermally treated at 350 ◦ C and in a 25 MPa semi-open system during 50 h in the cyclic regime. The inﬂuence of the artiﬁcial maturation on geochemical parameters, elastic and microstructural properties was studied. Rock-Eval pyrolysis of non-extracted and extracted oil shale samples before and after hydrothermal exposure and SARA analysis were employed to analyze bitumen and kerogen transformation to mobile hydrocarbons and immobile char. X-ray computed microtomography (XMT) was performed to characterize the microstructural properties of pore space. The results demonstrated signiﬁcant porosity, speciﬁc pore surface area increase, and the appearance of microfractures in organic-rich layers. Acoustic measurements were carried out to estimate the alteration of elastic properties due to hydrothermal treatment. Both Young’s modulus and Poisson’s ratio decreased due to kerogen transformation to heavy oil and bitumen, which remain trapped before further oil and gas generation, and expulsion occurs. Ultimately, a developed kinetic model was applied to match kerogen and bitumen transformation with liquid and gas hydrocarbons production. The nonlinear least-squares optimization problem was solved during the integration of the system of differential equations to match produced hydrocarbons with pyrolysis derived kerogen and bitumen decomposition.


Introduction
Due to the worldwide increase in energy consumption caused by rapid industrial development, the request for alternative hydrocarbon resources has emerged. However, despite the energy problem, oil is a strategic source of raw material for the petrochemical industry today. Driven by the growing demand, the petroleum industry innovated its way into the unconventional deposits of bitumen, heavy oil, and shale oil, which have changed the oil and gas industry's structure during the last decade [1]. The application of thermal enhanced oil recovery (EOR) methods today is one of the most efficient ways to significantly increase the recovery factor, especially in oil shale, heavy oil, and bitumen reservoirs [2,3].
A major concern of heating agent injection into oil shale reservoir is related to heat losses to surrounding formation. However, recently, it was shown numerically that in high pressure and high temperature conditions, initial steam quality, fluid flow rate, and insulation's thermal conductivity are vital parameters influencing the further water phase state along the tubing [29]. Hydrostatic and frictional pressure gradients are significant contributors to changes in thermodynamic properties and temperature profile along the well. However, cost optimization includes not only accounting for heat losses but also finding optimal technological parameters (temperature, injection and soaking time duration, etc.), which are dependent on kinetics of organic matter decomposition in the presence of water.
Hydrous pyrolysis or pyrolysis in the presence of water was first used by Lewan et al. to describe oil-generation kinetics of source rocks [30]. Lewan showed that in the case of pyrolysis of type II kerogen, water promotes thermal cracking reactions and inhibits carbon-carbon bond cross-linking compared to anhydrous conditions. It was observed that pyrolysate collected from the water surface at the end of the experiment has a composition of saturates, aromatics, resins, and asphaltenes (SARA) similar to natural petroleum [31]. Later, a comparative analysis of type III kerogen pyrolysis in open-system, closed-system anhydrous, and hydrous conditions was carried out [32]. The reaction products of closedsystem hydrous pyrolysis contained less hydrocarbon gas than products of the anhydrous open-system pyrolysis. The yield of liquid hydrocarbons was maximum after open-system anhydrous and closed-system hydrous pyrolysis, whereas it was three times lower under closed-system anhydrous conditions. The observed effect was attributed to the extra-high yield of aromatic and polar compounds compared to closed-system anhydrous pyrolysis [32]. Water as a heat carrier also has an important disadvantage: at very high temperatures and high pressure, there will be a lot of mutual solubility of water and oil even in heavy oils and bitumen [33]. Industrial application of the approach requires taking into account the problem of separating the fluids.
Many studies on numerical modeling are devoted to oil shale in situ retorting procedure development, including the following: acceleration of simulations using fast upscaled dynamic models [34], parameter space reduction and sensitivity analysis [35], reaction parameters estimation using inverse modeling [36], simulation with fully coupled thermal-hydraulic-mechanical properties [37], application of various operator splitting techniques [38], different kinetic models and heat delivery methods application [39][40][41].
However, to specify the optimal in situ oil shale retorting conditions, a few critical effects should be accounted for, such as developing a reliable kinetic model associated with organic matter conversion, natural porosity, and permeability, and the evolution of geomechanical properties during in situ retorting. We report the results of cyclic subcritical water injection into Bazhenov oil shale samples at 350 • C and 25 MPa and provide a description of the alteration of geochemical, elastic, and microstructural properties after treatment. Additionally, a common kinetic model was tuned to match kerogen and bitumen transformation with the produced hydrocarbons.

Subcritical Water Injection
In the experiment, 460.8 g of non-extracted oil shale samples were introduced into a one-liter Hastelloy C-276 autoclave (Parr Instrument Company, Moline, IL, USA) to simu- late subcritical water injection at reservoir conditions. After the reactor was sealed, 900 mL of distilled and deionized water was injected to implement complete water contact with the oil shale chips and maintain pressure at 7 MPa. After checking for leaks, temperature was raised to 350 • C in 60 min with an average slope of 5 • C/min. During heating, the autoclave pressure was maintained at 25 MPa using a back-pressure regulator to simulate reservoir pressure (Figure 1). Exposure temperature was chosen following the results of previous studies, which showed significant liquid hydrocarbons yield in laboratory conditions, only starting from 350 • C [25].

Subcritical Water Injection
In the experiment, 460.8 g of non-extracted oil shale samples were introduced into a one-liter Hastelloy C-276 autoclave (Parr Instrument Company, Moline, Illinois, USA) to simulate subcritical water injection at reservoir conditions. After the reactor was sealed, 900 mL of distilled and deionized water was injected to implement complete water contact with the oil shale chips and maintain pressure at 7 MPa. After checking for leaks, temperature was raised to 350 °C in 60 min with an average slope of 5 °C/min. During heating, the autoclave pressure was maintained at 25 MPa using a back-pressure regulator to simulate reservoir pressure (Figure 1). Exposure temperature was chosen following the results of previous studies, which showed significant liquid hydrocarbons yield in laboratory conditions, only starting from 350 °C [25]. The back-pressure regulator allowed the experiment to be conducted in a self-purging regime to simulate true in situ retorting. The generated products remained in the reactor's heated part until the pressure buildup was sufficient to produce them. The test consisted of 9 five-hour cycles with a pressure drop by 7 MPa at the end of each stage to produce generated gas and liquid hydrocarbons. During each cycle, the exposure period at reservoir pressure and temperature was followed by pressure decrease and water flushing at a rate of 3 mL/min to produce generated hydrocarbons ( Figure 1). After the experiment, the autoclave was cooled down to ambient temperature with simultaneous gradual pressure decrease over 4 h to prevent water phase change. The experimental system was disassembled; lines were rinsed with toluene. Produced gas and liquids, spent shale samples were collected for further analysis.

Rock-Eval Pyrolysis
Programmed Rock-Eval pyrolysis was performed by HAWK (Wildcat Technologies, Humble, Texas, USA) on both original shale samples pre-and post-extraction, and after exposure to subcritical water pre-and post-extraction. Table 1 summarizes the properties derived from the Rock-Eval analysis of fifteen non-extracted oil shale samples before and after the experiment. Parameters S0 and S1 indicate gas and light liquid hydrocarbons adsorbed on the organic or mineral surface, while parameter S2 demonstrates the amount of hydrocarbons that could be generated through pyrolysis of heavy oil, bitumen, and kerogen. The amount of CO2 released during pyrolysis is denoted as S3. The back-pressure regulator allowed the experiment to be conducted in a self-purging regime to simulate true in situ retorting. The generated products remained in the reactor's heated part until the pressure buildup was sufficient to produce them. The test consisted of 9 five-hour cycles with a pressure drop by 7 MPa at the end of each stage to produce generated gas and liquid hydrocarbons. During each cycle, the exposure period at reservoir pressure and temperature was followed by pressure decrease and water flushing at a rate of 3 mL/min to produce generated hydrocarbons ( Figure 1). After the experiment, the autoclave was cooled down to ambient temperature with simultaneous gradual pressure decrease over 4 h to prevent water phase change. The experimental system was disassembled; lines were rinsed with toluene. Produced gas and liquids, spent shale samples were collected for further analysis.

Rock-Eval Pyrolysis
Programmed Rock-Eval pyrolysis was performed by HAWK (Wildcat Technologies, Humble, TX, USA) on both original shale samples pre-and post-extraction, and after exposure to subcritical water pre-and post-extraction. Table 1 summarizes the properties derived from the Rock-Eval analysis of fifteen non-extracted oil shale samples before and after the experiment. Parameters S 0 and S 1 indicate gas and light liquid hydrocarbons adsorbed on the organic or mineral surface, while parameter S 2 demonstrates the amount of hydrocarbons that could be generated through pyrolysis of heavy oil, bitumen, and kerogen. The amount of CO 2 released during pyrolysis is denoted as S 3 . Table 1. Rock-Eval pyrolysis parameters of the non-extracted oil shale samples before / after the experiment.

Three-Dimensional X-ray Computed Microtomography
X-ray computed microtomography is a non-invasive and non-destructive laboratory method to characterize the properties and structure of the pore space. Phoenix v|tome|x L240 (General Electric, Boston, Massachusetts, USA) equipped with unipolar 240 kV micro-focus and 180 kV high-power nano-focus X-ray source allows the structure of oil shale samples to be scanned and reconstructed with up to 1 µm resolution. The oil shale samples and thermally treated (green) non-extracted oil shales; (b) S 2 and S 1 Rock-Eval parameters before the experiment and extraction (red), before the experiment and after extraction (blue), after the experiment and before extraction (brown), after the experiment and after extraction (green).

Three-Dimensional X-ray Computed Microtomography
X-ray computed microtomography is a non-invasive and non-destructive laboratory method to characterize the properties and structure of the pore space. Phoenix v|tome|x L240 (General Electric, Boston, MA, USA) equipped with unipolar 240 kV micro-focus and 180 kV high-power nano-focus X-ray source allows the structure of oil shale samples to be scanned and reconstructed with up to 1 µm resolution. The oil shale samples were mounted on a rotary stage, irradiated by micro-focus polychromatic X-ray source, and imaged on a detector.

Elastic Properties
The samples' elastic properties were studied using pulse transmission acoustic experiments on PIK-UZ (Geologika, Novosibirsk, Russian Federation). Acoustic wave velocities were measured at ambient conditions and then recalculated to Young's module and Poisson's ratio. Measurements were performed on non-extracted oil shale samples before and after the experiment with axial stress applied parallel to bedding to study the variation after artificial thermal maturation by subcritical water injection.

Product Analysis
The volume of produced gases was monitored by the gas meter. Helium with a 100 mL/min flow rate was used in the experiment as a carrier gas to deliver produced gases to the gas chromatograph (Agilent 7890B, Santa Clara, CA, USA). Chromatograms were recorded in 20 min intervals. The gas vent between the condenser and gas meter (Figure 1) was opened only during the pressure drop at the end of each cycle.
The produced hydrocarbons and water were separated using chloroform and a separating funnel. The SARA composition of hydrocarbons was determined by Iatroscan Thin-Layer Chromatography with Flame Ionization Detector (TLC-FID) using solvents: n-heptane, toluene:n-heptane 80:20 mixture, dichloromethane:methanol 95:5 mixture according to IP 469. The boiling point distribution of produced hydrocarbons was obtained by Simulated Distillation (SimDis) gas chromatography according to ASTM D2887.
After the experiment, the spent shale samples were dried in an oven at 105 • C for 24 h. Bitumen was extracted with chloroform from several samples to study the contribution of kerogen and bitumen to the generation of hydrocarbons separately. Rock-Eval pyrolysis was used to characterize the organic matter transformation after the experiment.

Laboratory Investigation
The open-system Rock-Eval pyrolysis results showed how organic matter distribution changed after the test (Table 1, Figure 2). Parameter S 0 decreased on average by 89.4%, S 1 and S 2 also lost 86.3% and 73.2%, respectively. T max mean value increased from 440 • C to 456 • C after subcritical water injection. Hydrogen index mean value reduced from 294 to 104, and the production index mean value dropped 2.2 times from 0.163 to 0.075.
Gas was generated along with oil from bitumen, kerogen, and various minerals within a sample. The GC results showed the dominant production of non-hydrocarbon gases (carbon dioxide and hydrogen sulfide). CO 2 was possibly generated due to organic matter decarboxylation reactions [42,43]. However, some studies show the possibility of CO 2 formation from carbonate mineral decomposition [43,44]. Hydrogen sulfide was conceivably produced due to organic compounds desulfurization reactions and pyrite to pyrrhotite transformation [45][46][47][48]. As for hydrocarbon gases, methane, ethane, propane, and butane each showed about 120 mg cumulative mass yield ( Figure 3). It should be noted that according to Figure 3, maximum molar fraction peak shift is observed depending on the compound: methane and ethane have a maximum generation peak at 20 h, propane and butane at 15 h, pentane and hexane at 5-10 h. The effect could be explained due to more energy (exposure time) needed to crack large molecules into smaller fragments. The non-monotonic gas production behavior could also be partially caused by the overlay effect of gas desorption, microfractures appearance with further gas release from non-isolated pores, and thermal cracking of generated hydrocarbons.  Liquid hydrocarbons yield also showed non-monotonic behavior with a maximum at the fifth cycle (25 h). This could be attributed to the overlay effect of hydrocarbons desorption, appearance of microfractures due to thermal expansion, continuous hydrocarbon release from non-isolated pores, and thermal cracking of kerogen, bitumen, and generated hydrocarbons. However, a substantial amount of hydrocarbons was obtained after rinsing the autoclave with toluene ( Figure 4). Results of the simulated distillation of the produced liquid hydrocarbons (Figure 4) showed the insignificant contribution of >C35 fraction of hydrocarbons, which did, however, slightly increase during the experiment. Hydrocarbons C18-C34 showed gradual fraction growth, and fraction of light hydrocarbons <C17 significantly reduced. Previously mentioned intervals C6-C17 and >C17 were used to split produced hydrocarbons on light oil and heavy oil pseudo-components to construct kinetic model of the organic matter transformation.   Liquid hydrocarbons yield also showed non-monotonic behavior with a maximum at the fifth cycle (25 h). This could be attributed to the overlay effect of hydrocarbons desorption, appearance of microfractures due to thermal expansion, continuous hydrocarbon release from non-isolated pores, and thermal cracking of kerogen, bitumen, and generated hydrocarbons. However, a substantial amount of hydrocarbons was obtained after rinsing the autoclave with toluene ( Figure 4). Results of the simulated distillation of the produced liquid hydrocarbons (Figure 4) showed the insignificant contribution of >C 35 fraction of hydrocarbons, which did, however, slightly increase during the experiment. Hydrocarbons C 18 -C 34 showed gradual fraction growth, and fraction of light hydrocarbons <C 17 significantly reduced. Previously mentioned intervals C 6 -C 17 and >C 17 were used to split produced hydrocarbons on light oil and heavy oil pseudo-components to construct kinetic model of the organic matter transformation.  Liquid hydrocarbons yield also showed non-monotonic behavior with a maximum at the fifth cycle (25 h). This could be attributed to the overlay effect of hydrocarbons desorption, appearance of microfractures due to thermal expansion, continuous hydrocarbon release from non-isolated pores, and thermal cracking of kerogen, bitumen, and generated hydrocarbons. However, a substantial amount of hydrocarbons was obtained after rinsing the autoclave with toluene ( Figure 4). Results of the simulated distillation of the produced liquid hydrocarbons (Figure 4) showed the insignificant contribution of >C35 fraction of hydrocarbons, which did, however, slightly increase during the experiment. Hydrocarbons C18-C34 showed gradual fraction growth, and fraction of light hydrocarbons <C17 significantly reduced. Previously mentioned intervals C6-C17 and >C17 were used to split produced hydrocarbons on light oil and heavy oil pseudo-components to construct kinetic model of the organic matter transformation.   SARA analysis of the produced liquid hydrocarbons was performed to understand the mechanism of the organic matter decomposition reactions during the experiment. The results showed that there was a significant decrease in the fraction of saturates (−7%) and a slight reduction in aromatic compounds (−3%) yield during the experiment ( Figure 5). As for polar compounds, the yield of resins and portion of asphaltenes soluble in dichloromethane:methanol mixture significantly increased (+8%), in contrast to the reduction in saturates and aromatics. The fraction of the remaining part of asphaltenes also slightly increased (+2%). The predominance of the aromatic compounds over saturates could be explained by the fact that almost all the samples are related to aromatic-rich type II and III kerogen (Figure 2). Tissot and Vandenbroucke showed that kerogen pyrolysates and natural extracts of these types of organic matter contain large amounts of aromatics [49,50]. As for the significant yield of polar compounds, a similar effect was observed by Behar in a comparative study of open-system and hydrous/anhydrous closed-system pyrolysis of lignite samples [32]. The authors explained the observed effect by accelerated primary cracking of kerogen compared with anhydrous conditions resulting in the generation of additional polar compounds. In our experiment, the strong solvent power of water possibly allowed these high-molecular-weight compounds to be produced. However, since the fluids were produced each five hours, the secondary cracking did not occur enough to transform resins and asphaltenes to low-molecular-weight compounds. Behar et al. studied thermal cracking of high-molecular-weight compounds. They performed C 14+ saturates, aromatics, and polar compounds closed-system pyrolysis to account for the role of secondary cracking during the maturation process [51]. Results showed that C 14+ aromatic and polar compounds generate a significant proportion of saturates and new heavy polyaromatics. It implies that corrected exposure time could improve the quality of the final liquid product. SARA analysis of the produced liquid hydrocarbons was performed to understand the mechanism of the organic matter decomposition reactions during the experiment. The results showed that there was a significant decrease in the fraction of saturates (−7%) and a slight reduction in aromatic compounds (−3%) yield during the experiment ( Figure 5). As for polar compounds, the yield of resins and portion of asphaltenes soluble in dichloromethane:methanol mixture significantly increased (+8%), in contrast to the reduction in saturates and aromatics. The fraction of the remaining part of asphaltenes also slightly increased (+2%). The predominance of the aromatic compounds over saturates could be explained by the fact that almost all the samples are related to aromatic-rich type II and III kerogen (Figure 2). Tissot and Vandenbroucke showed that kerogen pyrolysates and natural extracts of these types of organic matter contain large amounts of aromatics [49,50]. As for the significant yield of polar compounds, a similar effect was observed by Behar in a comparative study of open-system and hydrous/anhydrous closed-system pyrolysis of lignite samples [32]. The authors explained the observed effect by accelerated primary cracking of kerogen compared with anhydrous conditions resulting in the generation of additional polar compounds. In our experiment, the strong solvent power of water possibly allowed these high-molecular-weight compounds to be produced. However, since the fluids were produced each five hours, the secondary cracking did not occur enough to transform resins and asphaltenes to low-molecular-weight compounds. Behar et al. studied thermal cracking of high-molecular-weight compounds. They performed C14+ saturates, aromatics, and polar compounds closed-system pyrolysis to account for the role of secondary cracking during the maturation process [51]. Results showed that C14+ aromatic and polar compounds generate a significant proportion of saturates and new heavy polyaromatics. It implies that corrected exposure time could improve the quality of the final liquid product. X-ray computed microtomography of original and spent oil shale samples showed that extra pore space was generated by hydrothermal treatment (Figure 6). Microstructural properties of oil shale samples were characterized based on comparing the same sections of samples before and after the treatment. Samples BF13 and BF15 were scanned with 18 µm resolution. The porosity of the BF13 sample increased 3.3 times from 0.37% to 155.6 times from 0.0056% to 0.8716%, and specific pore surface area raised from 10.9 m 2 /m 3 to 897.8 m 2 /m 3 . The 3D-rendering of the pore space ( Figure 6) shows an interbedding effect, which could be explained by the observation that pores and microfractures are generated primarily in the organic-rich layers [52][53][54][55]. It should be noted that visible fractures appeared in almost all the samples. According to Kim et al., microfractures generally exist, and the cracks' dilation expands after the thermal treatment [54]. However, the alteration of microstructural properties is dependent on organic matter type because of pore structure connection with organic matter activation energy [56]. Saif et al. investigated the porosity evolution of source rock from room temperature to 500 • C using the dynamic X-ray synchrotron imaging technique. They showed that during step-like pyrolysis, first disconnected pores appear, then the poorly connected pore network transforms into a well-connected one accordingly with maturation and expulsion process, which was proved using TGA [57]. The maturation process also affects the surface roughness of the pores, which is demonstrated by means of multifractal analysis [58].
Energies 2021, 14, x FOR PEER REVIEW 9 of 16 1.23%, specific pore surface raised from 149.6 m 2 /m 3 to 423.4 m 2 /m 3 . Porosity of BF15 increased 155.6 times from 0.0056% to 0.8716%, and specific pore surface area raised from 10.9 m 2 /m 3 to 897.8 m 2 /m 3 . The 3D-rendering of the pore space ( Figure 6) shows an interbedding effect, which could be explained by the observation that pores and microfractures are generated primarily in the organic-rich layers [52][53][54][55]. It should be noted that visible fractures appeared in almost all the samples. According to Kim et al., microfractures generally exist, and the cracks' dilation expands after the thermal treatment [54]. However, the alteration of microstructural properties is dependent on organic matter type because of pore structure connection with organic matter activation energy [56]. Saif et al. investigated the porosity evolution of source rock from room temperature to 500 °C using the dynamic X-ray synchrotron imaging technique. They showed that during step-like pyrolysis, first disconnected pores appear, then the poorly connected pore network transforms into a well-connected one accordingly with maturation and expulsion process, which was proved using TGA [57]. The maturation process also affects the surface roughness of the pores, which is demonstrated by means of multifractal analysis [58].
The elastic properties were determined only for sample BF14, which retained its shape and did not fracture. Acoustic wave propagation velocities were measured, axial stress was applied parallel to the bedding plane before and after the experiment. Young's modulus and Poisson's ratio were calculated from the acoustic velocities. The results showed a 7% decrease in the modulus from 39.56 to 36.76 GPa and a 12% decrease in the Poisson's ratio, from 0.27 to 0.24, measured parallel to bedding. This behavior of elastic properties is expected for oil shale hydrous pyrolysis experiments [59][60][61]. However, it is shown in the literature that Young's modulus increases with maturity in natural conditions. The decrease in the macroscopic Young's modulus observed in the experiments is explained by transformation of kerogen to bitumen and heavy oil, which have lower Young's modulus (~2 GPa) and could stay trapped in the pore network before further transformation and expulsion of light oil and gas [62]. Ultimately, the elastic properties of oil shale strongly depend on the organic matter state in the rock. The elastic properties were determined only for sample BF14, which retained its shape and did not fracture. Acoustic wave propagation velocities were measured, axial stress was applied parallel to the bedding plane before and after the experiment. Young's modulus and Poisson's ratio were calculated from the acoustic velocities. The results showed a 7% decrease in the modulus from 39.56 to 36.76 GPa and a 12% decrease in the Poisson's ratio, from 0.27 to 0.24, measured parallel to bedding. This behavior of elastic properties is expected for oil shale hydrous pyrolysis experiments [59][60][61]. However, it is shown in the literature that Young's modulus increases with maturity in natural conditions. The decrease in the macroscopic Young's modulus observed in the experiments is explained by transformation of kerogen to bitumen and heavy oil, which have lower Young's modulus (~2 GPa) and could stay trapped in the pore network before further transformation and expulsion of light oil and gas [62]. Ultimately, the elastic properties of oil shale strongly depend on the organic matter state in the rock.
The results of the laboratory studies showed a significant change in the oil shale properties: organic matter and partial mineral decomposition, thermal expansion resulted in micropores and microfractures appearance. Microstructural and elastic properties were found to be in a strong correlation with geochemical properties and organic matter state. Produced hydrocarbons comprise a substantial fraction of polar compounds, which could be attributed to the high solvent power of subcritical water and lack of secondary cracking due to the short one-cycle exposure time (5 h), respectively.

Numerical Investigation: Optimization of the Kinetic Model
Bulk kinetic analysis was performed using HAWK workstation and Kinetics2015 software on the original and chloroform extracted oil shale samples from well D of the same reservoir. Three different heating rates (3, 10 and 30 • C/min) from 300 • C to 650 • C were used to obtain E a distribution. Discrete activation energy spectra with 1 kcal/mol spacing were obtained for fixed frequency factor A = 10 14 s −1 (Figure 7) to make the results of non-extracted and extracted samples pyrolysis comparable. The value of the frequency factor was chosen according to the observation that in most kerogen pyrolysis cases, the A value is between 3 10 13 s −1 and 8 10 14 s −1 , with a dominance near 10 14 s −1 . The results of the laboratory studies showed a significant change in the oil shale properties: organic matter and partial mineral decomposition, thermal expansion resulted in micropores and microfractures appearance. Microstructural and elastic properties were found to be in a strong correlation with geochemical properties and organic matter state. Produced hydrocarbons comprise a substantial fraction of polar compounds, which could be attributed to the high solvent power of subcritical water and lack of secondary cracking due to the short one-cycle exposure time (5 h), respectively.

Numerical Investigation: Optimization of the Kinetic Model
Bulk kinetic analysis was performed using HAWK workstation and Kinetics2015 software on the original and chloroform extracted oil shale samples from well D of the same reservoir. Three different heating rates (3, 10 and 30 °C/min) from 300 °C to 650 °C were used to obtain Ea distribution. Discrete activation energy spectra with 1 kcal/mol spacing were obtained for fixed frequency factor A = 10 14 s −1 (Figure 7) to make the results of non-extracted and extracted samples pyrolysis comparable. The value of the frequency factor was chosen according to the observation that in most kerogen pyrolysis cases, the A value is between 3 10 13 s −1 and 8 10 14 s −1 , with a dominance near 10 14 s −1 . The kinetic spectra were utilized to calculate logarithmic mean values of the activation energy for the original (bitumen + kerogen) and extracted (kerogen) samples. The following model was applied to determine the bitumen logarithmic mean value of the activation energy: 35 is a mass fraction of bitumen in S2 peak, = 0.65 is a mass fraction of kerogen in S2 peak, = 53.3 kcal/mol and = 53.1 kcal/mol were calculated from the spectra. Ultimately, the logarithmic mean value of the bitumen decomposition activation energy was determined from the equation: = 53.2 kcal/mol, which is extremely close to the kerogen activation energy. The obtained kinetic parameters from the open-system pyrolysis were used as an initial guess for further kinetic model optimization to derive relevant values of activation energy and frequency factor, that could be applied in fieldscale reservoir simulations using commercial software (CMG, Eclipse, etc.). The kinetic spectra were utilized to calculate logarithmic mean values of the activation energy for the original (bitumen + kerogen) and extracted (kerogen) samples. The following model was applied to determine the bitumen logarithmic mean value of the activation energy: where p bit = 0.35 is a mass fraction of bitumen in S 2 peak, p ker = 0.65 is a mass fraction of kerogen in S 2 peak, E LM a ker = 53.3 kcal/mol and E LM a = 53.1 kcal/mol were calculated from the spectra. Ultimately, the logarithmic mean value of the bitumen decomposition activation energy was determined from the equation: E LM a bit = 53.2 kcal/mol, which is extremely close to the kerogen activation energy. The obtained kinetic parameters from the open-system pyrolysis were used as an initial guess for further kinetic model optimization to derive relevant values of activation energy and frequency factor, that could be applied in field-scale reservoir simulations using commercial software (CMG, Eclipse, etc.).
In this study, a sequential reaction mechanism was utilized, where solid and immobile pseudo-components, kerogen and bitumen, crack to mobile and producible in the experiment heavy oil, light oil, hydrocarbon gas, and immobile char. This approach is widely used in the literature [36,63] and allows all real chemical species and reactions to be lumped into a few pseudo-components and pseudo-reactions. Adequate computational cost and moderate complexity make the model useful for petroleum reservoir simulations. Hydrocarbon gas cracking was neglected due to the high activation energy of the reaction and relatively low pyrolysis temperature of the experiment (Equations (1)-(5)): kerogen→a1 heavy oil + a2 light oil + a3 hydrocarbon gas + a4 char, (1) bitumen→b1 heavy oil + b2 light oil + b3 hydrocarbon gas + b4 char, heavy oil→c1 light oil + c2 hydrocarbon gas + c3 char, light oil→d1 hydrocarbon gas + d2 char.
Kerogen (K) and bitumen (B) masses were calculated cumulatively for all samples (BF1-BF15 in Table 1) using S 2 values of Rock-Eval pyrolysis of non-extracted and extracted oil shale samples before and after the experiment. Produced hydrocarbon gas (HCG) mass was determined using gas chromatography and gas meter data. Liquid hydrocarbon mass was weighted and divided into light (LO) and heavy oil (HO) using simulated distillation measurements (C 6 -C 17 for light oil, and >C 17 for heavy oil). Generated char mass was estimated as a remaining part of organic matter transformation. All mentioned masses are listed in Table 2.
The gas and liquid hydrocarbons produced after the final gradual pressure decrease to ambient conditions and the rinsing of the lines after the experiment were measured separately (row "after" in the left part of Table 2). These masses were further distributed among the nine cycles by multiplying each cycle hydrocarbon yield by (m 1-9 cum + m after )/m 1-9 cum for each of the group: HCG, LO, and HO. Such renormalization does not disturb the relationship between produced hydrocarbons at each cycle ( Figure 8) and allows the hydrocarbons to be accounted for, which left the rock but were produced only after final pressure release. The difference between Rock-Eval S 0 and S 1 peaks before and after the experiment (0.52 g and 1.8 g, respectively) was further deducted from the ultimate amount of the produced hydrocarbons in order to account for only hydrocarbons obtained due to pyrolysis reactions, not physical desorption process. However, such mechanisms also contribute to hydrocarbon production, and should be taken into account in the case of reservoir simulation. Under the assumption that all mentioned reactions belong to the first order type, the following system of ordinary differential equations (ODE) was obtained. The ODE system was solved numerically using Matlab software and assuming Arrhenius law dependence of reaction rate constants (Equations (5)-(11)).
Stoichiometry coefficients were calculated based on the mass conservation law and molecular weights of pseudo-components (Tables 3 and 4). The best match of the experimental data was reached by utilizing the least-squares approach based on the Levenberg-Marquardt algorithm. The objective function for this problem is the sum of squares of the differences between the ODE solution with parameters A, E a, and the experimental data (Equation (12)).
The optimized activation energy and frequency factors for Equations (6)-(11) are listed in Table 3. The results show (Figure 8) that the developed kinetic model describes organic matter evolution adequately. However, it should be noted that the short cycle durations (5 h) with further mobile hydrocarbons production allowed the secondary reaction contribution to be diminished and the mainly primary cracking regime of the experiment to be provided.
To improve the quality of the produced hydrocarbons, more attention should be paid to secondary reactions and their kinetic model development.   The results demonstrated in this study indicate that the evolution of oil shale com-   The results demonstrated in this study indicate that the evolution of oil shale composition during hydrothermal exposure could be adequately described using a number of pseudo-components (kerogen, bitumen, heavy oil, light oil, hydrocarbon gas, char) and a sequential reaction model. The relevant kinetic parameters were determined, especially for kerogen and bitumen decomposition. The obtained parameters might be applied for reservoir modeling of the in situ retorting process in the oil shale reservoir with comparable geochemical and petrophysical properties.

Conclusions
Cyclic subcritical water injection into the oil shale samples was studied as a promising in situ retorting technology for the first time. Artificial maturation noticeably alters geochemical parameters, elastic and microstructural properties of the oil shale. Rock-Eval pyrolysis of pre-and post-extraction oil shale samples before and after hydrothermal exposure showed bitumen and kerogen transformation to mobile hydrocarbons and immobile char. SARA analysis of the produced hydrocarbons showed a significant fraction of aromatic compounds due to the predominance of aromatic-rich type II and type III kerogen. Additionally, a large fraction of polar compounds was observed because of the high solvent power of subcritical water and lack of secondary cracking due to the short one-cycle exposure time. Hydrocarbons produced during the experiment were matched using nonlinear least-squares approach and a developed kinetic model, which adequately described kerogen and bitumen transformation with the production of liquid and gas hydrocarbons.
The X-ray computed microtomography derived porosity and specific surface area notably rose after the experiment due to thermal expansion and microfractures generation primarily in the organic-rich layers, kerogen and bitumen transformation to mobile hydrocarbons with further expulsion. It implies high potential efficiency of the cyclic subcritical water injection for the development of Bazhenov Formation and other oil shale reservoirs with similar geochemical and petrophysical properties.