Stability analysis of underground water-sealed oil storage caverns in China：A case study

The stability of underground water-sealed oil storage caverns is of great importance for safe excavation and operation. To analyze the scope of the failure zone and underground cavern stability accurately, a new method was developed that integrates the rock tunneling quality index Q-system and stability graph method with 3D laser scanning and numerical simulation. The point cloud data were obtained by 3D laser scanning, and the cavern model was built by using DIMINE software, which directly shows the 3D shape of the cavern. The rock mass physical and mechanical parameters and the corresponding stability coefficients were obtained based on Q-system and stability graph method. The plastic zone distribution and deformation characteristics of surrounding rock were analyzed through numerical simulation. Then, the corresponding relationship between caving zone and plastic zone was determined by comparing the numerical simulation results with the 3D laser scanning contour. The new method provides a reliable way to analyze the stability of the underground water-sealed oil storage cavern and also will helpful to design or optimize the subsequent support.


Introduction
Global energy shortages have prompted many countries to initiate or expand oil reserves. According to the proposal of the International Energy Organization, oil reserves of more than 90 days should be maintained. To ensure energy security, China needs to build underground oil reserves with an annual consumption of 10%-15% (Ding and Xie, 2006;Yang et al., 2000;Yuan et al., 2006). The use of underground water-sealed oil storage caverns in which oil and its products are reinjected into a selected underground space is a reasonable and important development trend for world oil reserves (Aberg, 1977;Goodall et al., 1988). Underground water-sealed storage cavern has the advantages of less land occupation, higher security, lower cost, and environmental benefits compared with traditional storage method, which is a widely used oil storage method . This approach is also recognized as a reserve bank with high global strategic security (State Council of the People's Republic of China, 2013). The basic principle of underground oil storage in unlined rock caverns is that the groundwater pressure surrounding the caverns should be higher than the storage cavern pressure, and the groundwater should flow into caverns to prevent oil leakages (Hassanpour et al., 2019;Lin et al., 2016). The water-sealing effect and rock mass stability are two key scientific problems for the construction of underground water-sealed oil storage caverns (Shi et al., 2018;Zhuang et al., 2017), especially the relation between rock mass stability and the caverns' stability.
Understanding the deformation characteristics of the caverns is particularly important for accurately analyzing the surrounding rock stability (Liu et al., 2020;Wang et al., 2017;Kar et al., 2017;Zhang et al., 2019). To further analyze the cavern stability, an appropriate technology must be adopted to monitor deformation and obtain basic information. Due to the complexity of in-situ condition, many monitoring methods, such as conventional electrical measurements, electromagnetic, ground penetrating radar (GPR) etc., were used to monitor the failure characteristics of underground caverns in the initial stage. However, the model accuracy only meets the basic engineering requirements and does not rigorously and accurately portray the true shape of the underground caverns (Chen et al., 2012).
3D laser scanning has been developed to monitor rock mass deformation in some engineering projects, such as metro tunnel (Andersson and Martin, 2009;Ghabraie et al., 2015;Shi et al., 2013;Vanneschi et al., 2014;Wang et al., 2009). The monitoring of tunnel-related deformation using 3D laser scanning has been studied from two aspects of data acquisition and processing (Xie et al., 2013). 3D laser scanning is also used to obtain the cross-section of caverns (Groccia et al., 2013). Compared with the traditional monitoring method, 3D laser scanning has the advantages of high efficiency, safety, accuracy, and abundant monitoring information.
Empirical methods are the widely used approach for the evaluation of rock mass stability. Based on the study of many underground excavation cases, Barton et al. (1974) proposed the rock tunneling quality index Q-system, which is one of the most common used empirical methods for stability analysis in mining and geotechnical engineering. Based on the work of Barton et al. (1974), Mathews et al. (1980) introduced the stability graph. Potvin (1988) revised the stability graph using case studies from underground Canadian mines. The stability graph is now widely used in underground hard rock mines. Although Q-system and stability graph methods are widely used, both fail to reflect the scope of the cavity failure zone in the stability analysis.
In this study, a new method, integrated empirical methods with 3D laser scanning and numerical simulation, was proposed and used in the underground water-sealed oil storage cavern in China. 3D laser scanning was used to monitor the cavern and obtain laser point clouds, combined the Q-system, and numerical simulation with the cavity model constructed by the data collected from 3D laser scanning, which shows the shape and depth of the damage and cavity of the underground water-sealed oil storage caverns.

Engineering background
The underground water-sealed oil storage cavern is located in Liaoning Province, the northeast of China, and which is in the northern temperate and semi-humid monsoon climate zones ( Figure 1). The geomorphology includes low and broad denudated hills with an altitude of 12.70-42.83 m, and mainly forest, farmland, and wasteland.
The design storage capacity of the project is 300 Â 10 4 m 3 for crude oil. Figure 2 (Lin et al., 2018) shows the layout of the underground facility, which is located more than 100 m below the surface. The underground facility is composed of eight storage caverns (two adjacent oil storage caverns are used as a group), namely North 1 -North 4 and South 1 -South 4, which are aligned parallel in the East-West direction. This project is equipped with two construction roadways that lead to oil storage caverns and water curtain roadways, which can be used as passes for large trackless equipment to transport slag. For each group of oil storage cavern, an oil inlet shaft with a diameter of 3 m and outlet shaft with a diameter of 6 m are set. Two oil storage caverns in each group of oil storage cavern are connected through tunnels to ensure that the oil is stored at the same height. During the cavern operation period, concrete sealing plugs are used to separate the connecting roadway between each group of caverns and tanks to form independent oil storage units. The crosssection of each cavern has a dimension of 19 m width and 24 m height ( Figure 3). The caverns are 934 m long along the East-West direction. The designed cross-section area is 436.1 m 2 . The elevation at the bottom of the cavern is À80 m. The study area in this paper labeled as the red boxed area in Figure 2.

Geology and hydrogeology
The surface layer of the cavern is residual soil. The rock mass in the reservoir area is mainly composed of medium coarse-grained granite, locally consists of diabase, fine-grained rock, hornblende porphyrite, granite porphyry, and other dykes. The slightly weathered and unweathered granite belong to hard rock, which is a whole massive structure. The joints  are generally in one to two groups, and the joint spacing is generally greater than 0.4 m, which is generally closed type. There are two types of weak zones in the deep unweathered and slightly weathered granite bodies. There is no apparent regularity of the two types with depth, which require attention during construction. There are no active or regional faults near the engineering site. Faults are exposed near the field area, but there is no indication that the Quaternary system faulted. There are no surface faults, which indicate minimal activity since the Quaternary and that the faults are inactive. The field geophysical exploration and drilling do not show any large-scale active faults that affect oil storage safety.
The underground water consists of Quaternary pore water and bedrock interstitial water. The underground water flow direction is generally southwest, and the water table is 13.33-38.31 m with an annual variation of 1-3 m. The hydrochemical type is sodium bicarbonate calcium chloride. The pH is between 6.68 and 9.15. Water salinity is 138.14-232.37 mg/L with good quality. The permeability coefficient is mostly less than 1.00 Â 10 À10 m/s and 1.55 Â 10 À9 to 3.50 Â 10 À7 m/s at the fractured zones with a mean value of 5.83 Â 10 À8 m/s.

Core statistics and rock mass quality classification
Borehole drilling is a commonly used geological exploration method in engineering. This approach can be used to record structural characteristics, such as joint groups, fillings, roughness, and alteration coefficient (Zhao et al., 2007). Drilling is carried out in the cavern to obtain more detailed information with regard to the geological conditions. Figure 4 shows part of the drill core images.
In addition to core provided by the construction site, we catalogued core samples and obtained a more detailed core data series. Q was calculated in the depth range of study according to the obtained data. The obtained Q was divided into five segments and the corresponding areas are marked in Figure 5.
The corresponding information was statistically analyzed according to the Q, and the proportion of Q in each section is shown in Figure 6.

Laboratory rock mechanics tests
Typical rock samples were collected at the construction site. Samples were taken at different positions, 200 m in front of the four caves including À50 to À60 m, À60 to À70 m, and À70 to À80 m sections. Intact rocks with minimal original cleats and fractures were selected. The physical and mechanical parameters of the rocks were determined from laboratory testing on intact rock samples following ISRM (1981) recommended methods and the test results are given in Table 1.

In-situ stress
According to the in-situ stress statistics (Table 2) of three boreholes (ZK8, ZK11, and ZK8), the maximum horizontal stress (r H ) is 6.19-11.50 MPa, the average dominant direction is NE74.3 , the minimum horizontal stress (r h ) is 3.63-9.02 MPa, and the vertical main stress (r V ) is 1.81-3.61 MPa. According to Engineering Rock Mass Classification Standard (GB 50218-94), R C =d max is generally >7.

Discontinuity survey and statistics
The investigation of discontinuity was mainly concentrated in the 501-515 area in the North 4 cavern. Detailed geological sketches and discontinuity statistics were carried out in the three zones selected for experiments to effectively determine the characteristics of discontinuities in the surrounding rock. Caverns of 15 m were measured and 121 discontinuities were recorded. The statistical data of discontinuity parameters are listed in Table 3.

3D laser scanning
Testing principle of 3D laser scanning 3D laser scanner was used to accurately and quickly scan the entire or part of the target object and directly obtain the point coordinates to generate the point cloud data (Dong et al., 2015;Ma, 2005). These data can reflect the 3D shape of measured objects (Abellan et al., 2006;Azarafza et al., 2019;Fardin et al., 2004;Han et al., 2013;Li et al., 2018).
The surveying principle of 3D laser scanning is that narrow laser pulses are emitted to the target through the scanning head of the 3D scanner, and the direction angle (a; h) of each laser beam is measured and controlled by an internal scanning control module. The time difference between laser transmitting and receiving is calculated by the internal measuring module, and the scanning center-to-target is calculated. The distance S between objects can also match the gray value of the target points according to the reflected laser intensity. The coordinate system used for data acquisition is a local coordinate system and the center of the    scanner is O. The X-and Y-axes are located on the horizontal plane of the local coordinate system. The X-axis is the scanning direction of the scanner, the Y-axis is the direction perpendicular to the scanning direction on the horizontal plane, and the Z-axis is the direction perpendicular to the scanning direction, as shown in Figure 7. The calculation formula (Feng et al., 2016) of 3D laser scanning target point p coordinates (X S , Y S , Z S ) is Field data collection of 3D laser scanning The scanning head of the laser range finder is used in the 3D laser scanning system. The laser scanner rotated 360 after reaching the cavern to continuously collect distance and angle data. After each 360 scan, the scanning head automatically raises its elevation according to the operator's pre-set angle for a new round of scanning, collecting data points on the larger rotating ring until the detection work is completed. The detection principle is shown in Figure 8. The specific methods of detecting caverns using 3D laser scanning are as follows.
1. Considering the factors of safety, equipment installation and convenience for measuring the coordinates of two points on the support, one or more equipment installation sites were selected to erect the detection equipment. For complex cavern structures that could not be scanned at one time, the general shape of the cavern should be analyzed, and several measuring points should be arranged according to the basic structure. 2. The scanning head was positioned in the cavern to be measured and should pay attention to avoid occlusion of the laser. 3. After the equipment was properly installed, a controller was used to carry on side detection process. 4. To incorporate the detection data into the coordinate system, the total station was used to measure the coordinates of the two points on the scanner head support and obtain basic data of the survey points. 5. Once scanning was complete, the data were transmitted to the computer and processed.

Data processing of 3D laser scanning
The original detection data obtained by 3D laser scanning are in "*.txt" format, which consists of angles and distance data. The original detection data were converted into "*. dtm" format to meet the requirement of DIMINE. The validity of the "*.dtm" file was verified in DIMINE and edited to import into SURPAC to generate and modify the entity model. After validation of the modified entity model, the three-dimensional cavern model was successfully constructed. Using the DIMINE software, a 3D cavern model was generated according to the detected data, as shown in Figure 9.
By using the 3D section interception function of DIMINE software, a 3D model of the cavern was intercepted with equal density sectioning to obtain a section cavern map ( Figure 10). The precision of the processing section is 0.1 m. In the effective scanning model (40 m), 4000 pieces of open area section maps were made with equal distances.
The 3D model can be de-noised to remove abnormal points caused by hardware conditions and field factors during the scanning process. Thus, the actual situation of the compound site and regional model required for the experiment was obtained. After slice processing was done in DIMINE, the sliced model was derived into "*.dwg" format and then imported into AUTOCAD (Figure 11).
According to Figure 11, the output cross-section diagrams were very irregular. Because there were many noise points and missing points caused by hardware and environment problems, a second processing was carried out by CAD to remove noise points and the missing points to form more complete isodensity cross-section diagrams. According to the test requirements, a three-segment 3D model of the cavern in Figure 10 was obtained according to the actual and experimental requirements of the compound site in Figure 12.
According to the actual engineering situation of the scanning zone, three parts of the scanning zone were selected to represent three different caving conditions, and the scanning data of these three sections were further processed. After denoising and supplementation in the later stage in AutoCAD, the cross-section contour of Figure 13 was obtained.
The midpoint data of Figure 13(a) to (f) were less incomplete. The point cloud density beyond 60 m of the scanner decreases owing to equipment reasons, so that a small  interference can greatly interfere with the accuracy of the point cloud data. Therefore, part of the section data in Figure 11 belongs to invalid point cloud data, and the obtained crosssection diagrams can no longer be completely restored to the cross-section diagrams of cavern and was abandoned. Finally, the isodense cross-section of Figure 14 and the three-dimensional model composed of isodense cross-sections were obtained.

Data analysis of 3D laser scanning
The repaired isodensity section of cavern was compounded with the design cavern outline as shown in Figure 3, and a compounded model (Figure 15) was obtained. The analysis area and caving scale were counted by designing the composite cross-section and scanning crosssections.
On the basis of the compounded model established in Figure 15(a), the falling area beyond the design model can be observed and further processing can be carried out to   obtain the falling area outline shown in Figure 16. Through this contour, the caving area of three areas can be identified.
The caving area can be obtained by taking the difference between the red area and blue area in Figure 16. Table 4 shows the specific values of the caving area.

Empirical methods
Rock tunneling quality index Q-system. The Q-system is essentially a weighting process in which the positive and negative aspects of the rock mass are assessed. The rock mass quality is  described using six parameters, which are the rock quality designation (RQD), joint set number (J n ), joint roughness (J r ), joint alteration (J a ), joint water reduction factor (J w ), and stress reduction factor (SRF) and Q was derived from the following formula (Barton et al., 1974) A description of Q is given in Table 5.  According to the evaluation index of Q-system, the corresponding caverns were investigated and evaluated. The Q values of three sections are summarized in Table 6.
Stability graph method. The stability graph method uses the hydraulic radius (HR) of the critical face and stability number (N 0 ) to estimate the stability of an underground opening (Hani et al., 2011). The HR is defined as the surface area divided by the perimeter for the vertical stope face. N 0 is defined as where Q 0 is the modified rock tunneling quality index Q with the stress reduction factor (SRF) and joint water reduction factor (J w ) equal to one, as these parameters accounted for separately within the analysis (Potvin, 1988).
The parameters A, B, and C are the stress factor, joint orientation factor, and gravity factor, respectively. The calculated parameters for all three sections are listed in Table 7.
The calculated N 0 and HR values are plotted in Figure 17. The zone A belongs to the transition zone but is close to the stable zone. Zone B is in the middle of the transition zone, which is a typical stable value. Zone C falls in the caving zone but is close to the boundary between the caving zone and transition zone.  According to the above stability numbers, small-scale caving is proposed in zone A. Zone B belongs to the transitional zone, but the caving scale belongs to the typical general Q area. For large-scale caving, the caving depth is not very deep and belongs to slight caving. Zone C belongs to a large-scale caving area.

2268
Energy Exploration & Exploitation 38 (6) Numerical analysis To verify the results of the empirical analyses, FLAC3D (Ding and Liu, 2018;Wang et al., 2019;Yang et al., 2019;Zhu et al., 2018) was used to calculate the deformation and depth of the plastic zones around caverns, owing to the staged excavation and rock support. Plane strain analysis was used because of the large ratio of length to cross-section dimension of the underground caverns. A single oil storage cavern was considered in this paper. The model is 100 m long, 100 m high and 24 m wide ( Figure 18). The rock mass was assumed to be an ideal elastic-plastic material. The generalized Mohr-Coulomb strength criterion was used to identify the elements undergoing yielding and the plastic behavior in rock masses. The rock mass was simplified as an equivalent isotropic medium. In this study, we use the Georgi method and Hoek-Brown method to determine the physical and mechanical parameters of the rock mass (Lin et al., 2011). The in-situ stress field was generated by applying stress on the right and back boundary faces of the model with the opposite faces fixed. The upper boundary was free, and all lateral boundary surfaces and the bottom face were foxed during the construction phase. The recommended values of physical and mechanical parameters of different regions of the rock mass were obtained after reduction (Table 8).
Vertical in-situ stress was generated to increase linearly with depth owing to the overburden weight. Rock stress measurements indicate that the maximum horizontal stress makes an angle of approximately 15 with the cavern axis. In this simulation, the maximum horizontal stress was assumed to be parallel to the cavern axis. The ratios of the maximum     and minimum horizontal stress to vertical stress were assumed to be 3.1 and 1.8, respectively. According to the experimental requirements and field investigation results, the model was divided into five parts along the Z-axis direction. The excavation area starts from Z ¼ 0, which was the predicted transition zone: Z ¼ 0-6; the predicted caving zone: Z ¼ 6-12; surrounding rock zone: Z ¼ 12-19; and stable zone: Z ¼ 19-24.
The main purpose was to simulate the mechanical state and change of rock mass during construction. Through simulation of the excavation process in the excavation area, the stress, displacement, and plasticity of nine sections in the stable area, transition area, and caving area selected by laser scanning were analyzed. The stress, displacement, and plasticity of each zone were analyzed according to the experimental requirements (Z¼ 0.5, 2.5, 4.5, 7, 9.5, 11.5, 19.5, 21.5, 23.5). Figures 19 to 21 show the distribution of stress, displacement, and plastic zone in each cross-section.
By simulating the stress, displacement, and plasticity distribution in different zones, the cross-section maps of the two sides and middle of each zone were selected for analysis and Table 9 was obtained. The stress, displacement, and plasticity distribution in different zones were preliminarily counted.
From the 3D laser scanning, it can be seen that the falling severity is C < A < B, and the most serious falling condition is in zone A. According to the above law, zones A, B, and C were all subject to tensile stress in the direction of the axial line of the cavity. The tensile stress in zone C was substantially greater than that in zones A and B. This likely owing to different scales of caving that occur in the latter, resulting in a stress relief effect in corresponding zones. The stresses in zones A and B were therefore relatively small.
Similarly, zone C had the smallest X-axis displacement. The strain in zones A and B was relatively similar, as well as their stress difference, so the stress and displacement in these zones were consistent with each other. The effect of stress relief and energy relief in the three zones was verified from another aspect.
According to the numerical simulation, plastic zone in the different zones was compounded and compared with the cavern scanning. The corresponding relationship between the falling zone and plastic zone was obtained by analyzing the compounded model. The plastic deformation and scanning section are shown in Figure 22.
From the comparison between the scanning results and numerical results, a certain relationship was observed to exist between the caving depth and plastic zone depth (Table 10).
Plastic zones A, B, and C (Table 10) fully conforms to the actual situation in the field. At the same time, on the basis of the scanning model and according to the numerical simulation, the plastic zone in the surrounding rock, which cannot be explored using scanning equipment in these three zones, can be inferred and used as a reference for the following support operation.

Conclusions
1. It is convenient to obtain more accurate and detailed geological data through detailed logging of geological boreholes in the study area, laboratory test, in-situ stress, and joint investigation, so as to provide a more detailed basis for the subsequent judgment of the collapse area and stability analysis of underground caverns. 2. Using 3D laser scanning to analyze the stability of the cavern, after de-noising and simplifying the collected point cloud data, using relevant software to carry out 3D model constructing, comparing the difference between the reconstructed section model and the designed section model, the collapse area of the specific position of the cavern can be obtained, so as to analyze the stability of the cavern. From the 3D laser scanning, it can be seen that the falling condition is C < A < B and the most serious falling condition is in zone A. 3D laser scanning is then used to carry out on-site measurements and followup processing to obtain the contour changes of excavated caverns. 3. According to empirical methods, small-scale caving is proposed in zone A. Zone B belongs to the transitional zone, but the caving scale belongs to the typical general Q area. By comparing the designed section model with the scanning model and simulation model, the relationship between caving depth and plastic zone depth is obtained, which provides a good reference for subsequent support work.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Project of NSFC Shandong United Fund (U1806208), National Key Research and Development Program (2016YFC0600803, 2018YFC0604401, 2018YFC0604604) and the Fundamental Research Funds for the Central Universities (N2001033).