*Article* **Salt Cavern Thermal Damage Evolution Investigation Based on a Hybrid Continuum-Discrete Coupled Modeling**

**Kai Feng 1,2,3, Wenjing Li 3,\*, Xing Nan <sup>4</sup> and Guangzhi Yang <sup>5</sup>**


**Abstract:** The integrity and stability of salt caverns for natural gas storage are subjected to a gas cycling loading operation. The coupled effect of confining pressure and temperature on the response of the salt cavity surrounding the wall is essential to stability analysis. In this study, a hybrid continuum-discrete model accounting for the thermal-mechanical process is proposed to investigate the thermal-damage evolution mechanism towards a field case with blocks falling off the salt cavity. The salt cavity is modeled by continuum zones, and the potential damage zones are simulated by discrete particles. Three specimens at different locations around the surrounding wall are compared in the context of severe depressurization. The dynamic responses of rock salt, including temperature spatiotemporal variation, microscopic cracking patterns, and energy evolution exhibit spatial and confinement dependence. A series of numerical simulations were conducted to study the influence of microproperties and thermal properties. It is shown that the evolution of cracks is controlled by (1) the thermal-mechanical process (i.e., depressurization and retention at low pressure) and (2) the anomalous zone close to the brim of the salt cavity surrounding the wall. The zone far away from the marginal surrounding wall is less affected by temperature, and only the mechanical conditions control the development of cracks. This continuum/discontinuum approach provides an alternative method to investigate the progressive thermal damage and its microscopic mechanism.

**Keywords:** salt cavern; underground gas storage; continuum-discrete coupled method; thermalmechanical coupling; thermal damage

### **1. Introduction**

Salt rock is an attractive candidate for hosting energy storage, due to its favorable low permeability [1–4]. The salt caverns, which are constructed by the solution mining process, have been used as gas storage for several decades. However, when the salt cavern is subjected to gas-cycling loading, there are potential risks of fractures generation, block fall, and even collapse on the cavern roof [5–7]. The thermal damage evolution induced by the cycling loading process is still unknown and is worth investigating by innovative methods.

### *1.1. Problem Statement*

A field case with blocks falling off from the roof of a salt cavern in Jintan, Jiangsu province of China has been demonstrated by Li et al. [8]. Compared with the sonar results between the four years, there is a 4-m displacement at the shoulder of the salt cavity, indicating the cavern geometry has been altered (see Figure 1). The field engineers assumed

**Citation:** Feng, K.; Li, W.; Nan, X.; Yang, G. Salt Cavern Thermal Damage Evolution Investigation Based on a Hybrid Continuum-Discrete Coupled Modeling. *Sustainability* **2023**, *15*, 8718. https://doi.org/10.3390/ su15118718

Academic Editors: Mahdi Hasanipanah, Danial Jahed Armaghani and Jian Zhou

Received: 19 April 2023 Revised: 24 May 2023 Accepted: 24 May 2023 Published: 28 May 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

that it was the consequence of the thermal effects induced by the gas-cycling loadings. They suspected that the thermal stress led to spalling and resulted in the collapse of the roof eventually, after a 4-year operation. A proposed thermal-mechanical modeling in FLAC3D has been established, and the numerical work confirmed the thermal effect had impacted the stability of the salt cavern to some extent. Additionally, the operation conditions triggering the roof collapse are concluded. Despite the previous research achievements, some fundamental mechanisms of the thermal-dynamic response of rock salt are yet to be understood, particularly the onset and propagation of thermal micro-cracks, i.e., the thermal damage process induced by the thermal effect during gas-cycling loading in a salt cavern. that it was the consequence of the thermal effects induced by the gas-cycling loadings. They suspected that the thermal stress led to spalling and resulted in the collapse of the roof eventually, after a 4-year operation. A proposed thermal-mechanical modeling in FLAC3D has been established, and the numerical work confirmed the thermal effect had impacted the stability of the salt cavern to some extent. Additionally, the operation conditions triggering the roof collapse are concluded. Despite the previous research achievements, some fundamental mechanisms of the thermal-dynamic response of rock salt are yet to be understood, particularly the onset and propagation of thermal micro-cracks, i.e., the thermal damage process induced by the thermal effect during gas-cycling loading in a salt cavern.

between the four years, there is a 4-m displacement at the shoulder of the salt cavity, indicating the cavern geometry has been altered (see Figure 1). The field engineers assumed

*Sustainability* **2023**, *15*, 8718 2 of 27

**Figure 1.** A field case in Jintan, sonar results of Cavern L in the years 2009 and 2015, respectively [8]. (**a**) Sonar monitoring results of the cavity; (**b**) Displacement difference before and after cavity collapse. **Figure 1.** A field case in Jintan, sonar results of Cavern L in the years 2009 and 2015, respectively [8]. (**a**) Sonar monitoring results of the cavity; (**b**) Displacement difference before and after cavity collapse.

#### *1.2. Micromechanism of Thermal Damage*

*1.2. Micromechanism of Thermal Damage*  Thermo-mechanical responses of containment rocks are critical to the design and safe operation of underground energy storage [9]. The thermal effect induced by a gas injection-and-withdrawal process in a salt cavern was discussed comprehensively [10] (2019), who indicated that a tensile crack is possibly created at the surrounding rock of the salt cavity. Following Bérest's work, both experimental and numerical investigations were conducted to investigate the thermo-mechanical response of salt caverns during rapid cooling [11–17]. To ensure the integrity and stability of salt caverns, fractures, and rock damage should be avoided [18–21]. Rock damage is defined as the degradation of the macroscopic properties, such as strength, stiffness, etc [22–24]. The damage is the consequence of microcracks propagation, coalescence [25–29], and stiffness degradation [30]. The damage mechanics in rock engineering studies the evolution of damage that starts from microcracks and results in rupture failure in the macroscale of the structure [31]. Creep, one of the features of rock salt, is accompanied by microfractures [2]. Under cyclic loading, the fatigue-induced damage of salt rock at first is relatively small and then increases rapidly when it is close to failure [32]. Quick cyclic loading is prone to damage [33]. Ding et al. [34] investigated the grain-scale micromechanisms of the deformation of salt rock and concluded that viscoelastic and hysteretic behaviors are associated with the microprocesses at grain boundaries. Li et al. [35] investigated the damage pattern of rock Thermo-mechanical responses of containment rocks are critical to the design and safe operation of underground energy storage [9]. The thermal effect induced by a gas injectionand-withdrawal process in a salt cavern was discussed comprehensively [10] (2019), who indicated that a tensile crack is possibly created at the surrounding rock of the salt cavity. Following Bérest's work, both experimental and numerical investigations were conducted to investigate the thermo-mechanical response of salt caverns during rapid cooling [11–17]. To ensure the integrity and stability of salt caverns, fractures, and rock damage should be avoided [18–21]. Rock damage is defined as the degradation of the macroscopic properties, such as strength, stiffness, etc [22–24]. The damage is the consequence of microcracks propagation, coalescence [25–29], and stiffness degradation [30]. The damage mechanics in rock engineering studies the evolution of damage that starts from microcracks and results in rupture failure in the macroscale of the structure [31]. Creep, one of the features of rock salt, is accompanied by microfractures [2]. Under cyclic loading, the fatigue-induced damage of salt rock at first is relatively small and then increases rapidly when it is close to failure [32]. Quick cyclic loading is prone to damage [33]. Ding et al. [34] investigated the grain-scale micromechanisms of the deformation of salt rock and concluded that viscoelastic and hysteretic behaviors are associated with the microprocesses at grain boundaries. Li et al. [35] investigated the damage pattern of rock salt subjected to cycling loadings for CAES and concluded that different responses of the internal structure to fatigue and creep lead to the interaction between creep and fatigue. The micromechanisms of thermal damage in rock salt are even more complicated and challenging.

#### *1.3. Development of Hybrid Modeling*

Numerical simulation is considered an important method to study the stability performance during the operation's full life cycle and the associated mechanisms of salt caverns for energy storage [36]. The finite element method (FEM) has been utilized extensively for the salt-cavity integrity analysis. Many constitutive damage-mechanics models have been developed for salt rock. However, those constitutive models are not able to predict damage and other dynamic behaviors associated with microfractures.

The discrete element method (DEM) [37] was proposed later than FEM, as an alternative method, and has its own advantages over other methods. DEM treats rock materials as an assembly of rigid particles bonded with certain segregated contact modes. Discrete elements are independent and allow departing from the rock mass. When forces acting on the particles exceed their bond strength, the contact bond breaks. In DEM, the fracture is deciphered explicitly, the damage progress is reproduced as microcracks coalesce into macrofractures, and the dynamic process can be simulated simultaneously. When compared with FEM, the damage mechanism in DEM is not based on complex damage constitutive correlations; instead, the breakage between particles is simple, while the macroscopic damage is the consequence of individual "breakage", the assembling of discrete dynamic behavior, and the properties of the individual particle. Zhao et al. [38] proposed a grain texture model (GTM) with DEM, and, for the first time, this model can capture the major macromechanical characteristics of textured rock, including the failure process. Despite the unique features of DEM in dealing with the dynamic process, the computation efficiency is limited owing to the huge number of particles. Usually, computation efficiency depends on the number of particles and the size of the domain.

In order to overcome the shortcoming of computational efficiency in DEM, and meanwhile to investigate the dynamic behavior of rock engineering problems, hybrid models based on the continuum-discrete method are adopted. Hu et al. [39] (2021) employed a 3D continuum-discrete coupled method to establish the triaxial Hopkinson bar system, in which the steel bars and a cubic specimen were modeled by continuum zones and bondedparticle sections, respectively. This model was able to simulate the dynamic responses of the rock under different load conditions. Zhang et al. [40,41] (2017, 2019) demonstrated the capability of hybrid discrete-continuum modeling to simulate hydraulic fracturing propagation and interactions with natural fractures. The continuum-discrete coupled model has advantages in reproducing confinement loading via continuum and, meanwhile, accounts for the microstructure and the dynamic microbehavior of the target specimen via a discrete method. In general, the hybrid method is promising in addressing dynamic deformation, the fracturing behaviors of rock, and the related dynamic problems. Particularly, the continuum-discrete hybrid model established a correlation between macroscopic performance and the microscopic mechanism in the damage of rock.

The field case with blocks falling off from the roof of a salt cavern in Jintan, China was investigated with a thermal-mechanical model in FLAC3D (Li et al., 2021) [8]. However, the thermal effects on the micromechanism and microcracking evolution are lacking, and the progressive damage mechanism during the operation is still unknown. In this study, a 3D continuum-discrete coupled hybrid model is established to investigate the thermalmechanical dynamic behavior of the surrounding rock of the salt cavity subjected to gas cycling loadings. The salt cavern is represented by continuum zones, while the rock specimens on the roof with potential damage risks are simulated by discrete-element modeling. The hybrid model accounts for the influence of temperature variation. The thermal-mechanical coupling mechanism in the hybrid model is described in Section 2 of this paper. In Section 3, three specimens at different locations around the surrounding wall of the salt cavern are selected to understand the thermal damage process in the context of severe depressurization during operation. A parametric study is performed in Section 4 to discuss the influence of confining loading and the effects of rock properties parameters on cracking development. Based on this, the conclusions are presented in Section 5.

#### **2. Numerical Method 2. Numerical Method**

A 3D continuum-discrete coupled hybrid model is presented in this section. This model is an improved one after Li et al. (2021) [8]. The salt cavern is represented by continuum zones using fast Lagranian analysis of continua (FLAC3D), while the rock specimens on the roof where the collapse occurs are represented by discrete element modeling, implemented by particle flow code (PFC3D), shown in Figure 2. The coupling methodology is described as follows in detail. A 3D continuum-discrete coupled hybrid model is presented in this section. This model is an improved one after Li et al. (2021) [8]. The salt cavern is represented by continuum zones using fast Lagranian analysis of continua (FLAC3D), while the rock specimens on the roof where the collapse occurs are represented by discrete element modeling, implemented by particle flow code (PFC3D), shown in Figure 2. The coupling methodology is described as follows in detail.

the salt cavern are selected to understand the thermal damage process in the context of severe depressurization during operation. A parametric study is performed in Section 4 to discuss the influence of confining loading and the effects of rock properties parameters on cracking development. Based on this, the conclusions are presented in Section 5.

*Sustainability* **2023**, *15*, 8718 4 of 27

**Figure 2.** Hybrid continuum-discrete model for a salt cavern subjected to gas-cycling loading. (**a**) Salt cavern surrounding rock continuum part; (**b**) Selected continuous-discrete coupling part; (**c**)The size of the discrete part of the ball ranges from 0.04 to 0.06 m. **Figure 2.** Hybrid continuum-discrete model for a salt cavern subjected to gas-cycling loading. (**a**) Salt cavern surrounding rock continuum part; (**b**) Selected continuous-discrete coupling part; (**c**) The size of the discrete part of the ball ranges from 0.04 to 0.06 m.

#### *2.1. Model Description 2.1. Model Description*

tively.

To investigate the thermal damage induced by the operation condition, a 3D continuum-discrete coupled model is adopted using FLAC3D-PFC3D. As shown in Figure 3, the salt cavern surrounding rock is a continuous medium and FLAC3D is employed. To improve the computational efficiency, the simulation objective is one quarter of the entire salt cavity. The geometry of the FLAC3D domain is a length of 60 m, a width of 30 m, and a height of 110 m. The FLAC3D zone face and PFC3D wall are the interface, the geometry is 0.5m × 1m × 2m. The rock properties used in the FLAC3D continuum-based method are listed in Table 1. To investigate the thermal damage induced by the operation condition, a 3D continuumdiscrete coupled model is adopted using FLAC3D-PFC3D. As shown in Figure 3, the salt cavern surrounding rock is a continuous medium and FLAC3D is employed. To improve the computational efficiency, the simulation objective is one quarter of the entire salt cavity. The geometry of the FLAC3D domain is a length of 60 m, a width of 30 m, and a height of 110 m. The FLAC3D zone face and PFC3D wall are the interface, the geometry is 0.5 m × 1 m × 2 m. The rock properties used in the FLAC3D continuum-based method are listed in Table 1.

The PFC3D domain is embedded in a FLAC3D domain, consisting of 10,607 particles, and is located at the cavity shoulder. The location is selected due to its potential damage risk, where the block falls apart from the surrounding rock. For those selected sections, PFC is employed to mimic the dynamic and irreversible damage behavior of the rock salt, aiming at reproducing the evolution of microfractures. From the view of the discrete element method, the rock specimen is regarded as assemblies of discrete rigid particles connected with certain contacts. The movements of particles are governed by Newton's second law. The contact bond breaks when the contact force exceeds the tensile or shear strength of the contact bonds, caused by motion between adjacent particles. The input parameters of the FLAC3D and PFC3D models are listed in Table 1 and Table 2, respec-The PFC3D domain is embedded in a FLAC3D domain, consisting of 10,607 particles, and is located at the cavity shoulder. The location is selected due to its potential damage risk, where the block falls apart from the surrounding rock. For those selected sections, PFC is employed to mimic the dynamic and irreversible damage behavior of the rock salt, aiming at reproducing the evolution of microfractures. From the view of the discrete element method, the rock specimen is regarded as assemblies of discrete rigid particles connected with certain contacts. The movements of particles are governed by Newton's second law. The contact bond breaks when the contact force exceeds the tensile or shear strength of the contact bonds, caused by motion between adjacent particles. The input parameters of the FLAC3D and PFC3D models are listed in Tables 1 and 2, respectively.

**Figure 3.** The workflow for simulation of the hybrid continuum-discrete model.

**Table 1.** Rock mass properties for rock salt in FLAC3D.


cient ℃−1 5 × 10−<sup>5</sup> **Table 2.** Model parameters in PFC3D.

Linear thermal expansion coeffi-


#### Thermal expansion coefficient 1/℃ 1 × 10−5/0.7 × 10−5/0.3 × 10−5/1 × *2.2. Coupling Mechanism of FLAC-PFC*

10−<sup>6</sup> Figure 3 illustrates the workflow of the DEM-PFC coupled simulation for the thermalmechanical process of a salt cavern subjected to gas cycling loading. The continuum behavior of the salt cavity is simulated with FLAC, and the DEM model of the rock salt specimen enclosed by a surrounding wall is established by using the commercial code PFC3D. The thermal mode is coupled with the mechanical mode in each computation step. The thermal-mechanical interactive interface between FLAC and PFC is developed to account for the temperature-boundary settings in the FLAC zone, the wall–zone interface, and the particle wall in PFC. The thermal-mechanical coupling process occurs only in one direction (Figure 4): the changes in temperature induce the thermal strains which lead to the change of the mechanical stress, while the influence of mechanical changes on

*Sustainability* **2023**, *15*, 8718 7 of 27

the heat conduction calculation is not considered. The contact forces between particles, displacement, and the distribution of cracks are updated. The model is solved to a predetermined ratio.

**Figure 4.** Coupled continuum-discrete model for thermal-mechanical coupling process.

**Figure 4.** Coupled continuum-discrete model for thermal-mechanical coupling process. The DEM model accounts for the thermal expansion of particles with linear parallelbond contacts (Itasca, 2017; Li et al., 2016; Li et al., 2017) [26,42,43]. Figure 4 illustrates the heat conduction in a network composed of thermal reservoirs and pipes: yellow disks represent particles; red dots indicate the heat source; blue lines passing through the contact points of two circular particles are the active thermal pipes. For a bonded linear parallelbond contact, its mechanical contact is associated with thermal contact. Two consequences can be induced by the thermal-mechanical coupling: first, the particle size is modified due to thermal strain; second, the normal component of the contact force is affected by the temperature changes. The corresponding increment of particle radius ∆*R* induced by a temperature increment ∆*T* is:

$$
\Delta R = \alpha R \Delta T \tag{1}
$$

where *α* is the coefficient of linear thermal expansion, in the unit of 1/◦C. It is a micro property associated with the particle material.

The normal component of the force vector carried by the bond is assumed to be affected by the change in temperature. The relationship between the present parallel bond and active thermal pipe is expressed as:

$$
\Delta \overline{F}^{\prime} = -\overline{k}^{\prime} A \Delta U^{\prime} = -\overline{k}^{\prime} A \left( \overline{\pi} \overline{L} \Delta T \right) \tag{2}
$$

**Figure 5.** Coupled thermal-mechanical interaction at the FLAC-PFC interface. where *k n* is the bond's normal stiffness, *A* is the area of the bond's cross-section, *a* is the expansion coefficient of bond material, *L* is the bond length, and ∆*T* is the temperature increment, which equals the average temperature change of the two particles at two ends of the pipe associated with the contact bond.

Figure 5 demonstrates the interaction between the continuum FLAC3D zone and discrete particles. The interface (wall zone) consists of FLAC3D zone surfaces and PFC3D walls, which are created coinciding with the zone faces. The PFC walls are composed of edge-connected triangular faces, and the balls are in contact with the wall facet, wrapping the zone face. The coupling mechanism for FLAC-PFC works by updating the force system at facet vertices in PFC, which is determined by contact forces and moments at each ball– facet contact. The forces along with stiffness are communicated at grid points/nodes. The acting contact force and movement are distributed at grid points/nodes by equivalent forces. The grid points/nodes at the coupling wall zone and grids of zones move synchronously, and the updating force involves continuum zone computation in FLAC. Similarly, the deformation of the continuum zone leads to the movement of the coupling wall zone. In

response to the forces and velocities acting at the coupling wall zone, in DEM the particles displace and generate cracks when the stress at contact bonds exceeds the prescribed tensile or shear stress. **Figure 4.** Coupled continuum-discrete model for thermal-mechanical coupling process.

*Sustainability* **2023**, *15*, 8718 7 of 27

**Figure 5.** Coupled thermal-mechanical interaction at the FLAC-PFC interface. **Figure 5.** Coupled thermal-mechanical interaction at the FLAC-PFC interface. **3. Thermal Progressive Damage Evolution** 

#### **3. Thermal Progressive Damage Evolution** Figure 6 illustrates the cycling gas-loading process over the 5-year period. Five years

Figure 6 illustrates the cycling gas-loading process over the 5-year period. Five years is the time period for sonar monitoring for salt-cavity volume convergence. Particularly, at the time of 3.14 years, the temperature and pressure drop abruptly (from 16 MPa to 8 MPa) when there is gas withdrawal. However, not all the stages of the process of gas injection and withdrawal are xposure to the risk of thermal damage. Li et. al. (2021) [8] investigated cavern L and found that when the gas withdrawal is fast and followed by retaining low pressure, thermal cracking or even fractures occur. Therefore, we focus on the cracking development after 3.14-years of operation by comparing three different locations of surrounding wall of the salt cavity. is the time period for sonar monitoring for salt-cavity volume convergence. Particularly, at the time of 3.14 years, the temperature and pressure drop abruptly (from 16 MPa to 8 MPa) when there is gas withdrawal. However, not all the stages of the process of gas injection and withdrawal are 8xposure to the risk of thermal damage. Li et. al. (2021) [8] investigated cavern L and found that when the gas withdrawal is fast and followed by retaining low pressure, thermal cracking or even fractures occur. Therefore, we focus on the cracking development after 3.14-years of operation by comparing three different locations of surrounding wall of the salt cavity.

**Figure 6.** Gas temperature and pressure variations during five years of gas-cycling loading operation (Li et al., 2021) [8]. **Figure 6.** Gas temperature and pressure variations during five years of gas-cycling loading operation (Li et al., 2021) [8].

#### *3.1. Thermal Effect at Three Observed Locations*

Color Scale

*3.1. Thermal Effect at Three Observed Locations*  For the field case of cavern L (Li et al., 2021) [8], a sharp pressure drop occurs due to gas withdrawal at 3.14-years of operation, before and after, the pressure drop reaches For the field case of cavern L (Li et al., 2021) [8], a sharp pressure drop occurs due to gas withdrawal at 3.14-years of operation, before and after, the pressure drop reaches

7 MPa (Figure 7a). The three selected monitoring areas and corresponding locations in the hybrid model are (1) at the knee point of the cavity shoulder, i.e., the convexity; (2) the

Monitoring Area Location

7 MPa (Figure 7a). The three selected monitoring areas and corresponding locations in the hybrid model are (1) at the knee point of the cavity shoulder, i.e., the convexity; (2) the right part of the knee point; and (3) 11.8 m away from the cavity surrounding wall brim.

The initial temperature of rock-salt formation changes linearly with the depth at a temperature gradient of 2.55 ◦C/100 m. The salt cavern is located from −1030 m to −1080 m, assuming the ground temperature is 20 ◦C and the temperature of the cavity is approximately 45 ◦C. In response to the gas withdrawal and the consequent sharp pressure drop at 3.14 years of operation, the temperature decreases correspondingly. The temperature decreases gradually from the surface to the neighboring deeper domain in the formation. Overall, the temperature distribution of the monitoring areas after sharp gas depressurization exhibits obvious confinement dependence on the locations. Location 1 is the closest spot to the brim of the surrounding wall with the lowest temperature. Location 2 has the highest temperature; however, the temperature variation range is small. Locations 2 and 3 are less affected by the thermal effect. The closer to the salt cavity surrounding the wall brim, the more distinct the thermal influence it receives.

To assess the thermal effect on the response of the salt cavern surrounding wall, we proposed a thermal-mechanical factor *TMF*, which is defined as:

$$TMF = \frac{\text{Crack}(\text{Thermal} - \text{Mechanical}) - \text{Crack}(\text{Mechanical})}{\text{Total Crack}(\text{Thermal} - \text{Mechanical})} \tag{3}$$

*TMF* is the ratio of the difference between the crack number induced by the coupled thermal-mechanical effect and the crack number only induced by the mechanical effect to the total crack number. The higher *TMF* indicates that the thermal effect is dominant. *TMF* is close to 0.5, indicating both the thermal and mechanical working together. If *TMF* is approximately 0, the cracking is simply induced by the mechanical effect and is irrelevant to the thermal effect.

Figure 8 illustrates the crack development comparison between the thermal-mechanical coupled effect and the mechanical effect only for three different locations. The different particle colors represent the different detached broken particles formed as a result of cracking. At locations 1, 2, and 3, the *TMF* = 0.65, 0.18, and 0.58, respectively. Location 1, nearest to the cavity, is affected by the thermal effect most. Location 1 is sensitive to the thermal effect and the microcracking is controlled by the ball heat capacity, coefficient of thermal expansion, zone conductivity, conductivity coefficient of thermal contact mode, etc.

1

2

3

1

**3. Thermal Progressive Damage Evolution** 

tions of surrounding wall of the salt cavity.

*3.1. Thermal Effect at Three Observed Locations* 

tion (Li et al., 2021) [8].

Figure 6 illustrates the cycling gas-loading process over the 5-year period. Five years is the time period for sonar monitoring for salt-cavity volume convergence. Particularly, at the time of 3.14 years, the temperature and pressure drop abruptly (from 16 MPa to 8 MPa) when there is gas withdrawal. However, not all the stages of the process of gas injection and withdrawal are 8xposure to the risk of thermal damage. Li et. al. (2021) [8] investigated cavern L and found that when the gas withdrawal is fast and followed by retaining low pressure, thermal cracking or even fractures occur. Therefore, we focus on the cracking development after 3.14-years of operation by comparing three different loca-

**Figure 6.** Gas temperature and pressure variations during five years of gas-cycling loading opera-

For the field case of cavern L (Li et al., 2021) [8], a sharp pressure drop occurs due to gas withdrawal at 3.14-years of operation, before and after, the pressure drop reaches 7 MPa (Figure 7a). The three selected monitoring areas and corresponding locations in the hybrid model are (1) at the knee point of the cavity shoulder, i.e., the convexity; (2) the right part of the knee point; and (3) 11.8 m away from the cavity surrounding wall brim.

Location Temperature Distribution

**Figure 7.** *Cont*.

1

2

3

**Figure 7.** Three selected locations for observation of microcrack development in PFC (**a**); and the corresponding temperature (**b**) after sharp pressure drop due to gas withdrawal at 3.14 years of operation. **Figure 7.** Three selected locations for observation of microcrack development in PFC (**a**); and the corresponding temperature (**b**) after sharp pressure drop due to gas withdrawal at 3.14 years of operation.

3 are less affected by the thermal effect. The closer to the salt cavity surrounding the wall

To assess the thermal effect on the response of the salt cavern surrounding wall, we

*TMF* is the ratio of the difference between the crack number induced by the coupled thermal-mechanical effect and the crack number only induced by the mechanical effect to

(ℎ − ℎ) − (ℎ)

(ℎ − ℎ) (3)

The initial temperature of rock-salt formation changes linearly with the depth at a temperature gradient of 2.55℃/100 m. The salt cavern is located from −1030 m to −1080 m, assuming the ground temperature is 20℃ and the temperature of the cavity is approximately 45℃. In response to the gas withdrawal and the consequent sharp pressure drop at 3.14 years of operation, the temperature decreases correspondingly. The temperature decreases gradually from the surface to the neighboring deeper domain in the formation. Overall, the temperature distribution of the monitoring areas after sharp gas depressurization exhibits obvious confinement dependence on the locations. Location 1 is the closest

proposed a thermal-mechanical factor , which is defined as:

brim, the more distinct the thermal influence it receives.

=

vant to the thermal effect.

etc.

**Figure 8.** Comparison of crack development at three locations for coupled thermal-mechanical effect and mechanical effect only. **Figure 8.** Comparison of crack development at three locations for coupled thermal-mechanical effect and mechanical effect only.

the total crack number. The higher indicates that the thermal effect is dominant. is close to 0.5, indicating both the thermal and mechanical working together. If is approximately 0, the cracking is simply induced by the mechanical effect and is irrele-

Figure 8 illustrates the crack development comparison between the thermal-mechanical coupled effect and the mechanical effect only for three different locations. The different particle colors represent the different detached broken particles formed as a result of cracking. At locations 1, 2, and 3, the = 0.65, 0.18, and 0.58, respectively. Location 1, nearest to the cavity, is affected by the thermal effect most. Location 1 is sensitive to the thermal effect and the microcracking is controlled by the ball heat capacity, coefficient of thermal expansion, zone conductivity, conductivity coefficient of thermal contact mode,

### *3.2. Dynamic Response to Depressurization and Progressive Damage Characteristics*

Crack increments at three different locations are compared (Figure 9). The crack numbers at locations 1 and 3 are much higher than the number at location 2. The closer to the center of the cavern, the more affected by the thermal effect, and the more cracks formed. Particularly for location 1, close to the section of the most irregular geometrical of the caverns, which is the concentration of high stress, and has the potential local failure zones resulting from microcracking and might lead to damage. The tension crack is the dominant crack mode, as the tension cracking is the failure consequence of tension force (Equation (2)) induced by the thermal effect. Only the normal component of contact force between particles is affected by thermal expansion, resulting in tension failure. However, the influence of temperature is limited to some extent. The thermal effects are negligible beyond 10 m away from the brim of the surrounding wall of the salt cavity. At location 2, it is noticed that the increasing rate of tension crack decreases with time. Figure 10 illustrates the comparison of cracks development due to severe pressure drop, before and after. It can be seen clearly from the side view of the DEM specimen. After the depressurization, it is observed at location 1 that the cracks are generated and propagated starting from the lower left corner close to the cavity brim, the shear cracks are dominant. The sideview shows that the microcracking induced the discontinuities and was prone to particles falling off.

Location 3 is much closer to the brim of the cavity along the edge and cracks are formed and propagate on both sides of the DEM specimen. Location 2 is farther away from the cavity wall, hence there are fewer cracks, and the cracks propagate from the left lower corner dominated by tension cracks. sideview shows that the microcracking induced the discontinuities and was prone to particles falling off. Location 3 is much closer to the brim of the cavity along the edge and cracks are formed and propagate on both sides of the DEM specimen. Location 2 is farther away from the cavity wall, hence there are fewer cracks, and the cracks propagate from the left lower corner dominated by tension cracks.

*Sustainability* **2023**, *15*, 8718 12 of 27

*3.2. Dynamic Response to Depressurization and Progressive Damage Characteristics* 

Crack increments at three different locations are compared (Figure 9). The crack numbers at locations 1 and 3 are much higher than the number at location 2. The closer to the center of the cavern, the more affected by the thermal effect, and the more cracks formed. Particularly for location 1, close to the section of the most irregular geometrical of the caverns, which is the concentration of high stress, and has the potential local failure zones resulting from microcracking and might lead to damage. The tension crack is the dominant crack mode, as the tension cracking is the failure consequence of tension force (Equation (2)) induced by the thermal effect. Only the normal component of contact force between particles is affected by thermal expansion, resulting in tension failure. However, the influence of temperature is limited to some extent. The thermal effects are negligible beyond 10 m away from the brim of the surrounding wall of the salt cavity. At location 2, it is noticed that the increasing rate of tension crack decreases with time. Figure 10 illustrates the comparison of cracks development due to severe pressure drop, before and after. It can be seen clearly from the side view of the DEM specimen. After the depressurization, it is observed at location 1 that the cracks are generated and propagated starting from the lower left corner close to the cavity brim, the shear cracks are dominant. The

**Figure 9.** Temporal and spatial evolution of crack increments: A—Location 1; B—Location 3; C— Location 2; D—The ratio of tensile cracks to total cracks at Location 1; E—The ratio of tensile cracks to total cracks at Location 3; F—The ratio of tensile cracks to total cracks at Location 2. **Figure 9.** Temporal and spatial evolution of crack increments: A—Location 1; B—Location 3; C—Location 2; D—The ratio of tensile cracks to total cracks at Location 1; E—The ratio of tensile cracks to total cracks at Location 3; F—The ratio of tensile cracks to total cracks at Location 2.

In the PFC3D model, the fracture mass density is defined as the total fracture surface area per unit volume. The definition of fracture area is as followings:

If the domain is cubic, *L* is the length of the side of the cube. The number of fractures with sizes between *l*<sup>1</sup> and *l*<sup>2</sup> is given by:

$$n(l\_1 \le l \le l\_2) = \int \limits\_{l\_1}^{l\_2} n(l) \cdot L^3 dl = a \left( \frac{l\_2^{1-a} - l\_1^{1-a}}{1-a} \cdot L^3 \right) \tag{4}$$

The term *a* fixes the total fracture density by a range of fracture sizes.

$$d\_m(l\_c) \cong \int\_{l\_c}^{\infty} n(l) \cdot l^2 \cdot L^{-3} \cdot dl \tag{5}$$

The damage process is associated with cracking induced by cycles of gas pressurization and depressurization in the salt cavern.

To better evaluate the distribution of microcracks in discrete modeling, a normalized relative index *c<sup>d</sup>* is employed. The base point is 1.25 <sup>×</sup> <sup>10</sup>−<sup>3</sup> <sup>m</sup>2/m<sup>3</sup> , as it is the value of the fracture surface area per unit volume. The absolute index *c<sup>d</sup>* is the fracture mass density, which is defined as the fracture surface area per unit volume. The normalized concentration degree of crack index *c<sup>d</sup>* ranging from 0 to 1 is defined as the absolute index normalized by the maximum to minimum value among all the cracks. Figure 11 shows the comparison of three different locations of the surrounding wall of the salt cavity, which are: (a) location 1—at the knee point of the cavity shoulder, i.e., the convexity closest to the rim; (b) location 2—right part off the knee point, deep in the rock salt formation; (c) location 3—11.8 m away from the cavity surrounding the wall edge, the most farthest from the centerline of the cavity. The microcrack forms and distributes into the entire specimen in location 1. The crack distribution is more concentrated compared with the other two locations far away from the convexity. The microcracks become less intensive in the region that is not close to the inner cavity with no convexity–concavity. The damage

failure exhibits operation confinement and shape dependence. Figure 12 indicates the trend of the crack number using the normalized relative index *c<sup>d</sup>* . It shows that most of the cracks are with *c<sup>d</sup>* = 0.2. To investigate the heterogeneities of thermal cracking in rock salt, Figure 13 illustrates the rose diagrams of tensile and shear cracks at location 1, location 2, and location 3, respectively. The radial length of each bin indicates the number of shear or tensile cracks oriented within the corresponding angles. It shows that tensile cracks tend to initiate in the horizontal orientation at location 1 which is close to the cavity surrounding the wall. While the orientation of shear cracks in all three selected locations is uniformly distributed. *Sustainability* **2023**, *15*, 8718 13 of 27

**Figure 10.** Comparison of crack development at three locations before and after the sharp pressure drop. **Figure 10.** Comparison of crack development at three locations before and after the sharp pressure drop.

మ

The term *a* fixes the total fracture density by a range of fracture sizes.

area per unit volume. The definition of fracture area is as followings:

(ଵ ≤≤ଶ) = ∫భ

In the PFC3D model, the fracture mass density is defined as the total fracture surface

If the domain is cubic, *L* is the length of the side of the cube. The number of fractures

2 3 ( ) () *c*

relative index ௗ is employed. The base point is 1.25 × 10−3 m2/m3, as it is the value of the fracture surface area per unit volume. The absolute index ௗ is the fracture mass density, which is defined as the fracture surface area per unit volume. The normalized concentration degree of crack index ௗ ranging from 0 to 1 is defined as the absolute index normal-

The damage process is associated with cracking induced by cycles of gas pressuriza-

 () ⋅ ଷ = ቆଶ

ଵି − ଵ ଵି

*m c <sup>l</sup> d l n l l L dl* <sup>∞</sup> <sup>−</sup> ≅ ⋅⋅ ⋅ (5)

1− ⋅ ଷቇ (4)

with sizes between *l*1 and *l*2 is given by:

1

2

formly distributed.

**Figure 11.** Concentration degree of distribution of microcracks at (**a**) Location 1; (**b**) Location 2; (**c**) Location 3. **Figure 11.** Concentration degree of distribution of microcracks at (**a**) Location 1; (**b**) Location 2; (**c**) Location 3. *Sustainability* **2023**, *15*, 8718 15 of 27

ized by the maximum to minimum value among all the cracks. Figure 11 shows the comparison of three different locations of the surrounding wall of the salt cavity, which are: (a) location 1—at the knee point of the cavity shoulder, i.e., the convexity closest to the rim; (b) location 2—right part off the knee point, deep in the rock salt formation; (c) location 3—11.8 m away from the cavity surrounding the wall edge, the most farthest from the centerline of the cavity. The microcrack forms and distributes into the entire specimen in location 1. The crack distribution is more concentrated compared with the other two locations far away from the convexity. The microcracks become less intensive in the region that is not close to the inner cavity with no convexity–concavity. The damage failure exhibits operation confinement and shape dependence. Figure 12 indicates the trend of the crack number using the normalized relative index ௗ. It shows that most of the cracks are with ௗ = 0.2. To investigate the heterogeneities of thermal cracking in rock salt, Figure 13 illustrates the rose diagrams of tensile and shear cracks at location 1, location 2, and location 3, respectively. The radial length of each bin indicates the number of shear or tensile cracks oriented within the corresponding angles. It shows that tensile cracks tend to initiate in the horizontal orientation at location 1 which is close to the cavity surrounding the wall. While the orientation of shear cracks in all three selected locations is uni-

**Figure 12.** Crack number variation with concentration degree normalized relative index ௗ. **Figure 12.** Crack number variation with concentration degree normalized relative index *c<sup>d</sup>* .

Shear cracks Tensile cracks

**Figure 12.** Crack number variation with concentration degree normalized relative index ௗ.

 Location 1 Location 3 Location 2

0.0 0.5 1.0

Cd

0

10000

Number

20000

**Figure 13.** Orientation distribution of shear/tension cracks at different locations. **Figure 13.** Orientation distribution of shear/tension cracks at different locations.

#### *3.3. Energy Tracking 3.3. Energy Tracking*

2019).

l

<sup>μ</sup> Δ δ

When rock is subjected to loading, the energy dissipates with the elastic process. Numerical simulation is used to track the development of the strain energy, the slip energy, and the kinetic energy. Those energy partitions reflect the energy evolution during fractures propagation subjected to the thermal-mechanical coupling process. When the thermal effects are taken into account, the total stress is the summation of the mechanical stress and the thermal stress. When rock is subjected to loading, the energy dissipates with the elastic process. Numerical simulation is used to track the development of the strain energy, the slip energy, and the kinetic energy. Those energy partitions reflect the energy evolution during fractures propagation subjected to the thermal-mechanical coupling process. When the thermal effects are taken into account, the total stress is the summation of the mechanical stress and the thermal stress.

2

*K*

*E*

μ

( )<sup>l</sup> Fs *<sup>o</sup>* : linear shear force at the beginning of the timestep;

 μ

Strain energy, *EK*, is defined as the energy stored in the linear springs (Itasca Inc.,

( )<sup>2</sup> <sup>2</sup> <sup>l</sup> <sup>1</sup> <sup>F</sup> <sup>s</sup>

*K K*

, is defined as the total energy dissipated by frictional slip.

<sup>1</sup> : FF <sup>2</sup> *<sup>o</sup> E E*

: shear displacement component decomposed into slip component;

kinetic

= − + ⋅Δ

( ) ( )l lμ s s

2

1 2 δ

*E mv* = (8)

*n S*

(6)

(7)

 = + 

*l n*

*F*

*<sup>l</sup> Fn* : linear normal force;

Fs : linear shear force;

*Kn*: Normal stiffness;

<sup>l</sup> Fs : linear shear force;

*m*: particle mass; *V*: particle velocity.

*E*kinetic : kinetic energy of the particle.

μ

*KS*: Shear stiffness;

Slip energy *E*

Strain energy, *EK*, is defined as the energy stored in the linear springs (Itasca Inc., 2019).

$$E\_K = \frac{1}{2} \left( \frac{\left( {}^{F\_n}\right)^2}{{}^{K\_n}} + \frac{\left\| {}^{F\_s^1}\right\|^2}{{}^{K\_S}} \right) \tag{6}$$

*Fn l* : linear normal force;

F l <sup>s</sup>: linear shear force;

*Kn*: Normal stiffness;

*KS*: Shear stiffness;

Slip energy *Eµ*, is defined as the total energy dissipated by frictional slip.

$$E\_{\mu} := E\_{\mu} - \frac{1}{2} \left( \left( \mathbf{F\_s}^1 \right)\_o + \mathbf{F\_s}^1 \right) \cdot \Delta \delta^{\mu} \tag{7}$$

 Fs l *o* : linear shear force at the beginning of the timestep;

F l <sup>s</sup>: linear shear force;

∆*δ* <sup>µ</sup>: shear displacement component decomposed into slip component; *E*kinetic: kinetic energy of the particle.

$$E\_{\text{kinetic}} = \frac{1}{2}mv^2\tag{8}$$

*m*: particle mass;

*V*: particle velocity.

The evolution of these energy partitions is significantly impacted by the sharp pressure drop, as shown in Figure 14. At location 1, the strain energy decreases slightly followed by a dramatic rise and an increase of crack numbers. The initiation of cracks is positively correlated to the dissipation of strain energy, as the stain energy is stored in the contacts of neighboring particles. The higher rate of energy dissipation means that more cracks initiate. Slip energy rises gradually, indicating the increase of the relative deformation between particles. There is fluctuation in kinetic energy at all the three locations. The particle motion leads to the accumulation of energy to create new cracks. At first, the kinetic energy along the slip energy starts to rise rapidly, probably due to the thermal effect. Meanwhile, the strain energy, *Estrain*, assumed stored majority input work gradually increases. The breakage of contact bonds indicates more microcracks are initiated. With the reduction of kinetic energy, the kinetic energy is converted to other forms of energy, which leads to the consistent growth of the strain energy. Then, fewer new micro-cracks are developed, and the relative particle motion reduces. As a result, kinetic energy is decreasing gradually. The energy evolution can be attributed to the confinement condition (i.e., the pressure drop in this case), and is clearly associated with fracture generation and propagation.

To evaluate the thermal effect on the energy of the microscopic particle system, we proposed a thermal-mechanical energy factor. *TMK*, *TML*, *TMS*, and *TMK* are defined as:

$$TMK = \frac{A\_K - B\_K}{A\_K} \tag{9}$$

where *A<sup>K</sup>* is the kinetic energy of particles under thermal-mechanical force and *B<sup>K</sup>* is the kinetic energy of particles under mechanical force only.

and propagation.

The evolution of these energy partitions is significantly impacted by the sharp pressure drop, as shown in Figure 14. At location 1, the strain energy decreases slightly followed by a dramatic rise and an increase of crack numbers. The initiation of cracks is positively correlated to the dissipation of strain energy, as the stain energy is stored in the contacts of neighboring particles. The higher rate of energy dissipation means that more cracks initiate. Slip energy rises gradually, indicating the increase of the relative deformation between particles. There is fluctuation in kinetic energy at all the three locations. The particle motion leads to the accumulation of energy to create new cracks. At first, the kinetic energy along the slip energy starts to rise rapidly, probably due to the thermal effect. Meanwhile, the strain energy, ௦௧, assumed stored majority input work gradually increases. The breakage of contact bonds indicates more microcracks are initiated. With the reduction of kinetic energy, the kinetic energy is converted to other forms of energy, which leads to the consistent growth of the strain energy. Then, fewer new microcracks are developed, and the relative particle motion reduces. As a result, kinetic energy is decreasing gradually. The energy evolution can be attributed to the confinement condition (i.e., the pressure drop in this case), and is clearly associated with fracture generation

**Figure 14.** Energy evolution at three different locations. **Figure 14.** Energy evolution at three different locations.

kinetic energy of particles under mechanical force only.

sliding energy of particles under mechanical force only.

strain energy of particles under mechanical force only.

influence of the thermal effect gradually attenuates.

netic energy of particles is affected by dynamic mechanical force only.

To evaluate the thermal effect on the energy of the microscopic particle system, we proposed a thermal-mechanical energy factor. , , , and are defined as:

*<sup>A</sup> <sup>B</sup> TMK*

where *AK* is the kinetic energy of particles under thermal-mechanical force and *BK* is the

*<sup>A</sup> <sup>B</sup> TML*

where *AL* is the sliding energy of particles under thermal-mechanical force and *BL* is the

*<sup>A</sup> <sup>B</sup> TMS*

where *AS* is the strain energy of particles under thermal-mechanical force and *BS* is the

From Figure 15, it is not difficult to see that the *TMLs* of position 1 and position 2 are near the red line, indicating that the thermal effect has an influence on the sliding energy. According to the slope of each line, it can be seen that from location 1 to location 2 and then to location 3, which is, the farther away from the sensitive brim of the salt cavity, the

The larger the *TMK* is, the greater the influence of the thermal effect on the kinetic energy of particles, and the closer *TMK* is to one, indicating that the kinetic energy of particles is dominated by the thermal effect. *TMK* is close to zero, indicating that the ki-

*K K K*

*L L L*

*S S S*

<sup>−</sup> <sup>=</sup> (9)

*<sup>A</sup>* <sup>=</sup> <sup>−</sup> (10)

*<sup>A</sup>* <sup>=</sup> <sup>−</sup> (11)

*A*

The larger the *TMK* is, the greater the influence of the thermal effect on the kinetic energy of particles, and the closer *TMK* is to one, indicating that the kinetic energy of particles is dominated by the thermal effect. *TMK* is close to zero, indicating that the kinetic energy of particles is affected by dynamic mechanical force only.

$$TML = \frac{A\_L - B\_L}{A\_L} \tag{10}$$

where *A<sup>L</sup>* is the sliding energy of particles under thermal-mechanical force and *B<sup>L</sup>* is the sliding energy of particles under mechanical force only.

$$TMS = \frac{A\_S - B\_S}{A\_S} \tag{11}$$

where *A<sup>S</sup>* is the strain energy of particles under thermal-mechanical force and *B<sup>S</sup>* is the strain energy of particles under mechanical force only.

From Figure 15, it is not difficult to see that the *TMLs* of position 1 and position 2 are near the red line, indicating that the thermal effect has an influence on the sliding energy. According to the slope of each line, it can be seen that from location 1 to location 2 and then to location 3, which is, the farther away from the sensitive brim of the salt cavity, the influence of the thermal effect gradually attenuates. *Sustainability* **2023**, *15*, 8718 19 of 27

**Figure 15.** Thermal effect on energy of microscopic particle system at three different locations. **Figure 15.** Thermal effect on energy of microscopic particle system at three different locations.

#### **4. Results and Discussion 4. Results and Discussion**

1.0

1.1

1.2

Crack Number (×105

)

1.3

1.4

1.5

*4.1. Influence of Confining Pressure 4.1. Influence of Confining Pressure*

The thermodynamic response of the surrounding wall in a salt cavern subjected to a gas-cycling loading process is further investigated in the context of different confinement conditions. The operation pressure range of cavern L (Li et al., 2021) [8] is 16 MPa to 8 MPa. Three different confining pressures are schemed, which are 16 MPa, 12 MPa, and 8 The thermodynamic response of the surrounding wall in a salt cavern subjected to a gas-cycling loading process is further investigated in the context of different confinement conditions. The operation pressure range of cavern L (Li et al., 2021) [8] is 16 MPa to 8 MPa. Three different confining pressures are schemed, which are 16 MPa, 12 MPa, and 8 MPa, respectively.

MPa, respectively. Figure 16 shows the influence of the confinement pressure at location 1 (see Figure 7), which is at the knee point of the cavity shoulder (i.e., the closest location to the surrounding wall brim). The reduction of confining pressure (from 16 MPa to 8 MPa) augments the cracks generation during the depressurization process. The propagation and accumulation of microcracks depend on the variation of the confining pressure. Lower Figure 16 shows the influence of the confinement pressure at location 1 (see Figure 7), which is at the knee point of the cavity shoulder (i.e., the closest location to the surrounding wall brim). The reduction of confining pressure (from 16 MPa to 8 MPa) augments the cracks generation during the depressurization process. The propagation and accumulation of microcracks depend on the variation of the confining pressure. Lower confinement enhances the thermal damage of the salt cavity.

(**a**)

2.0 2.2 2.4 2.6 2.8 3.0

Time Step (×10<sup>3</sup>

)

 confinement condition 16 MPa confinement condition 12 MPa confinement condition 8 MPa

confinement enhances the thermal damage of the salt cavity.

confinement enhances the thermal damage of the salt cavity.

**Figure 15.** Thermal effect on energy of microscopic particle system at three different locations.

The thermodynamic response of the surrounding wall in a salt cavern subjected to a gas-cycling loading process is further investigated in the context of different confinement conditions. The operation pressure range of cavern L (Li et al., 2021) [8] is 16 MPa to 8 MPa. Three different confining pressures are schemed, which are 16 MPa, 12 MPa, and 8

Figure 16 shows the influence of the confinement pressure at location 1 (see Figure 7), which is at the knee point of the cavity shoulder (i.e., the closest location to the surrounding wall brim). The reduction of confining pressure (from 16 MPa to 8 MPa) augments the cracks generation during the depressurization process. The propagation and accumulation of microcracks depend on the variation of the confining pressure. Lower

**4. Results and Discussion** 

MPa, respectively.

*4.1. Influence of Confining Pressure* 

**Figure 16.** Influence of confining pressure on (**a**) crack increments; (**b**) failure mode; Yellow represents the intact particles; Blue represents the detached particles induced by tension force; Red represents the detached particles induced by shear force. **Figure 16.** Influence of confining pressure on (**a**) crack increments; (**b**) failure mode; Yellow represents the intact particles; Blue represents the detached particles induced by tension force; Red represents the detached particles induced by shear force.

The pressure is first increased to 8 MPa and then is retained at 8 MPa for a period. As shown in Figure 17, more microcracks initiate when the pressure is retained at 8 MPa. Tensile cracking is the dominant crack mode during the lower-pressure operation. However, the proportion of tensile cracks reduces with time due to the limitation of temperature extension. The more damage accumulates at that lower pressure and increases with the retaining time. Low pressure coupled with retaining time enhances the development of microcracks, ultimately resulting in a damage zone. The pressure is first increased to 8 MPa and then is retained at 8 MPa for a period. As shown in Figure 17, more microcracks initiate when the pressure is retained at 8 MPa. Tensile cracking is the dominant crack mode during the lower-pressure operation. However, the proportion of tensile cracks reduces with time due to the limitation of temperature extension. The more damage accumulates at that lower pressure and increases with the retaining time. Low pressure coupled with retaining time enhances the development of microcracks, ultimately resulting in a damage zone.

Crack Number

 Tension Number/Crack Number Shear Number/Crack Number

42.0 44.0 46.0 48.0 50.0 52.0 54.0 56.0 58.0

%

1.3

1.4

Crack Number (×105

)

1.5

(**a**)

Before retain Retain 1-month

16 MPa 12 MPa 8MPa (**b**)

resents the detached particles induced by shear force.

of microcracks, ultimately resulting in a damage zone.

**Figure 16.** Influence of confining pressure on (**a**) crack increments; (**b**) failure mode; Yellow represents the intact particles; Blue represents the detached particles induced by tension force; Red rep-

The pressure is first increased to 8 MPa and then is retained at 8 MPa for a period. As shown in Figure 17, more microcracks initiate when the pressure is retained at 8 MPa. Tensile cracking is the dominant crack mode during the lower-pressure operation. However, the proportion of tensile cracks reduces with time due to the limitation of temperature extension. The more damage accumulates at that lower pressure and increases with the retaining time. Low pressure coupled with retaining time enhances the development

**Figure 17.** Temporal and spatial evolution of microcracking at operation pressure of 8 MPa: (**a**) crack increment; (**b**) damage propagation process with retain time. **Figure 17.** Temporal and spatial evolution of microcracking at operation pressure of 8 MPa: (**a**) crack increment; (**b**) damage propagation process with retain time.

#### *4.2. Effect of Particle Microproperties 4.2. Effect of Particle Microproperties*

The development of cracks is affected by the microproperties of the particle, including the normal-to-shear stiffness ratio and the coefficient of interparticle friction. As shown in Figure 18, with the increase of the normal-to-shear stiffness ratio, more cracks initiate, leading to the lower yield strength. The development of cracks is affected by the microproperties of the particle, including the normal-to-shear stiffness ratio and the coefficient of interparticle friction. As shown in Figure 18, with the increase of the normal-to-shear stiffness ratio, more cracks initiate, leading to the lower yield strength.

(**a**)

0.0 0.5 1.0 1.5 2.0

 kratio 1.0 kratio 1.2 kratio 1.4 kratio 1.6

Strain (%)

0.0

1.0

2.0

Differential Stress (MPa)

3.0

generated.

**Figure 17.** Temporal and spatial evolution of microcracking at operation pressure of 8 MPa: (**a**) crack

The development of cracks is affected by the microproperties of the particle, includ-

shown in Figure 18, with the increase of the normal-to-shear stiffness ratio, more cracks

Retain 2-month Retain 3-month (**b**)

increment; (**b**) damage propagation process with retain time.

*4.2. Effect of Particle Microproperties* 

initiate, leading to the lower yield strength.

**Figure 18.** Effect of normal-to-shear ratio on (**a**) stress–strain curves; (**b**) microcrack increments. **Figure 18.** Effect of normal-to-shear ratio on (**a**) stress–strain curves; (**b**) microcrack increments.

Figure 19 illustrates the effect of the interfacial friction coefficient. In general, the friction coefficient has a slight effect on stress–strain curve. With the increase of the friction coefficient, the deformation resistance of particles is enhanced. It is shown that when the friction coefficient increases from 0.5 to 0.8, the larger friction between the particles makes it difficult to exceed the tensile or shear strength of the bonds. Therefore, fewer cracks are Figure 19 illustrates the effect of the interfacial friction coefficient. In general, the friction coefficient has a slight effect on stress–strain curve. With the increase of the friction coefficient, the deformation resistance of particles is enhanced. It is shown that when the friction coefficient increases from 0.5 to 0.8, the larger friction between the particles makes it difficult to exceed the tensile or shear strength of the bonds. Therefore, fewer cracks are generated.

> fric 0.3 fric 0.4 fric 0.5 fric 0.6

0.0

1.0

2.0

Differential Stress (MPa)

3.0

4.0

0.0 0.5 1.0 1.5 2.0

Strain (%)

(**a**)

generated.

ments.

wall of a salt cavern as well.

1.0

1.5

Crack Number (×105

)

1.0

1.5

Crack Number (×105

)

(**b**) **Figure 18.** Effect of normal-to-shear ratio on (**a**) stress–strain curves; (**b**) microcrack increments.

it difficult to exceed the tensile or shear strength of the bonds. Therefore, fewer cracks are

Figure 19 illustrates the effect of the interfacial friction coefficient. In general, the friction coefficient has a slight effect on stress–strain curve. With the increase of the friction coefficient, the deformation resistance of particles is enhanced. It is shown that when the

2.0 3.0

Time Step (×10<sup>3</sup>

)

 kratio 1.0 kratio 1.2 kratio 1.4 kratio 1.6

**Figure 19.** Effect of interfacial friction coefficient on (**a**) stress–strain curves; (**b**) microcrack incre-**Figure 19.** Effect of interfacial friction coefficient on (**a**) stress–strain curves; (**b**) microcrack increments.

### *4.3. Effect of Thermal Properties*

*4.3. Effect of Thermal Properties*  Simulations with various thermal conductivities, particle-specific heats, and thermal expansion coefficients are carried out to understand the effect of the thermal properties. All the other parameters are set to be the same as the case with pressure equal to 16 MPa at operation. As shown in Figure 20, the increase in thermal conductivity leads to an increment in crack number. On the other hand, a larger particle-specific heat enhances the energy stored between particles and slows down the heat conduction and the thermal expansion of the particles. Thus, fewer cracks are generated. The thermal-expansion coefficient is positively correlated with the number of cracks by augmenting the force carried by the bond. Hence, the thermal properties of the particles have an important impact on the crack generation and propagation, and the macro-behavior of the surrounding rock Simulations with various thermal conductivities, particle-specific heats, and thermal expansion coefficients are carried out to understand the effect of the thermal properties. All the other parameters are set to be the same as the case with pressure equal to 16 MPa at operation. As shown in Figure 20, the increase in thermal conductivity leads to an increment in crack number. On the other hand, a larger particle-specific heat enhances the energy stored between particles and slows down the heat conduction and the thermal expansion of the particles. Thus, fewer cracks are generated. The thermal-expansion coefficient is positively correlated with the number of cracks by augmenting the force carried by the bond. Hence, the thermal properties of the particles have an important impact on the crack generation and propagation, and the macro-behavior of the surrounding rock wall of a salt cavern as well.

(**a**)

2.0 3.0

Time Step (×10<sup>3</sup>

)

 Thermal Conductivity 6.5 W/(m·K) Thermal Conductivity 7.5 W/(m·K) Thermal Conductivity 8.5 W/(m·K) Thermal Conductivity 9.5 W/(m·K) ments.

*4.3. Effect of Thermal Properties* 

1.0

1.5

Crack Number (×105

)

wall of a salt cavern as well.

**5. Conclusions** 

spatiotemporal confinement dependence.

(**b**) **Figure 19.** Effect of interfacial friction coefficient on (**a**) stress–strain curves; (**b**) microcrack incre-

2.0 3.0

Time Step (×10<sup>3</sup>

)

 fric 0.5 fric 0.6 fric 0.7 fric 0.8

Simulations with various thermal conductivities, particle-specific heats, and thermal expansion coefficients are carried out to understand the effect of the thermal properties. All the other parameters are set to be the same as the case with pressure equal to 16 MPa at operation. As shown in Figure 20, the increase in thermal conductivity leads to an increment in crack number. On the other hand, a larger particle-specific heat enhances the energy stored between particles and slows down the heat conduction and the thermal expansion of the particles. Thus, fewer cracks are generated. The thermal-expansion coefficient is positively correlated with the number of cracks by augmenting the force carried

**Figure 20.** Microcrack increments with variations of (**a**) thermal conductivity; (**b**) particle-specific heat; and (**c**) thermal expansion coefficient. **Figure 20.** Microcrack increments with variations of (**a**) thermal conductivity; (**b**) particle-specific heat; and (**c**) thermal expansion coefficient.

via continuum modeling; meanwhile, the thermal-mechanical effects and processive damage were systematically investigated with the proposed micromechanical numerical

framework of the hybrid model. The main conclusions are as the followings:

A thermal-mechanical model was proposed to study a field case with blocks falling

Three specimens at different locations around the surrounding wall of a salt cavern were selected to understand the thermal-damage evolution process under the severe depressurization condition. The crack microscopic patterns were subjected to the cavern's anomalous geometry. The energy evolution associated with fracture creation exhibited

Cycling loading at a lower confining pressure and the longer retaining time of lowpressure led to progressive microcracking, which further resulted in the development of macrofractures and the formation of the damage zone. The tensile crack was the dominant crack mode. The anisotropic distribution of crack orientation was observed in the zone close to the edge of the surrounding wall, which was the consequence of tension failure

#### **5. Conclusions**

A thermal-mechanical model was proposed to study a field case with blocks falling off from the salt cavern roof using a 3D hybrid continuum-discrete method. This hybrid model has advantages in allowing for implementing the gas-loading operation conditions via continuum modeling; meanwhile, the thermal-mechanical effects and processive damage were systematically investigated with the proposed micromechanical numerical framework of the hybrid model. The main conclusions are as the followings:

Three specimens at different locations around the surrounding wall of a salt cavern were selected to understand the thermal-damage evolution process under the severe depressurization condition. The crack microscopic patterns were subjected to the cavern's anomalous geometry. The energy evolution associated with fracture creation exhibited spatiotemporal confinement dependence.

Cycling loading at a lower confining pressure and the longer retaining time of lowpressure led to progressive microcracking, which further resulted in the development of macrofractures and the formation of the damage zone. The tensile crack was the dominant crack mode. The anisotropic distribution of crack orientation was observed in the zone close to the edge of the surrounding wall, which was the consequence of tension failure induced by the thermal effect. The thermal effect induced by gas injection and withdrawal was limited. The zone far away from the marginal surrounding wall was less affected by temperature and only the mechanical conditions controlled the development of cracks.

The proposed 3D continuum-discrete coupled thermal-mechanical hybrid model overcomes the limitation of the continuum method and was capable of analysis of rock heterogeneity by accounting for the thermal-mechanical effect and naturally capturing the microcrack initiation and propagation by discrete modeling. The hybrid model was successfully applied in the assessment of thermal-damage evolution. This study can shed light on the analysis of thermal-mechanical coupling and the understanding of the influence of microcrack development on the macroscopic behavior of salt rock. In the future, the validity of the continuum-discrete coupled thermal-mechanical hybrid model should be verified by more field data and the thermal effect of the gas injection and withdrawal process needs to be investigated for energy storage in the deep stratum.

**Author Contributions:** Conceptualization, W.L.; Investigation, K.F. and W.L.; Methodology, X.N.; Software, K.F.; Supervision, G.Y.; Writing—original draft, W.L.; Writing—review and editing, K.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study is funded by the National Natural Science Foundation of China project (Grant NO. 42102304). Open Research Fund of State Key Laboratory of Shale Oil and Gas Enrichment Mechanisms and Effective Development—State Energy Center for Shale Oil Research and Development.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors also thank Jintan Underground Gas Storage, China Oil and Gas Piping Network Corporation for their advice on this study.

**Conflicts of Interest:** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
