Numerical Simulation and Experimental Studies of Karst Caves Collapse Mechanism in Fractured-Vuggy Reservoirs

Karst cave collapse usually occurs during the process of drilling and reservoir development. Karst cave collapse increases drilling risk and seriously a ﬀ ects oil production. In order to reveal the collapse mechanism of karst caves in fractured-vuggy reservoirs in Tahe oil ﬁ eld, rock mechanic experiments and numerical simulation studies were conducted. In rock mechanic experiments, three types of rock mechanics tests (uniaxial compression, triaxial compression, and Brazilian splitting) were carried out to obtain the basic rock mechanics parameters. In numerical simulation studies, the collapse conditions of a single karst cave and a cave system were studied. Numerical simulation results indicate that for a single karst cave, the size and geometrical features of karst caves are the main factors in ﬂ uencing collapse condition. For a cave system, a combination pattern and distance between caves are two main factors a ﬀ ecting collapse condition. In order to facilitate the application of these research results, the authors present formulas for calculating the critical conditions of karst cave collapse by means of multivariate linear regression method. The calculation formulas and prediction chart were veri ﬁ ed as consistent with actual ﬁ eld results. The research results in this paper have great theoretical and practical value for revealing the mechanism of karst cave collapse in fractured-vuggy carbonate reservoirs.


Introduction
The reserves of carbonate reservoirs account for half of the world's total oil reserves. Oil production of carbonate reservoirs accounts for 60% of the world's total oil production [1][2][3]. In China, carbonate reservoirs have been discovered in Bohai Bay, Ordos, Sichuan, and Tarim basins [4][5][6]. Tahe oil field served as the research target of this paper. It is located in the north of the Tarim basin. By 2012, the proven reserves of Tahe oil field have reached 12 × 10 8 t, and the annual oil production has reached 735 × 10 4 t. The fractured-vuggy reservoir is an important type of carbonate reservoir [7,8]. Karst caves play an important role in fractured-vuggy reservoirs as the main storage spaces [9]. The existence of karst caves usually leads to drilling tools being dropped and drilling fluids being lost. Furthermore, karst caves are likely to collapse because of changing wellbore pressure from fracturedvuggy reservoir development [10]. Cave collapse seriously affects oil well production [11]. Therefore, studying the collapse mechanism of karst cave is of great theoretical and practical value.
In order to study, the mechanism of cave collapse, the formation, evolution, and distribution of karst caves should be first understood. Previous scholars have performed some research on Tahe karst caves. The overall burial depth of karst caves in Tahe oil field is from 4,300 to 5,800 m. Tahe oil field undergoes the transformation of karstification in Caledonian and early Hercynian [12]. Because of multistage karstification, carbonate reservoirs develop multiple sets of cave systems. The size of caves ranges from a few meters to tens of meters [13]. Xiao et al. [14] found that the evolution of karst caves generally moves from cave to cave layer to cave system. Xu et al. [15] statistically sorted out the depths of 323 caves in Tahe oil field. Results showed 28% of karst caves being distributed within 20 m below the unconformity surface. The number of karst caves encountered in a single well is usually one to three. Loucks' cave-system distribution model shows vertical and horizontal as the two main distribution modes [16]. Vertical karst caves are mainly multilayered falling-water caves. Horizontal karst caves are widely developed along karst channels because of the influence of underground channels and permeable zones. Li et al. [17] studied the distribution law of karst caves in the fourth block of Tahe oil field. They proposed that cave-system spatial distribution could be divided into horizontal and vertical states. They divided the vertical cave system into three karst zones, namely, the surface karst zone, the vertical percolation karst zone, and the subsurface karst zone. Generally speaking, karst caves in Tahe oil field mainly exist in the form of a cave system and are distributed near the unconformity.
Karst cave collapse directly affects drilling safety and oil production. The critical nature of the potential for karst cave        3 Geofluids collapse needs to be discussed urgently. Lucia et al. [18][19][20][21][22] discussed the collapse condition of karst caves. They thought that given the influence of overburden pressure, the cave roof is subjected to tensile stress, and the cave wall is subjected to shear stress. With an increase in tensile stress and shear stress, cracks would appear on the cave roof and cave wall, eventually leading to cave collapse. Tang [23] also pointed out that karst cave collapse is caused mainly by cracks. Qian et al. [24] studied cave collapse in the central and northern Tarim basin. They found that the collapse time of karst caves can be divided into early collapse and late collapse. All of the above understandings represent qualitative knowledge of karst cave collapse. The critical condition of karst cave collapse is seldom described quantitatively. Zheng et al. [25] studied the mechanism of cave collapse and deduced the formula for calculating critical conditions of cave collapse using rock mechanics methodology. They assumed the geometrical feature of a cave to be square. As the overlying strata continue to deposit, the burial load goes beyond the rock compressive strength. Then, the cave begins to collapse. There are three independent variables in the formula including rock compressive strength, cave span, and distance away from unconformity surface. The value of rock compressive strength is replaced by the average compressive strength of limestone. This model can be used to judge whether a karst cave has collapsed. However, this model has some shortcomings. First of all, according to geological knowledge, there are many different shapes of karst caves deep underground. It is inappropriate to use a square to indicate geometrical feature of a cave in this model. Secondly, the mechanic parameters of carbonate rock in different areas are obviously different. Using an average value for the mechanic parameters of limestone is inappropriate for calculating the critical condition of karst cave collapse.
Presently, the critical condition of karst cave collapse can be studied by physical model testing and numerical simulation. The factors influencing the collapse condition are geometrical features, geologic characteristics, combination pattern, rock mechanics parameters, and so on. The mentioned variables present great challenges for physical model testing. Physical model testing consumes a lot of time when adjusting parameters. It is not a practical method for wide use. However, numerical simulation can solve the problem effectively [26][27][28]. Numerical simulation can analyze the collapse condition of a single cave, as well as a cave system. There are many numerical simulation methods for the stress field, such as the boundary element method, finite element method (FEM), discrete element method, and finite difference method [29]. Among these numerical simulation methods, FEM is the most widely used. Wang et al. [30] assembled numerical models for gas and oil storage caverns using FEM and discussed the special mechanical characteristics of bedded salt mass. The damage potential of the mudstone interbeds and rock salt surrounding the gas storage cavern were presented. Wang et al. [31] built a three-dimensional (3D) geomechanical model of a salt cavern in bedded rock salt formations and calculated the critical collapse length. Wang et al. [32] built a 3D geomechanical model of JK-A cavern to find the cause of cavern collapse. They showed that large-span flat roof and rapid drop of the internal gas pressure were the two main causes of cavern roof failure. The numerical simulation method, based on FEM, has been used widely to find out the cause of cavern collapses. As a result, this work applied FEM to study the collapse mechanism of karst caves.
In this study, with the assistance of rock mechanics experiments and numerical simulation methods, three innovative accomplishments were achieved. First, single cave models with different geometrical features were established, and the stress variation of a single cave was simulated. The influence of roof thickness, cave span, and burial depth on collapse condition was analyzed and discussed here. Next, numerical models of a cave system were designed for the first time. Two influencing factors, combination pattern, and distance between caves were analyzed and discussed here. The differences in the collapse condition between a single cave and a cave system were explored. Finally, formulas were developed for quantitatively calculating the stress of a single cave and cave system. They are presented here, as is a chart for judging critical collapse condition. The formulas and charts can be applied in Tahe oil field to improve the success rate of drilling and increase oil production.

Experimental Section
2.1. Materials. In order to obtain the rock mechanics parameters of the carbonate rock in Tahe oil field, nine cores were obtained from eight wells: well T314, well T607, well T615, well S86, well Tashen 2, well TK427, well S72-11, and well S80. Two samples were collected from well S80, and one was collected from each of the others. Photos for core samples are shown in Figure 1.
The core samples were observed carefully and described, and their lithology and burial depth were recorded. The collected samples were cut and polished. Six standard specimens with a diameter of 50 mm and a height of 100 mm were included, as well as three standard specimens with a diameter of 50 mm and a height of 50 mm. The specific information is shown in Table 1.

Experimental Procedures.
In this study, rock mechanics testing consisted of three tests-uniaxial compression test, triaxial compression test, and Brazilian splitting test. Nine specimens were numbered and divided into three groups. The uniaxial compression test was carried out on specimens 1 to 3, the triaxial compression test on 4 to 6, and the Brazilian splitting test on 7 to 9.
For the uniaxial compression test, the instrument used was the GAW-2000 servo rock rigid pressure tester. The pressure of rock specimens was measured by the pressure sensor of the tester, and the compression deformation was 4 Geofluids measured by a static resistance strain gauge. The specific test steps were as follows: (a) put the specimen in the center of the bearing plate with a spherical seat of the press. (b) Adjust the bearing plate to make the specimen uniformly loaded. (c) Add axial pressure at a speed of 0.5 MPa/s until the specimen is destroyed. The triaxle compression test was carried out with the RLW-1000 rock automatic triaxle servo rheometer. The test steps were as follows: (a) install the resistors in the middle of the specimen. Fit each specimen with two resistors in the axial and circumferential directions. (b) Place the specimen on the press and connect it to the data acquisition device in the way of the full bridge. (c) Apply axial pressure and confining pressure simultaneously to the specimen at a speed of 0.05 MPa/s. When confining pressure reaches 10, 20, and 30 MPa, the confining pressure remains stable. Add axial pressure at a speed of 0.7 MPa/s until the specimen is destroyed. (d) Record the axial load and rock failure form.
The Brazilian splitting test was carried out with the WDW-100E universal testing machine. The test procedures were as follows: (a) put the specimen between the upper and lower bearing plates of the machine. (b) Adjust the bearing plate to make the specimen uniformly loaded. (c) Add pressure at a rate of 0.5 MPa/s until the specimen is destroyed. (d) Record the maximum failure load and observe the failure process. According to the elastic-plastic characteristics of carbonate rock, the classic Mohr-Coulomb constitutive model was adopted in this analysis. The yield surface of the Mohr-Coulomb model is a hexagon in the deviatoric stress plane [33], as shown in Figure 2.

Numerical Simulation Model
The Mohr-Coulomb model yield criterion assumes that damage occurs when the shear stress at any point in the rock reaches a certain value. The model is based on the maximum principal shear stress yielding theory. For the σ~τ coordinate system, the failure line is shown in Figure 3 [34].
The governing equation for the Mohr-Coulomb failure criterion is shown as: In the Mohr-Coulomb failure model, τ and σ can be replaced: After transformation, the equation becomes where σ 1 , σ 2 , and σ 3 represent the first, second, and third stress, respectively; φ is internal friction angle; c is cohesive force.

Initial
Conditions of the Model. Figure 4 shows the numerical model and boundary conditions of the karst cave. Vertical stress (P v ) was applied to the top of the model. Horizontal stress (P H ) was applied to the surrounding boundary of the model. The bottom boundary of the model was fixed. Vertical stress was caused by the gravity of the overlying strata. P v was replaced by the value multiplied by the average bulk density of the strata and the actual burial depth of the cave. The model size of the cave was consistent with the actual size. Considering the boundary effect, the distance between the cave and the boundary was more than five times the cave span. The overall size of the model was 100 × 100 m.
Under the condition of high geostress, the deeply buried rock was in a state of compression. There were two analysis steps in the simulation process: setting geostress balance in the original state and loading the rock under compression according to the change in burial depth. In the calculation process, horizontal stress gradient was T h = 0:0155 MPa/m, vertical stress gradient was T v = 0:025 MPa/m, lateral pressure   Step time = 0. 16 Step time = 0.32 Step time = 0.48 Step time = 0.64 Step time = 0.80 Step time = 1.00 (a) Cubic cave S, Mises (avg: 75%) 89.73E+06 Step time = 0. 16 Step time = 0.32 Step time = 0.48 Step time = 1.00 Step time = 0.80 Step   In order to gain a better understanding of the mechanism of cave collapse in fractured-vuggy reservoirs, numerical simulation studies were conducted on the basis of the experimental results. The whole numerical simulation was divided into three parts. In the first part, three different shapes of karst cave were designed to simulate the stress variation process. In the second part, a large number of models of a single cave were designed, and the factors affecting the critical conditions of cave collapse were analyzed. In the third part, 14 models of cave system were designed, and the differences of critical collapse conditions between a single cave and cave system were analyzed. Detailed information on numerical simulation models is shown in Table 2.

Experimental Results.
In the experimental section of this work, rock mechanics tests were carried out. Compressive strength, elastic modulus, and Poisson's ratio parameters of the rock were measured by the uniaxial compression test, and the calculation formulas are shown in Equation (5)-(7). The results of uniaxial compression testing are shown in Table 3. Through the triaxle compression test, compressive strength and elastic modulus of the rock under different confining pressures were obtained, as well as cohesive force and internal friction angle. The results are shown in Table 4. The tensile strength of rock was obtained by the Brazilian splitting test. The calculation was carried out according to Equation (8).
The results are shown in Table 5.
where σc is rock compressive strength; P is maximum failure load; A is cross-sectional area perpendicular to the loading direction; E50 is elastic modulus; σ50 is the value of stress equivalent to 50% of the compressive strength; ε h50 is the longitudinal strain value when the stress is 50% of the compressive strength; μ is Poisson's ratio; εd50 is the transverse strain value when the stress is 50% of the compressive strength; σt is tensile strength; D is specimen diameter; H is specimen height.
The results of the above three tests were sorted out, and the mechanics parameters, such as compressive strength, elastic modulus, Poisson's ratio, cohesive force, and internal friction angle of the carbonate rocks in Tahe oil field, were obtained. The basic mechanics parameters were used in subsequent numerical simulation. Rock mechanics parameters are shown in Table 6.

Stress Variation Process of a Single Cave.
In the numerical simulation, geometric features of the karst cave were processed abstractly. This work simplified them into three shapes: cube, sphere, and cylinder. The karst cave model was built to analyze the stress distribution of cave surrounding rock. In simulation, roof thickness of the karst cave was 20 m, cave span was 20 m, and burial depth was 6,000 m. The model was established according to Figure 4. The stress variation processes of cubic cave, spherical cave, and cylindrical cave are shown in Figure 5.
The numerical simulation method can calculate the whole process of three shapes of karst cave from initial condition to final collapse. As shown in Figure 5, collapse position was obtained.
For the cubic cave, the stress concentration first appeared around the cave wall and then extended to the roof and bottom of the cubic cave. The maximum stress emerged in the upper and lower parts of the cave wall. Finally, the cubic cave moved downward. The cave roof was concave, and the cave wall bent inward to a certain extent. When the step time was 0.64, Mises stress reached 78.08 MPa, which is greater than 74.2 MPa. It was concluded that the cubic cave began to collapse. The influence scope of stress was two times the cave span in the horizontal direction and 4.2 times the cave span in the vertical direction.
For the spherical cave, the stress variation process was similar to that of the cubic cave. The stress first concentrated on the middle part of the sphere. Then, the area of stress concentration expanded to the upper and lower parts of the cave. Upon reaching a stable state, the cave compressed into an ellipsoid. When the step time was 0.80, Mises stress reached 74.77 MPa. From the stress variation diagram, the most vulnerable position at which to collapse can be seen as the middle of the sphere. The stress influence scope was two times the cave span horizontally and four times the cave span vertically.
The cylindrical cave was different from the cubic cave and the spherical cave in the stress field. The stress first appeared on the surrounding cave wall. Then, the stress expanded to the two ends of the cave. When the step time was 0.64, Mises stress reached 74.72 MPa. Finally, the stresses at the upper and lower parts of the cylindrical cave wall were the largest. The two ends of the cave were shown to be the most unstable position. The stress influence scope was 2.9 times as long as the cave span horizontally and 3.9 times vertically.

Influencing Factors on Cave
Collapse. In this part of the paper, the influences of roof thickness, cave span, and burial depth on karst cave collapse condition are discussed. Considering that the overall burial depth of the karst caves in Tahe oil field is between 4,300 and 5,800 m, close to the unconformity surface, the calculation conditions were as follows: cave span was 2.5 to 30 m, roof thickness was 5 to 30 m, and burial depth was 4,000 to 6,000 m. The calculation results of the three karst cave shapes are shown in Figure 6.
In view of the results of the numerical simulation, a good linear relationship was found between Mises stress and cave span or roof thickness (see Figure 6). The critical condition of cave collapse could be characterized by Mises stress value. When roof thickness was constant, Mises stress of a karst cave increased with an increase in cave span. That is to say, the larger the cave span, the more likely the cave will collapse. When the cave span was constant, Mises stress of a karst cave decreased with an increase in roof thickness. That is, the smaller the roof thickness, the easier the cave will collapse.

Geofluids
When roof thickness and cave span were fixed, Mises stress of a karst cave increased with an increase in burial depth. That is, the deeper the cave, the more likely the cave will collapse. When cave span, roof thickness, and burial depth were the same, the Mises stress of a cubic cave was the largest, that of a cylindrical cave was the second, and that of a spherical cave was the smallest. After comparing and analyzing the calculation results of the three cavity shapes, the cubic cave was found to be the most likely to collapse, and the cylindrical cave and the spherical cave were found to be harder to collapse.

Analysis of Collapse Conditions of a Karst Cave System.
Karst caves deep underground usually exist in the form of a cave system. Fourteen models of cave system were designed. The calculation conditions were as follows: cave span was 10 m, roof thickness was 20 m, and burial depth of the cave was 5,000 m.
Among the all models, nine represented double caves. The combination patterns of double caves were horizontal, vertical, and oblique 45°. The Mises stress values were calculated as the distance between caves increased. The simulation results are shown in Figure 7.
As shown in Figure 7, Mises stress decreased with an increase in distance between caves. When the distance was larger than 1.5 times the cave span, the cave system basically reached a stable state. Among the three combination patterns, the Mises stress value of the horizontal cave system was the largest, followed by the vertical cave system, and the stress of caves combined obliquely at 45°was the smallest. Among the three types of karst cave, Mises stresses of the cubic cave and spherical cave were the largest, stresses of cylinder, and sphere were the second, and those of cube and cylinder were the smallest.
On the basis of double caves, this work further simulated the critical collapse condition of multiple caves. According to the geological background, five models, including two threecave models and three four-cave models, were designed. By changing the distance between caves, Mises stress values of cube, sphere, and cylinder cave system were calculated, respectively. The calculation results are shown in Figure 8.
The effects of combination pattern and distance between caves on critical collapse conditions were studied. As shown in Figure 8(a), when the three caves were combined horizontally, Mises stress increased rapidly with a decrease in distance. When the distance was about 15 m, the stress value was equal to that of a single karst cave. The smaller the distance, the greater the Mises stress, and the more likely the cave system to collapse. When the three caves were combined vertically (see Figure 8   13 Geofluids can be concluded that distance has little effect on the stability of a cave system when caves are combined vertically. As shown in Figures 8(c) and 8(d), the variational tendency of the stress curve was similar. When the distance was less than 10 m, the stress increased rapidly. When the distance was more than 10 m, the stress gradually began to stabilize. But the difference is that the stress value of (d) was smaller than that of (c) on the whole. In view of the last combination pattern (see Figure 8(e)), the cave system remained stable when the distance was greater than 5 m. The stress value was slightly larger than that of the vertical cave system.
According to the analysis of the above combination patterns, results showed a distance of 10 m to be a critical value. When the distance went beyond 10 m, Mises stress of a cave system did not change dramatically. The critical value was consistent with the stress influence scope of a single cave. On the other hand, when there were multiple caves combined in the vertical direction, the value of Mises stress decreased. When the cave systems were combined horizontally, the stresses influenced each other, and the stress value was greater than that of a single cave. By comparing the simulation results of a single cave and cave system, a significant difference was found in the collapse condition. For a single cave, a cylindrical cave was found easier to collapse, and a cubic cave and spherical cave were found harder to collapse. Correspondingly, a cubic cave system was the most likely to collapse, while a cylindrical cave system and spherical cave system were found harder to collapse.

Cave Collapse Prediction Study.
Numerical simulation was conducted to study the critical condition of cave collapse in fractured-vuggy reservoirs. According to the experimental and numerical simulation results, calculation formulas and a prediction chart are proposed.   and burial depth. Therefore, it can be assumed that the linear relation between Mises stress and burial depth, roof thickness, and cave span is as follows: where X 1 , X 2 , X 3 , and X 4 are undetermined coefficients, and F, D, H, and L, respectively, represent Mises stress, burial depth, roof thickness, and cave span.
According to the results of the numerical calculation, the calculation formulas of cave collapse were obtained through the multivariate linear regression method as follows: Cubic cave: The results of numerical calculation were compared with those of formula calculation. The results are shown in Figure 9.
In the study, data were summarized further for a single karst cave. It was found that the ratio of roof thickness to cave span can be used as a criterion for cave collapse at different depths. In order to show this result more intuitively and quantitatively, the results of the numerical calculation were put into a coordinate system, with the horizontal axis being the ratio of roof thickness of cave span and the vertical axis being the burial depth of the cave. As shown in Figure 10, the stress value reaching the rock compressive strength is expressed in red, while the stress value under compressive strength is expressed in blue. Through analysis, a good mathematical correspondence was found. Red dots are all located in the left side, with blue dots on the right. Therefore, the left side of the chart represents the collapsed area, and the right side represents the uncollapsed area. This chart can be used to judge quantitatively whether a karst caves have collapsed or not.

Collapse Prediction of a Cave
System. On the basis of a single karst cave, this work also performed linear regression analysis for a cave system. Five factors affect the critical condition of cave system collapse, including burial depth, roof thickness, cave span, distance, and combination pattern. The combination pattern of a cave system cannot be expressed by numerical value, so the numbers 1 to 5 were used for calculation. Numbers 1 to 5 represent the five combination patterns in Figure 8, respectively. When establishing the multifactor stress prediction formula for a cave system, a linear relationship was assumed between Mises stress and cave span, roof thickness, burial depth, cave distance, and combination pattern, as follows: where X 1 , X 2 , X 3 , X 4 , X 5 , and X 6 are undetermined coefficients, X is the distance between caves, and C is the combination pattern of karst caves. According to the results of the numerical simulation, the optimum solution was obtained by multivariate linear regression method. The predictive formulas and coefficient of correlation can be written as follows: Cubic cave system: F = 0:010D + 0:36H + 0:75L − 0:359 X − 1:52C + 13:10 , R = 0:882. The results of the numerical calculation were compared with those of formula calculation. The results are shown in Figure 11.

Engineering Verification.
In order to verify the accuracy of the formula and judgment chart for a single karst cave, six wells in Tahe oil field were selected from the literature [22]. Mises stress was calculated by using a single-cave prediction formula in Section 4.3. The calculation results of six wells were compared with the actual measured results in the project, as shown in Table 7.
The values of the actual burial depth and H/L are shown on the judgment chart ( Figure 12). In the chart, well TK 1 to TK 4 fell into the uncollapsed area, and well TK 5 and TK 6 fell into the collapsed area.
In order to verify the accuracy of the prediction formulas for the cave system, well TK481X in the fourth block of Tahe oil field was selected for analysis. Well TK481X was drilled in Ordovician formation at 5,503 m. Unfilled karst caves were  It can be seen that the numerical simulation method and model adopted in this paper are reasonable and reliable.

Conclusions
In the paper, rock mechanics tests and numerical simulation method were conducted to study the mechanism of karst cave collapse. The following conclusions can be drawn.
(i) Rock mechanics testing can obtain rock mechanics parameters in Tahe oil field. Compressive strength is 74.2 MPa, elastic modulus is 36.3 GPa, Poisson's ratio is 0.25, cohesive force is 12 MPa, and internal friction angle is 36°( ii) Numerical simulation can establish and calculate a large number of karst cave models. The critical collapse condition of a single cave and cave system can be proposed as follows (a) For a single karst cave, the larger the burial depth, the larger the cave span, and the smaller the roof thickness, the more likely the cave is to collapse. Geometrical features of a karst cave affect collapse conditions. A cylindrical cave is easier to collapse, and a cubic cave and spherical cave are harder to collapse. The three types of karst cave are roughly the same in the influence scope of Mises stress. The influence scope of stress is about 2 times the cave span in the horizontal direction and about 4.2 times the cave span in the vertical direction