Correction: Groundwater regime evolution in an excavation damaged zone of underground water-sealed petroleum storage caverns

Energy is getting increasingly serious with fast economic development. Petroleum storage is an effective approach to alleviate or avoid the influence of oil crisis on the economy, society and security, through balancing the demand and supply, stabilizing the oil price and strategic oil reserve for unexpected emergencies. At present, many countries have built petroleum strategic reserves system with a large and increasing amount of oil storage.Petroleum storage includes conventional storage on the ground and the emerging underground storage. Because of the large amount of storage with very little land usage, the underground storage plays more and more important role in oil reserves. Because of oil's lower density relative to water and its poor intermiscibility with water, underground water-sealed petroleum storage caverns have usually been built at a certain depth below the phreatic line, in which the oil can be sealed by underground water. Those deep rock caverns are the natural oil reservoirs supported by the rock bearing capacity. Compared with other oil reserve approaches, underground water-sealed petroleum storage caverns have many advantages, such as good safety performance, few lands used, low investment, long usage period and little environmental pollution (Zhang et al., 2019;Qian et al., 2025).In 1948, the first underground water-sealed petroleum storage caverns were built in Sweden. From then on, more and more underground water-sealed petroleum storage caverns have been developed and applied in other places of the world. At present, more than 300 underground water-sealed petroleum storages have been built in more than 30 countries, which are mainly located in Scandinavian Peninsula, South Korea, Japan, Germany, France, and Saudi Arabia (Chung et al., 2009;Xue et al., 2015). In China, underground water-sealed petroleum storage caverns project is still in its infancy and requires further comprehensive development.The excavation of underground rock caverns is accompanied by unloading of surrounding stresses and release of energy. The surrounding rock is unloaded in the direction normal to the cavern surface, and the EDZ forms gradually induced by the rebound deformation on the excavation boundary of surrounding rocks. In EDZ, the mechanical and hydraulic properties of rocks greatly change as a result of rock damage and disturbance. The EDZ controls the cavern stability. As a result of the external loading and unloading, the EDZ evolves, and the permeability coefficient in the EDZ grows sharply, which greatly endanger the stability of the rock cavern.Consequently, understanding the evolution of the mechanical and hydraulic properties of EDZ is of great significance in engineering practice, with respect to supporting optimization, safe construction and stability evaluation.From early 1980s, EDZ has received a lot of attention, especially for the safety evaluation of nuclear waste treatment plants. Many theoretical and experimental works have been conducted, in particular the classification of the EDZ. In 1982, experimental and monitoring study on the EDZ was carried out in the Underground Research Laboratory of Atomic Energy of Canada Limited (Read, 2004). Autio et al.(Autio et al., 1998) classified the EDZ into the crushed zone, the fractured zone and the fissured zone. Backblom et al. (Bäckblom and Martin, 1999) found the EDZ comprised failed zone, damage zone and disturbed zone. Sato et al. (Sato et al., 2000) conducted field tests about EDZ and classified the EDZ into fracturing zone, desaturation zone and stress redistribution zone. Sheng et al.(Sheng et al., 2002) concluded that there exist damage zone, affected zone and slightly affected zone in the EDZ. Maritino et al. (Martino and Chandler, 2004) believed that, after excavation, the surrounding rock consisted of excavation damage zone and excavation disturbed zone. Cai et al. (Cai and Kaiser, 2005) proposed quantitative evaluation of the degree and the range of EDZ with the micro-mechanical model, using parameters of crack density and crack distribution. Tsang et al.(Tsang et al., 2005) summarized some existing findings and proposed the concepts of the excavation influence zone and excavation damage zone based on water flow and seepage properties. The seepage flow state and its properties changed dramatically in the excavation damage zone, but not in the excavation influence zone. Based on the mesomechanics, Souley et al. (Souley et al., 2001) studied the evolution of damage variable and the seepage property in the EDZ. Kwon et al. (Kwon et al., 2006) investigated the engineering characteristic of the EDZ in an engineering project. Li et al.(Li et al., 2022;Qiao et al., 2022) developed an identification method for dominant seepage channels in fractured rock masses based on water curtain system tests and discrete element modeling, established discriminant criteria. Cao et al. established the multi-parameter nonlinear dynamic seepage equation of fractured rock masses and studied the structural stability of the seepage system in fractured rock masses (Cao et al., 2024;Jia et al., 2024;Cao et al., 2025). In addition, evolution of seepage properties of damaged rock has been extensively studied (Salari et al., 2004;Blümling et al., 2007;Levasseur et al., 2013;Deshpande et al., 2025;Sun et al., 2025).Existing studies have predominantly focused on the qualitative classification and static characterization of the EDZ, with limited attention paid to the dynamic coupling mechanisms between damage evolution and transient groundwater seepage during construction. Moreover, the quantitative relationship between damage variables and permeability evolution in the EDZ has not been systematically established, and coupled numerical models suitable for large-scale storage projects remain underdeveloped. Furthermore, the validation of existing models against real-world engineering data is insufficient, constraining their reliable application in practice. To address these gaps, in this paper, three-dimensional unstable seepage evolution with excavation induced damage surrounding underground water-sealed petroleum storage caverns was studied numerically, considering the initial permeability characteristics in surrounding rock and the rainfall process. The damage constitutive model of EDZ was also established to investigate its damage evolution during cavern excavation. The method and findings were applied to the large-scale on-going Huangdao underground water-sealed petroleum caverns storage project taking a project of underground water-sealed petroleum caverns storage as a case study. The seepage evolution law and its correlation with damage properties in the EDZ were determined. Comparison between numerical modeling results and field measurements has been performed to verify the numerical method. In addition to academic merits in rock mechanics, hydro-mechanical coupling and underground excavation, the findings are also of engineering significance to real-timely adjust the construction period and construction sequence during cavern construction and to effectively design the water-sealed parameters after construction. The objectives of this paper are to establish an efficient numerical method to investigate the seepage properties in the EDZ and its correlation with damage evolution during and after the excavation of water-sealed petroleum caverns, and to verify its applicability and reliability through comparison with field measurements of an on-going engineering project.As the national oil strategic reserve project, the project of underground water-sealed petroleum caverns storage studied in this paper is the first extra large scale project in domestically at present, which is 300×10 4 m 3 of design capacity. The location of caverns storage is very superior with the convenient traffic. It locates at the north part of Tailong region in south of Jiaodong Peninsula with the hills landform. The caverns are in Longque Mountain trending from the east to the west, with the altitude of 280 ~350m, There is the escarpment on the north side of ridge and the abrupt slope on the south side, which is about 35~55°. The gully on the mountain foot extends on north-south direction or on north-east direction. The main part of caverns storage locates on the southwest part of Longque Mountain, with the average altitude of 220m. The highest point is called Dadingzi, which is 350.9m altitude. The lowest point is at the position of mouth of the well of zk012 drill hole, which is 97.50m altitude. The surface of caverns storage location is covered by large area of forest, including pine tree and locust tree, and by certain area of exposed rock, and by a few area of farmland near the mouth of the well at the mountain foot.The area that caverns storage located is on joint part of two plates, the north China plate and the Yangtze plate, which is the Jiaonan-Weihai orogenic belt with ductile shear zone, brittle faulting structure as main characteristic and the poorly developed fold structure. It also locates on the south rim of Mouping-Jimo fault zone trending mainly from the north to the east or nearly from the east to the west. On the west side of location is the Laojunta faulting zone lying on the north-east direction, and on the east side is the Sunjiagou faulting zone at the same direction, while on the north side is the front Maliangou faulting zone trending from east to west direction.All the faulting zones mentioned above are out of the caverns storage location, which influences the siting of caverns storage only. The structural map of caverns storage is shown in Fig. 1. Based on the project overview and the survey data, the weathered surface area of the mountain is classified as Class IV and Class V rock masses according to the degree of weathering, and the excavation area of the cavern is classified as Class III2 rock mass. The specific parameters are detailed in Table 1. According to hydrogeology investigation, the groundwater at the storage location mainly consists of the pore water in the loose rock mass and the fissure water in the bedrock including the one in the shallow rock strata distributing in the mesh type and the one in deep rock strata distributing in the nervation type.The project of underground petroleum storage includes items as follows: 2 construction roadways, 9 main caverns, 6 vertical shafts and 4 roadways for the water curtain, which is illustrated in Fig. 2. The entrances of 2 construction roadways locate at the southwest of storage with the designed altitude of 70m (in the standard altitude system of Huanghai in 1956, the same below). The roadways continue to the north part of storage and crosses at the north end on both east and west sides of storage, on which there are 3 branches along the main caverns going to the south of storage, with the total length of about 5819m and the end designed altitude of -30m. The average slope is about 13.3%. The cavern span of the construction is 8.3m, and the height of the cavern is 7.5m. The 9 main caverns are parallel at the direction of 45°west by north, every 3 ones of which are connected to be a whole tank by the construction roadways between them, with 3 tanks in total. The designed bottom altitude of main cavern is -50m, with the designed length of 560~650m and the designed cavern span of 20m and the designed height of 30m, being 5537m length in total 9 caverns. The net distance between the adjacent construction roadways of the main caverns is 25m, with the net distance between the adjacent main caverns of 30m. There are 2 vertical shafts in every one of 3group-cavern tanks, with the diameters of 3m and 5m respectively and the mouth of the well designed altitude of 100m and the depth of 110m. On the location of 30m above top of the main cavern (10m altitude), there are 4 horizontal roadways for water curtain, laying vertically to the main cavern direction, with the cavern span of 5m and the height of 4.5m and the total length of 2835m. The excavation work in underground engineering results in stress redistribution, deformation induced by unloading, damage and failure which leads to the unrecoverable EDZ. Based on the damage mechanics, Shao et al. (Shao et al., 1999) present an expression of damage tensor as follows:k k N k k k n n a a a m D     1 3 0 3 0 3 ) ( (1)Where, 0 a is the average radius of initial random cracks (penny crack), k a is the extended crack radius, N is the total groups number of cracks, k m is the density of kth group cracks. k n is the unit vector of kth group cracks.The deformation energy of EDZ is expressed as a linear equation of damage variable as follow:) ( 2 ) ( ) ( ) ( 2 ) ( ) , ( 2 D tr D tr tr tr tr D tr g D w                 (2)Where  ,  are the lame parameters of the undamaged material, g is a constant which represents the unrecoverable deformation,  ,  are the correction factors induced by damage in free energy formula.  ， D are the deformation and damage tensor respectively. The five parameters related to the constitutive model (  ,  , g,  ,  ) can all be determined through standard laboratory tests. The testing methods are detailed in reference (Homand-Etienne et al., 1998;Shao et al., 1999) regarding the determination of these parameters.The constitutive equation of rock in damage zone can be derived from equation 2 as follow:] [ 2 ] ) ( [ 2 ) , (              D D D tr I D tr I tr D g D w           (3)Where,  represents the stress tensor.The numerical simulation of the EDZ evolution and the associated hydro-mechanical coupling process was conducted using FLAC3D. The model incorporated the damage constitutive model described in Section 3.1. For the boundary conditions, the bottom of the model was fixed, the lateral boundaries were constrained in the normal direction, and the ground surface was free. Hydraulically, the lateral and bottom boundaries were set as no-flow boundaries, while the ground surface was set as a flux boundary to account for rainfall infiltration. Based on this model, the EDZ formation and evolution resulting from the construction of the water-sealed petroleum storage caverns were simulated. The model of caverns is shown in Fig. 3, the model boundary is a square with a length of As illustrated in Fig. 5, the evolution of the EDZ during excavation exhibits distinct spatiotemporal characteristics. When the excavation face advanced to a distance of 10m from the profile (Fig. 5a), initial damage zones emerged at the cavern arch and sidewalls, with a maximum depth of approximately 1.5 meters. This is primarily attributed to excavation unloading, which leads to a rapid reduction in radial stress and a concentration of tangential stress around the opening, thereby generating tensile and shear-tensile microcracks in the shallow surrounding rock. After the face advanced 50m past the profile (Fig. 5b), stress redistribution was largely complete, and the EDZ expanded significantly and stabilized. The damage zone depths at the arch and sidewalls increased to 3.2 meters and 2.8 meters, respectively. This evolutionary process quantitatively demonstrates that EDZ development primarily occurs within the influential zone ahead of the excavation face. Once this zone is passed, the damage range stabilizes, indicating the surrounding rock has reached a new equilibrium under a secondary stress state. This mechanism clarifies that timely support installation within this "disturbed zone" behind the face is most effective in curbing further EDZ development. 6a (systematic rock bolting), it is evident that rock bolts Arched ceiling: φ25@L=4~6m, 1.5m×1.5m, arranged in an interlaced pattern. Side wall: φ25@L=4~5m, 1.5m×1.5m, arranged in an interlaced pattern), through their suspending and combined arch effects, significantly restrained the relaxation of the shallow rock mass. This reduced the average EDZ thickness around the cavern from approximately 3.5 m to 2.5 m, a reduction of 28.6%. Furthermore, as shown in Fig. 6b, the application of anchor cables between adjacent caverns effectively controlled the deformation and damage interaction within the rock pillar, reducing the EDZ extent at the pillar center by an additional 40% compared to the unsupported case. This quantitative reduction in damage has significant hydro-mechanical implications: a smaller, less developed EDZ results in a more restricted zone of enhanced permeability, thereby significantly mitigating the expansion of the groundwater depression cone and establishing a sound rock mass foundation for the effectiveness of the water-sealing system during the operational period. Because of the anisotropy property of crack rock, the seepage process induced by the underground cavern excavation is very complex. In this paper, the saturated-unsaturated seepage process of underground water is studied under the influence of excavation disturbing of crack rock mass with considering rainfall infiltration.The seepage governing equation of underground water is shown as follows:Q z g P g k S t n          )] ( [ ) (      (4)Where, S is the saturability, k is the permeability, P is the underground water pressure, g is the gravitational acceleration, Q is the source sink term. In the equation, it is assumed that the water is incompressible, meaning that the density  does not change with time t. Based on the way mentioned above, according to the data of underground water level in the drill holes before construction, the initial condition and the boundary condition of seepage flow are inverted, and the seepage field in the construction location is simulated under the rainfall condition. The hydraulic head of initial seepage field is shown in Fig. The change and process of seepage under stress or deformation are by many with the results that it is the nonlinear relationship between stress and seepage Based on the deformation of Souley et al. (Souley et al., 2001) proposed the as a a k k a a a a k k    0 0 0 3 3 0 3 0 ) ( ) is the 0 k is the initial seepage is the of the crack length the seepage can to to the initial The of 0 k , and have been determined and through methods (Homand-Etienne et al., 1998;Shao et al., change of seepage induced by the crack in EDZ can be quantitatively using the equation addition, numerical simulation can be in on the change of seepage under the rock deformation, the relationship between the seepage and Fig. 9 the of the coupling process within the rock and the relationship between seepage and stress is shown as is the seepage of rock mass under the is the initial seepage of rock mass. The relationship between the seepage and stress is shown in Fig. The simulated seepage field in Fig. the of the groundwater flow induced by the cavern group excavation. Fig. and a large-scale groundwater depression cone on the storage caverns. The mechanism for this is that the underground as a sink for while the surrounding EDZ as a which greatly groundwater the caverns. The cone is its with the of the main caverns, and the groundwater is above the of the water curtain This the role of the water curtain system in water above the caverns. This seepage field a and quantitative for water-sealing that excavation has the Fig. a that quantitatively demonstrates the influence of considering rock mass damage on the seepage field In the considering damage (Fig. the hydraulic around the caverns is and the water is the permeability of the EDZ is (Fig. the hydraulic around the caverns and the groundwater by resulting in a more depression The in Fig. a the phreatic level in the model considering damage is significantly lower around the cavern arch and sidewalls than in the model the engineering investigation, the underground water in in the storage location were which are depth with the altitude range of controlled by the The ground water in and around the storage location are at the altitude range of with the lowest level of which is as the same as the one of on the south side of the the construction the water in all the including ones and ones are which a lot of data for to the model in this Because the data from drill holes in this the is as an The evolution of underground water under the disturbing influence of cavern excavation is studied in the evolution period of which is shown in Fig. The between the and simulated groundwater not the reliability but also the dynamic mechanism of the groundwater to excavation disturbance. The initial in water level on the to the excavation face into the zone of influence of The to with the process of EDZ and the associated in permeability as excavation This demonstrates that the groundwater is not but is controlled by the damage evolution process. The of the water level after the that the damage and permeability in the surrounding rock in this area and the seepage field reached a new dynamic