Modeling and Parameter Optimization of Multi-Step Horizontal Salt Cavern Considering Heat Transfer for Energy Storage

Abstract

Highlights

What are the innovation of the article?

• A 3D multi-physics coupled flow model was established to consider the impact of heat transfer and turbulence flow.

• Heat transfer not only accelerates brine transport and diffusion, but also increases the salt rock dissolution rates with higher wall temperatures.

What are the main findings?

• As dissolution progresses, the previously formed cavity continues to expand at a gradually decreasing rate rather than remaining unchanged under the influence of brine flow.

• The final cavity shape is influenced by the injection flow rate, the concentration, the step distance, the number of steps, and the dissolution time of each stage. Each subsequent stage takes more time.

Abstract

Horizontal salt caverns represent a prime choice for energy storage within bedded salt formations. Constructing multi-step horizontal salt caverns involves intricate fluid and chemical dynamics, including salt boundary dissolution, cavern development, brine flow, heat transfer, and species transportation. In this paper, the influence of heat transfer and turbulent flow is considered in developing a 3D multi-physics coupled flow model for the construction of multi-step horizontal salt caverns. The feasibility and accuracy of the model are verified by comparisons with the field data of the Volgograd horizontal salt cavern. The effects of turbulent flow and heat transfer on the dissolution process are thoroughly analyzed. By analyzing the characteristics of the flow field, the brine concentration distribution, and cavern expansion, the results indicate a steady rise in cavity brine concentrations throughout the leaching phases, with the previously formed cavities continuing to enlarge during subsequent leaching stages, albeit at a diminishing rate of expansion. Furthermore, the results reveal that a larger injection flow rate results in a larger cavern volume, whereas higher injection concentrations result in smaller cavern volumes. While the step distance has a minimal impact on cavern volume, identifying the optimal step distance remains crucial. This analysis of construction parameters aims to provide valuable insights into the design and engineering practices involved in developing multi-step horizontal salt caverns for energy storage purposes.

1. Introduction

Despite the extremely low porosity and permeability of salt rock, underground salt caverns are ideal for various energy storage applications, such as natural gas, hydrogen, compressed air, strategic energy reserves, waste, and CO₂ [1,2,3]. Energy storage in salt rock offers significant advantages, including excellent sealing, high safety, and a low cost [4]. Salt caverns are extensively used for Underground Natural Gas Storage to regulate the regional and seasonal imbalances in gas demand [5,6] and Strategic Petroleum Reserves to ensure energy security [7,8,9]. Salt caverns can also be used for Compressed Hydrogen Energy Storage (CHES) and compressed air energy storage (CAES) [10,11,12]. These are two of the most cost-efficient large-scale energy storage technologies, which can buffer the electricity supply and demand cycles [13,14]. Currently, all the diabatic utility-scale CHESs and CAESs in the world have been built in salt caverns [15,16].

At present, there are only three salt caverns in China, and these are mainly used for storing natural gas [17,18]. As shown in Figure 1, the domestic production of natural gas has increased slowly, and the consumption of natural gas has risen rapidly in China since 2008 [19]. China’s natural gas imports and dependency are rising, which highlights the necessity for sufficient storage and peak shaving to ensure energy security. The construction of China’s three salt cavern gas storages employed the single-well-vertical method [20,21,22], a technique that requires a salt formation of substantial height to secure a large cavity volume. However, most of the salt rock formations in China deemed appropriate for constructing underground gas storage are continental deposits. These are characterized by thin single layers, many interlayers, and a high content of insoluble impurities [23]. The stability of salt caverns is affected by numerous unfavorable factors, resulting in more complex designs and construction processes for the salt caverns within layered salt rock formations [24,25]. Therefore, the bedded salt formations in China are unsuitable for large-scale straight salt caverns, making horizontal salt cavern construction a preferable alternative [26]. At present, the only successful case is the Volgograd horizontal salt cavern in Russia [27]. Compared with the two-well-horizontal (TWH) method [22,28], the multi-step horizontal construction method can be implemented in salt layers with a thickness of less than 80 m, can be utilize more salt formations, and does not require an antisolvent [10]. Figure 2 shows the construction process of a multi-step horizontal salt cavern. A directional well is drilled to connect with the vertical well. Then, a tube is put into the directional well. The process is to inject freshwater from the inlet and replace high-concentration brine from the outlet of the vertical well, and injection and discharge are cycled continuously. The injection tube is removed at a certain distance for a certain time, and then the next stage is leaching. After several steps, a multi-step horizontal cavern is formed finally. The modeling of multi-step horizontal salt caverns construction is less frequently conducted due to the complexity of their shape and flow fields. Presently, the theoretical studies on the construction of multi-step horizontal caverns primarily concentrate on numerical simulations and extensive indoor salt rock experimentation. This technology remains in its nascent stages, straddling between theoretical inquiry and preliminary field explorations, and has yet to reach maturity [29,30,31].

Many models and codes have been developed to simulate the construction of vertical caverns, including SANSMIC, WinUbro, and PCL [32]. However, there are few models for simulating horizontal salt caverns due to their complex shape and flow field. The multi-step horizontal construction method has been proposed [33], highlighting the importance of controlling the step distance to maintain the integrity of step caverns. Based on the floating plume model to describe the flow field, Saberian [34] initially simulated the expansion of a horizontal cavern’s middle segment using a horizontal cylindrical grid, disregarding the vertical concentration gradient of the brine. Russo [35] has developed the HORSMIC model to simulate the shape of the middle section of a horizontal cavern. Both the limits are that the flow area in the cavity is divided in advance and needs to be adjusted manually. Liu [28] conducted physical simulations of TWH salt cavern construction using a high-strength steel mold to shape large rock salt specimens, making the construction process visible and easily observable. The influences of the injection rate, the transferring of the injection well, the oil blanket, and the injection position on cavern development were discussed. Kang [36] performed a physical simulation experiment to investigate the characteristics of horizontal leaching, focusing on how the leaching parameters, like the leaching rate, the brine discharge position, and the back-step length, affect the cavern shape and the brine discharge concentration. This experiment is based on similar principles and uses small salt bricks for the cavity construction experiments. Wan [37] simulated the TWH salt cavern construction process using the TWHSMC code program based on a Navier–Stokes equation to model a U-shaped cavern from a U-shaped borehole. While this model accounts for the fluid flow influences, it overlooks the effects of turbulent flow and heat transfer. Its application to multi-step modeling is limited due to the fixed direction of boundary movements. Li [38] established an M-SHSCLS model using a time difference method and a finite volume method based on the balance principle, simulated the cavity shape of a salt cavern in multi-step horizontal construction, and optimized the cavity shape. However, the brine concentration is simplified, and the influence of flow field on brine distribution and the dissolution rate was ignored. The previously formed cavity did not expand at the next leaching stage. Zeng [32] established a TWHMPC model based on multiple coupled equations of dissolution, transport, and fluid flow to simulate a horizontal salt cavern. This model found the influence of turbulence to be not negligible at the initial stage of solution mining, without considering the influence of heat transfer. This model only simulated the first stage of the cavity leaching process and did not involve multi-step cavity construction.

At present, no one has considered the impact of temperature on dissolution and expansion for construction simulations of salt caverns, including vertical caverns, TWH caverns, and horizontal caverns. However, the construction of a multi-step horizontal salt cavern is a complicated process of fluid dynamics, chemical dynamics, and thermodynamics [39,40]. Heat transfer will affect the fluid flow and brine concentration distribution in the cavity, and also affect the expansion rate of the cavern boundary. The depth of a salt cavern built in a salt rock formation is about 1500~2500 m, and the cavern temperature is about 70~90 °C. The temperature of injection freshwater is about 20 °C. Therefore, there is a large temperature difference between the injected fluid and the brine in the cavity, so there is a thermal convection phenomenon between the forced convection fluid, as shown in Figure 3. The dissolution rate increases with the increase in temperature. The dissolution rate at 80 °C is 2–3 times larger than that at 20 °C [41,42,43]. In addition, there are few studies on multi-step horizontal cavern construction, and there is currently no criterion for the design of multi-step horizontal caverns.

In this study, two critical elements of heat transfer and turbulent flow are considered in developing a 3D multi-physics coupled flow model to simultaneously quantify salt rock dissolution, cavity expansion, heat transfer, brine convection, turbulence flow, and diffusion effects during the construction of multi-step horizontal salt caverns. Utilizing software COMSOL 5.6, the proposed model integrates the flow field, concentration field, temperature field, and deformation field, facilitating the simultaneous resolution of multiple physical fields. Adaptive meshing and automatic meshing methods are used to describe cavity boundary expansion. The feasibility and accuracy of the model are verified by comparisons with the field data of the Volgograd horizontal salt cavern. Analysis encompasses both the mechanism and the influence of heat transfer on the flow field, brine concentration, and cavity expansion, as well as the characteristics of the turbulent flow field, concentration field, and cavity expansion. Several groups of simulations have been executed for the optimization of the technological parameters, with a discussion on how these parameters influence the cavern volume and the outlet brine concentration.

2. Theory and Method

Salt cavern solution mining is a complicated process of fluid dynamics and chemical dynamics, including salt rock dissolution, cavern expansion, brine flow, heat transfer, and species transportation. The reaction processes occur simultaneously and interact with each other. Based on the theory of flow convection diffusion, fluid mechanics, chemical dynamics, and thermal dynamics, a multi-physical coupled flow model is established.

2.1. Assumptions

In the previous studies, there are some assumptions for simplifying the problem for simulation. And assumptions have been proposed as follows [32,37]:

• (1). The influence of anisotropy on the microstructure on the rock salt dissolution process is neglected, such as stratigraphic dipping, salt crystal direction, bedding, etc.

• (2). Insoluble is 100% precipitated, and pure salt rock is considered with no insoluble sediments.

• (3). There is no heat exchange, and the system remains isothermal.

In the simplified assumptions, each additional assumption encourages the deviation from reality. Considering the accuracy and efficiency of the calculations, it is essential to assume that salt rock is anisotropic. The coupled model takes account the effects of fluid flow, primarily investigating the impact of flow on leaching process. Including insoluble materials would make the simulation almost infeasible due to the significantly increased complexity. However, it is practical to incorporate the influence of temperature in addition to fluid flow, which can enhance the accuracy of the simulation.

In this study, a 3D multi-physics coupled flow model is established to describe the process of multi-step horizontal salt cavern construction. The assumptions were corrected as follows.

• (1). The salt rock is isotropic with no cracks and fractures, ignoring the influence of anisotropy during the dissolution process.

• (2). Insolubility is fully precipitated, and we ignored its influence on fluid flow and brine convection. The salt rock is pure with no insoluble sediments.

2.2. The Governing Equations

The multi-step leaching process is a multi-physical field coupled fluid flow process which includes complex fluid flow, mass convection and diffusion, heat transfer, and a boundary expansion process.

2.2.1. Dissolution Equation

The previous dissolution equation only considered the effect of concentration, but the effect of temperature is must be considered [41,43]. The dissolution rate of salt boundary can be expressed as Equation (1), which describes the process of cavern development.

(1)${\mathit{\Phi }=\mathrm{C}}_{dis}\left(c_{sat}-c\right)\cdot \ln {\displaystyle \frac{{\left(T-T_0\right)}^{0.44}}{e}}$

where $\mathit{\Phi }$ is the dissolution flux at the salt–brine interface, mol/(m²·s); ${\mathrm{C}}_{dis}$ is the coefficient of dissolving rate, m/s; $c_{sat}$ is the concentration of saturated brine, mol/m³; $c$ is the concentration of brine, mol/m³; T is the temperature of the salt cavern, K; and T₀ is the freezing temperature of water, K.

2.2.2. Heat Transfer Equation

The heat transfer equation describes the heat transfer processes in fluids and cavern walls. During the dissolution process, there is a temperature difference between the injected freshwater, the fluid in the cavern, and the salt solid wall. There will be heat convection exchange between the fluids of different temperatures and heat conduction between the fluids and the salt wall, expressed as Equation (2).

(2)$\rho {\mathrm{C}}_p{\displaystyle \frac{\partial T}{\partial t}}+\rho {\mathrm{C}}_p{u}_i{\displaystyle \frac{\partial T}{\partial x_i}}+{\displaystyle \frac{\partial }{\partial x_i}}\left(-k{\displaystyle \frac{\partial T}{\partial x_j}}\right)=Q$

where $\rho$ is the fluid density, kg/m³; $C_p$ is the fluid heat capacity at constant pressure, J/K; $Q$ is the heat source, W/m³; $u_i$ is the fluid velocity, m/s; and $k$ is fluid thermal conductivity, m/K.

2.2.3. Species Transport Equation

Equation (3) can express mass convection and diffusion [32].

(3)${\displaystyle \frac{\partial }{\partial t}}\left(\rho c\right)+{\displaystyle \frac{\partial }{\partial x_i}}\left(\rho c{u}_i\right)={\displaystyle \frac{\partial }{\partial x_i}}\left(D\rho {\displaystyle \frac{\partial c}{\partial x_i}}\right)+{\displaystyle \frac{\partial \left(\rho \mathrm{\Phi }\right)}{\partial x_i}}$

where D is the diffusion coefficient for NaCl, m²/s; $u_i$ is the brine velocity component in the i direction, m/s, i = 1, 2, 3.

2.2.4. Interface Movement Equation

The cavern expansion can be expressed as Equation (4). The interface moves toward the interior of the rock domain as dissolution goes on.

(4)$w_n=\mathit{\Phi }M/\rho$

where $w_n$ is the moving velocity of the interface in the normal direction, m/s; $\mathit{M}$ is the molar mass of sodium chloride, kg/mol.

2.2.5. Fluid Flow Equations

The injection rate is usually above 150 m³/h, and Reynold number ranges from 625 to 1.22 × 10⁶. In the previous studies, turbulence flows occur during the process of initial first stage construction [32]. The Navier–Stokes κ-ε equations describe the complex turbulence flow in the horizontal salt cavern.

The mass conservation equation [32,40]:

(5)${\displaystyle \frac{\partial \rho }{\partial t}}+{\displaystyle \frac{\partial }{\partial x_i}}\left(\rho u_i\right)=0$

The momentum equation:

(6)${\displaystyle \frac{\partial }{\partial t}}\left(\rho u_i\right)+{\displaystyle \frac{\partial }{\partial x_j}}\left(\rho u_i{u}_j\right)=-{\displaystyle \frac{\partial p}{\partial x_i}}+\mu {\displaystyle \frac{{\partial }^2{u}_i}{\partial x_j\partial x_j}}+{\displaystyle \frac{\partial }{\partial x_j}}\left(-\rho \overline{u_i^{\prime }{u}_j^{\prime }}\right)$

The turbulent kinetic energy equation:

(7)${\displaystyle \frac{\partial }{\partial t}}\left(\rho k\right)+{\displaystyle \frac{\partial }{\partial x_i}}\left(\rho k{u}_i\right)={\displaystyle \frac{\partial }{\partial x_j}}\left[\left(\mu +{\displaystyle \frac{{\mu }_t}{{\sigma }_k}}\right){\displaystyle \frac{\partial k}{\partial x_j}}\right]+G_k+G_b-\rho \epsilon -Y_M$

The turbulent dissipation rate equation:

(8)${\displaystyle \frac{\partial }{\partial t}}\left(\rho \epsilon \right)+{\displaystyle \frac{\partial }{\partial x_i}}\left(\rho \epsilon u_i\right)={\displaystyle \frac{\partial }{\partial x_j}}\left[\left(\mu +{\displaystyle \frac{{\mu }_t}{{\sigma }_{\epsilon }}}\right){\displaystyle \frac{\partial \epsilon }{\partial x_j}}\right]+C_{1\epsilon }{\displaystyle \frac{\epsilon }{k}}\left(G_k+C_{3\epsilon }{G}_b\right)-C_{2\epsilon }\rho {\displaystyle \frac{{\epsilon }^2}{k}}$

where $u_i$ is the brine velocity component in the i direction, m/s, i = 1, 2, 3; $u_j$ is the brine velocity component in the j direction, m/s, j = 1, 2, 3; $p$ is the pressure, Pa; ${\mu }_t$ is the turbulent viscosity, Pa·s; $u_i^{\prime }$ is the fluctuating velocity component in i direction, m/s, i = 1, 2, 3; $u_j^{\prime }$ is the fluctuating velocity component in the j direction, m/s, j = 1, 2, 3; $k$ is the turbulence kinetic energy, m²/s²; $\epsilon$ is the turbulence dissipation rate, m²/s³; $G_k$ is the turbulence kinetic energy due to the mean velocity gradients, kg/(m·s³); $G_b$ is the turbulence kinetic energy due to the buoyancy, kg/(m·s³); $Y_M$ is the contribution of the fluctuating dilatation to the dissipation rate, kg/(m·s³); $C_{1\epsilon }$ is the empirical coefficient measured by experiment and is dimensionless; $C_{2\epsilon }$ is the empirical coefficient measured by experimenters and is dimensionless; $C_{\mu }$ is the experimental coefficient and is dimensionless; ${\sigma }_k$ is the turbulent Prandtl number for kinetic energy and is dimensionless; ${\sigma }_{\epsilon }$ is the turbulent Prandtl number for the dissipation rate and is dimensionless; and $u_g$ is the velocity component in the gravity direction, m/s.

2.2.6. Auxiliary Equation

The relationship between the density and concentration of the brine can be expressed as Equation (9) [22]. The brine concentration changes with fluid flow and solute diffusion, and the density of brine also changes.

(9)$\rho ={\rho }_0+{\displaystyle \frac{c}{c_{sat}}}\left({\rho }_{sat}-{\rho }_0\right)$

2.3. Definite Conditions

The initial condition can be expressed as

(10)${\left.u\right|}_{t=0}=0,{\left.c\right|}_{t=0}=c_{sat},{\left.p\right|}_{t=0}=p_0,{\left.T\right|}_{t=0}=T_{wall}$

The boundary conditions can be expressed as

(11)${\left.u\right|}_{in}=u_0,{\left.c\right|}_{in}=0,{\left.p\right|}_{out}=p_0$

(12)${\left.T\right|}_{in}=T_{\mathrm{in}},{\left.T\right|}_{out}={\left.T\right|}_{wall}=T_{wall}$

where $p_0$ is the pressure of the salt formation, Pa; $T_{\mathrm{in}}$ is the inlet temperature, K; and $T_{wall}$ is the cavern wall temperature.

2.4. Geometry Model and Parameters

To simulate the leaching process of the multi-step horizontal salt cavern, an idealized geometry model is assumed, as shown in Figure 4a. The initial geometry is composed of a reflux cavity and solution cavity. The boundary is the normal direction deformation boundary. The values of the formation and geometry parameters are adopted [32,44,45], as shown in Table 1. The depth of the formation is 1500 m. Thus, the temperature and pressure of the cavern is the same as that of the formation. The concentration of the initial brine in the salt cavern is saturated. The leaching duration of every stage is 50 days. Whole leaching includes 5 stages.

2.5. Numerical Algorithms

The coupled flow model offers theoretical formulations for the 3D multi-step horizontal salt cavern leaching process. The adaptive meshing method addresses the moving boundary challenge within the coupled calculations, as shown in Figure 5. Initially, at time T₀, the geometric model undergoes adaptive meshing following the input of geological parameters, after which the mathematical model is initialized for simulation. A maximum mesh element is set, and the geometric model is meshed, followed by iterative calculations, with each step recording the mesh quality. At the subsequent moment T₀ + ∆T, the distribution of the brine concentration, temperature field, and flow field are determined, facilitating cavern boundary expansion via the brine concentration and temperature field. At the next time T₀ + ∆T, the brine concentration distribution, the temperature field, and the flow field can be obtained, and then cavern boundary expansion is achieved through the brine concentration distribution and the temperature field. The new geometric model is extracted after the cavern boundary expansion. At the time T₀ + ∆T′, the new geometric model is re-meshed adaptively based on the previously set maximum mesh element when the mesh quality is less than 0.2. If not, the calculation will continue.

The numerical algorithms of the coupled model are shown in Figure 6. The segregated algorithm is used to solve the equation sequentially. The flow field is solved first using the Navier–Stokes κ-ε equations to obtain the velocity field, pressure distribution, turbulent kinetic energy, and turbulent dissipation. Then, the temperature field is calculated using the heat transfer equation. The brine concentration field can be obtained using the species transport equation. Next, the displacement of the boundary is calculated using the interface movement equation, and the new cavern boundary is obtained by updating geometry and re-meshing at the same time. After iterative calculation, the cavity shape of one stage can be obtained.

Then, the brine concentration field is calculated using species transport equations based on the parameter of the flow field. Next, the displacement of the boundary is calculated based on brine concentration distribution, and the new cavern boundary is obtained by updating the mesh at the same time. After iterative calculation, the cavity shape of one stage can be obtained. After the calculation of one stage, the last cavern shape is exported, and the new initial and boundary conditions are defined and calculated iteratively until all the stages have ended. After multi-step cyclic iterative calculation, the cavern expansion, the brine concentration distribution, the temperature distribution, and the velocity distribution are finally obtained.

2.6. Comparison and Verification

The Volgograd underground gas storage represents the only successful implementation of horizontal multi-step salt cavern gas storage to have been reported. To verify the validity and accuracy of the coupled flow model, the process of constructing the Volgograd salt cavern is simulated with the same parameters of field [27]. By using sonar detection technology, a sonar detector is deployed after the cavity formation stage ends. This process captures the propagation and reflection time of sound waves in different media, subsequently obtaining data on the cavity volume and morphology. The outlet brine concentration and cavern volume data are monitored and compared. As shown in Figure 7 and Figure 8, the proposed model predictions are consistent with the Volgograd salt cavern monitored data with the same trend, and the deviations are within 9.4% for the outlet brine concentration and 6.3% for the cavern volume. During the process of multi-step cavity construction, the brine concentration initially experienced a rapid decline, which subsequently rose to approximately 240 g/L by day 90, and thereafter incremented slowly, stabilizing around 270 g/L.

For the cavern shape, due to the limitations of current sonar detection technology in mapping the shape of horizontal sections, the 3D shape of field horizontal caverns created through multi-step leaching is not documented in the existing literature. The side cavern shape of the proposed model is compared with a large-scale leaching experiment on a salt dome [46] and the M-SHLS model [38], as shown in Figure 9. The side of the cavern is crown-shaped above and cone-shaped below, as shown in both the experiment and model results.

Due to the limited data available on horizontal salt cavern construction, a degradation verification is performed by simulating the process of multi-step vertical salt cavern construction and comparing the data with software PCL 1.0 to further confirm the model’s reliability. As illustrated in Figure 10, the outlet brine concentration and cavity volume evolutions with the leaching time were analyzed and compared. The results demonstrate that the outlet brine concentration and the cavity volume data obtained from the coupled model are largely consistent with the data of the PCL 1.0. As dissolution progresses, the brine concentration initially increases rapidly, and then rises steadily at a slower rate, while the cavity volume initially increases slowly, and then rises steadily at a faster rate. The comparison results further verified the reliability of the coupled model.

Above all, through comparisons of the outlet brine concentration, the cavern volume, and the side shape, the proposed coupled flow model is reliable and has a high accuracy.

3. Results and Analysis

3.1. Analysis of Heat Transfer

During solution mining, thermal convection and conduction arise from the temperature difference between the injected freshwater, the fluid within the cavern, and the cavern walls, impacting fluid flow, the brine concentration, and the salt rock wall’s dissolution rate. Figure 11 illustrates the directions of thermal conduction and fluid flow. The injected fluid induces forced convection, creating a shearing effect between the fluids with temperature differences, thus aligning the forced thermal convection direction perpendicularly to the fluid flow direction in the cavern.

Forced convection impacts brine concentration distribution, initiating transportation and diffusion towards equilibrium. This process, accelerated by forced thermal convection, aligns with the brine diffusion direction, making it perpendicular to the brine concentration gradient, as shown in Figure 12. Such convection hastens brine transportation and dissolution, with higher wall temperatures increasing the salt rock dissolution rates, thus accelerating cavity construction. The cavity expansion caused by different formation temperatures were compared, as shown in Figure 13. A constant injection fluid temperature alongside higher formation temperatures and greater temperature differences leads to faster cavity construction rates and increased expansion heights.

3.2. Analysis of Flow Field and Flow Characteristics

Injecting freshwater into the cavity, which is akin to a submerged wall jet, as illustrated in Figure 14, creates complex and variable flow fields due to water injection, fluid convection, and boundary expansion. Initially, the cavern is categorized into the main flow, circulation, and reflux regions based on the flow characteristics, delineated by the null streamwise velocity line (v = 0) and the inflection point (∂2v/∂y2 = 0) of the velocity profile [32,40]. The main flow region exhibits a higher velocity due to momentum from the injected freshwater, which quickly diminishes post-injection. Meanwhile, brine flow forms a vortex due to forced convection; this initial flat vortex expands as the cavity grows upwards, though its intensity decreases, and its center shifts towards the outlet. Figure 15 displays the flow velocity, streamline distribution, and vortex shape during the second leaching stage. Forced convection at the injection point creates a single vortex from the first stage. As the flow momentum diminishes, another vortex begins to form in the dissolution section, eventually merging into a larger single vortex as the cavity expands.

During the later leaching stages, despite cavern expansion, the dimensions of flow and reflux regions stabilize, with a constant flow rate of 200 m³/h. The flow strength diminishes as the cavern enlarges, resulting in more vortices and a weaker flow intensity. From the third to the fifth stages, vortices consistently form in the dissolution section, merging with the previous vortices as expansion continues. After each stage, the vortex center shifts to the dissolution section’s adjacent cavity, away from the outlet. By the fifth stage, with the outlet over 200 m from the injection point, gravitational differentiation surpasses the forced convection’s impact on the brine, causing the denser brine to sink and the lighter brine to rise, forming a new weaker vortex.

The vortex cyclically circulates fluid within the cavity, enhancing the brine convective diffusion speed. This process replaces high-concentration brine near the cavity wall with lower concentration brine, leading to faster leaching rates and higher construction efficiency. Cavity expansion is primarily driven by dissolution and erosion under forced convection. Initially, erosion dominates in the dissolution cavern section, but dissolution becomes the primary factor over time. The expansion rate of the formed cavity section depends on the vortex strength, as erosion gives way to dissolution. As shown in Figure 16, the dissolution section expands the quickest, with expansion slowing in the sections closer to the outlet due to a diminished vortex intensity.

3.3. Analysis of Brine Concentration Distribution

The injected freshwater initially dissolves and erodes the salt wall, creating high-concentration brine in the early leaching stage, which then diffuses due to the high-speed fluid. Figure 17 illustrates the brine displacement concentration contours at the vertical section (X = 0 m) during the first stage, showing a relatively uniform brine concentration on the salt rock surface at the same leaching time, with higher concentrations in the reflux cavern than those in the dissolution cavern. As leaching progresses, the brine concentration in the salt cavern increases quickly, redistributing due to convection, diffusion, and vortex effects, with convection being predominant due to continuous high-speed freshwater injection. Due to the limited vertical extent of the horizontal cavern, the gravitational effects are minimal, making fluid erosion the leading factor and resulting in brine concentrations ranging from 60 to 120 g/L during the first stage.

Fluid erosion is initially predominant in the construction section, diminishing as leaching progresses, while dissolution takes precedence in the cavern formed during the first stage. At the onset of the second stage, the brine concentration decreases due to the low height of the dissolution cavern and the swift overall brine flow, as shown in Figure 18. The brine concentration in the dissolution cavern is notably lower than that in the cavern formed in the first stage. With cavity expansion, the vortices from both the sections gradually merge, leading to a gradual equalization of the brine concentration. In the reflux cavern, the brine concentration significantly rises compared to that in the first stage, reaching up to 270 g/L.

Figure 19 illustrates the brine concentration contours at the vertical sections (X = 0 m) for each stage, revealing a gradual increase in the overall brine concentration within the cavern, with the concentrations nearing saturation by the final leaching stage. Initially, the highest brine concentrations are in the reflux cavern, while the lowest ones are in the dissolution section. By the fourth and fifth stages, the outlet concentrations exceed those in the reflux cavern, yet remain the lowest in the dissolution section. Due to stable flow field distribution, the brine concentration is relatively uniform across stages, with even distribution along the cavern walls and stable cavity expansion. At the fifth stage, the injection section exhibits the lowest concentration, increasing towards the outlet, as the concentrations approach saturation, indicating high construction efficiency.

3.4. Analysis of Cavern Development

The cavern contours at the vertical section (X = 0 m) and the different cross sections are depicted after each leaching stage in Figure 20 and Figure 21. Initially, the solution cavity reaches a height of approximately 12 m. By the second stage, it expands upwards by about 7 m, roughly 2/3 of its initial height. In the third stage, the previously formed cavities continue to expand upwards, albeit at a diminishing rate. Continuing leaching results in no outward expansion for the first-stage cavity, reaching a maximum height influenced by the injection flow. Thus, the previously formed cavities expand at a gradually decreasing rate in the subsequent stages.

The results indicate that the solution cavity’s height in the first stage surpasses that in the subsequent four stages, with the latter stages maintaining a consistent height. The reflux cavity measures approximately 9 m in height and 15 m in length. Despite the variations in brine concentration and fluid flow, the reflux cavities at the inlet section retain a consistent height and shape after each leaching stage.

As shown in Figure 22, the outlet brine concentration and cavern volume evolve with the leaching time. In the first stage, the outlet brine concentration quickly increases, with the rate of rise decelerating over time. The concentration experiences a significant initial surge during the first leaching stage, followed by a sharp decrease, and then a gradual increase in the subsequent stages. The fluctuations in brine concentration become progressively milder. At start of each stage, there is a slight drop in concentration, which then slowly climbs, with the later stages showing higher concentrations than the earlier ones. As the cavern expands, the outlet brine concentration eventually stabilizes. Once the leaching rate stabilizes, the outlet concentration also gradually stabilizes, indicating that the changes in outlet concentration directly influence the leaching rate. For the cavern volume, despite rapid expansion initially, the cavern volume grows slowly due to the small contact area between the salt rock wall and brine.

4. Sensitivity Analysis and Discussions

A series of simulations with varied parameters were performed to optimize the construction of multi-step horizontal caverns. The key parameters include the borehole length, the step time, the injection flow rate, the brine concentration, and the step distance. In the simulation, the step distance, the step time, the injection brine concentration, and the flow rate were consistent during the steps. The coupled flow model kept the borehole length and the step time constant. The injection flow rate and the step time were correlated, limited by the salt formation’s thickness. The injection brine concentration was adjusted to recycle the discharged brine and manage the upward dissolution rates. Thus, this discussion focuses on the injection flow rate, the brine concentration, and the step distance, detailed in Table 2.

4.1. The Impact Analysis of Injection Flow Rate

Under varying flow rates, the outlet brine concentration patterns tend to be consistent, initially surging, and then gradually stabilizing. Higher injection flow rates yield lower brine concentrations with less fluctuation, as shown in Figure 23. The cavern volume exponentially increases with higher flow rates, leading to greater cavern heights and volumes, as shown in Figure 24. Cavern expansion enlarges the contact area between the salt rock wall and brine, thus enhancing the leaching rate.

Figure 25 illustrates that the final shape of the cavern formed under various injection flow rates remains largely consistent. Higher flow rates result in taller dissolution and reflux caverns, with uniform increases in the cavern contour height corresponding to uniform increases in the flow rate. The contour of the cavern shaped by corrosion is relatively flat, with an increase from 100 to 300 m³/h in flow rate raising the cavern contour height by about 3 m, suggesting that high flow rates are not advisable at the first leaching stage. Higher injection flow rates result in taller dissolution caverns and the quicker expansion of the previously formed caverns, as shown in Figure 26. In the later leaching stages, these caverns continue to dissolve and expand upwards. Increasing the injection flow rate enhances forced convection, intensifies fluid flow in the dissolution cavern, speeds up convective diffusion, and promotes continuous cavern expansion. Thus, increasing the injection flow rate and increasing the dissolution time in the subsequent stages is necessary to make the cavern smoother.

4.2. The Impact Analysis of Brine Concentration on Injection Fluid

Figure 27 shows that higher injection brine concentrations result in a narrower fluctuation range of outlet brine concentration in the first stage, starting at the injection concentration level. The pattern of brine concentration change, featuring an initial rapid increase followed by a gradual approach to saturation, remains similar across the different injection concentrations. Higher injection concentrations lead to higher overall outlet concentrations with smaller fluctuation ranges. With low injection concentrations, the outlet concentrations initially drop sharply at start of each stage, and then increase smoothly. This phenomenon will gradually disappear as the injection concentration rises. The cavern volume grows exponentially under various injection brine concentrations, mirroring the control group, as shown in Figure 28. Higher injection concentrations slow cavern volume expansion, with uniform increases in concentration leading to uniform decreases in the final volume. To maximize the cavern volume, injecting as much freshwater as possible, while minimizing the recycling of discharged brine, is advisable.

As shown in Figure 29, the injection brine concentration significantly impacts cavern expansion and the leaching rate in the first stage. Higher concentrations result in reduced heights and lengths of both the dissolution and return cavities, with the cavern height decreasing more rapidly as the concentration increases uniformly. At a concentration of 200 g/L, the height of dissolution cavity was only half compared to that when using freshwater. The dissolution cavity contour remains relatively smooth under varying concentrations. Higher concentrations lead to slower leaching rates and cavern expansion, reducing the height of the dissolution cavities, as shown in Figure 30. The cavern shape remains largely unchanged with concentration adjustments at the injection point. The effect of injection concentration on expansion significantly surpasses the temperature influence, highlighting the substantial impact of both the injection flow rate and the brine concentration on cavern development.

4.3. The Impact Analysis of Step Distance

With an increasing step distance, the trend and values of the brine discharge concentration closely align with the control group, as shown in Figure 31. The outlet brine concentration first rapidly increases, and then rises steadily to saturation. As the step distance lengthens in each stage, the overall concentration of discharged brine increases, the fluctuations lessen, and it nears saturation. At a 100 m step distance, the brine concentration in the last three stages stabilizes and achieves saturation. Thus, reaching saturation occurs earlier with larger step distances, highlighting the importance of optimizing the step distance in practical leaching processes. As the step distance grows, the increase in cavern volume is not obvious, as shown in Figure 32. While a larger step distance theoretically expands the contact area between the salt wall and brine, potentially enhancing the leaching rate and the cavern volume, in practice, the step distance does not significantly impact the leaching rate or the final cavern volume. Instead, the cavern volume is primarily influenced by the injection flow rate.

At a 50 m step distance, the cavern is relatively flat with a uniform height of approximately 12.5 m after the first stage. The height and length of the reflux cavern decrease as step distance increases, as shown in Figure 33, due to more fluid reaching the distal end and less returning to the reflux cavern. The cavern contour tends to be higher in the middle and lower areas of the outlet with an increased step distance, a trend that becomes more pronounced with larger distances, as shown in Figure 34. At certain intervals, the highest part of the cavern is just in front of the injection point, tapering off towards the outlet. A larger step distance results in a greater disparity between the maximum and minimum cavern heights. A step distance of 50 m results in a flat, smooth cavern, while at 100 m, the height difference reaches 4 m, with discrepancies forming during the later dissolution stages.

Cavern shape consistency across stages is maintained with specific step distances, injection flow rates, brine concentrations, and temperature differences. A shorter step distance under a given flow rate results in more irregular cavern shapes, indicating the minimum necessary step distance. For an injection flow rate of 200 m³/h, this minimum is 60 m, below which the cavern shape becomes uneven. The optimal step distance for achieving uniform dissolution lies between 70 and 80 m. Increasing the step distance requires extending the dissolution time of each stage to attain the desired cavern heights. In practice, a comprehensive assessment of various parameters is essential for optimizing the design parameters to maximize the volume and minimize the costs.

5. Conclusions

The main results and conclusions obtained in this study are summarized as follows:

• (1). A 3D multi-physics coupled flow model is established for constructing multi-step horizontal salt caverns to simultaneously quantify the dissolution, cavity expansion, heat transfer, brine convection, turbulence flow, and diffusion effects. The results are consistent with field data of Volgograd, verifying the accuracy of the proposed model.

• (2). The influence of heat transfer and turbulence flow is considered. The heat transfer accelerates the transportation and diffusion of brine. Temperature has a great influence on the dissolution rate of salt rock. The higher the temperature at the wall is, the faster the dissolution rate is. The direction of forced thermal convection is perpendicular to the direction of fluid flow and consistent with the direction of the brine diffusion in the cavity.

• (3). The expansion of the new cavity is affected by both erosion and dissolution, and the expansion of old formed cavity is mainly affected by dissolution, and the expansion rate depends on the strength of the vortex. Forced convection has a significant effect on brine concentration distribution.

• (4). The previously formed cavity will continue to expand at a decreasing rate gradually as dissolution proceeds in the next stage. In the first stage, cavity expansion outward measures approximately 12 m, 7 m, 4 m, 2 m, and 0.5 m, respectively. The reflux cavities at the inlet section have the same shape. The solution cavity has the same height at 9 m at the last four stage.

• (5). Increasing the injection flow rate results in a larger cavern volume and efficiency. Increasing the injection concentration results in a smaller cavern volume and efficiency. The step distance hardly affects the cavern volume, but it does affect the cavern shape. A larger injection flow rate, a lower injection concentration, and a middle step distance are preferable in the field.

Author Contributions

Conceptualization, Z.W.; methodology, J.W. (Jinchao Wang). and Q.Z.; software, J.W. (Jinchao Wang); validation, J.W. (Jinchao Wang); formal analysis, J.W. (Jinchao Wang); investigation, J.W. (Jinchao Wang), J.W. (Jun Wang) and B.L.; resources, Z.W.; data curation, J.W. (Jinchao Wang), J.W. (Jun Wang) and B.L.; writing—original draft preparation, J.W. (Jinchao Wang); writing—review and editing, Q.Z.; visualization, J.W. (Jinchao Wang); supervision, Z.W.; project administration, Z.W.; funding acquisition, Z.W. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author due to privacy.

Acknowledgments

The authors wish to acknowledge the financial support from the PetroChina Research Institute of Petroleum Exploration and Expansion (PIPED.CN-2f-2019-JS-359) and China Scholarship Council (CSC202306440111).

Conflicts of Interest

The authors declare no conflicts of interest.

List of Symbols

$c$	Concentration of brine
${\mathrm{C}}_{dis}$	Coefficient of dissolving rate
${\mathrm{C}}_p$	Fluid heat capacity at constant pressure
$c_{sat}$	Concentration of saturated brine
$C_{1\epsilon }$	Empirical coefficient measured in experiment
$C_{2\epsilon }$	Empirical coefficient measured in experiment
$C_{\mu }$	Experimental coefficient
D	Diffusion coefficient for NaCl
$G_k$	Turbulence kinetic energy due to the mean velocity gradients
$G_b$	Turbulence kinetic energy due to buoyancy
$k$	Fluid thermal conductivity
$M$	Molar mass of sodium chloride
$p$	Pressure, Pa; ${\mu }_t$ is the turbulent viscosity
$p_0$	Pressure of the salt formation
$Q$	Heat source
T	Temperature of the salt cavern
T 0	Freezing temperature of water
$T_{\mathrm{in}}$	Inlet temperature
$T_{wall}$	Cavern wall temperature
$u_i$	Brine velocity component in the i direction
$u_j$	Brine velocity component in the j direction
$u_i^{\prime }$	Fluctuating velocity component in i direction
$u_j^{\prime }$	Fluctuating velocity component in the j direction
$u_g$	Velocity component in the gravity direction
$w_n$	Moving velocity of the interface in the normal direction
$Y_M$	Contribution of the fluctuating dilatation to the dissipation rate
$\mathit{\Phi }$	Dissolution flux at the salt–brine interface
$\rho$	Fluid density
$k$	Turbulence kinetic energy
$\epsilon$	Turbulence dissipation rate
${\sigma }_k$	Turbulent Prandtl number for kinetic energy
${\sigma }_{\epsilon }$	Turbulent Prandtl number for dissipation rate

Footnotes

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.

Figures and Tables

Figure 1. Schematic diagram of the consumption of natural gas in China.

Figure 2. Schematic diagram of a multi-step horizontal salt cavern.

Figure 3. A diagram of thermal convection and thermal conduction during solution mining. (The yellow arrows represent thermal conduction, and the blue arrows represent thermal convection.)

Figure 4. Schematic diagram of geometry and meshing results.

Figure 5. The process of the adaptive meshing in the calculation.

Figure 6. Flow chart of the simulation of the coupled model.

Figure 7. Data comparison of outlet brine concentration evolutions with leaching time.

Figure 8. Data comparison of cavern volume evolutions with leaching time.

Figure 9. Comparison of cavern side shapes [46].

Figure 10. Data comparison of concentration and cavern volume evolutions with leaching time.

Figure 11. Diagram of the thermal convection direction and fluid flow direction of brine in the cavity. (The blue arrows represent thermal convection, and the red arrows represent fluid flow.)

Figure 12. Diagram of the thermal convection direction and the concentration contour of brine in the cavity. (The blue arrows represent thermal convection, and the curves represent the concentration contour line).

Figure 13. The cavern shape under different temperatures after 10 days of leaching.

Figure 14. Contours of velocity magnitude at vertical section (X = 0 m) during first stage. (The red arrow line represents the streamline of brine).

Figure 15. Contours of velocity magnitude at vertical section (X = 0 m) during second stage (The red arrow line represents the streamline of the brine).

Figure 16. Contours of velocity magnitude at vertical section (X = 0 m) of each stage (The red arrow line represents the streamline of the brine).

Figure 17. Contours of brine displacing concentration at vertical section (X = 0 m) during first stage (rainbow lines represent flow velocity lines).

Figure 18. Contours of brine displacing concentration at vertical section (X = 0 m) during second stage (rainbow lines represent flow velocity lines).

Figure 19. Contours of brine displacing concentration at vertical section (X = 0 m) of each stage (black lines represent brine concentration contour lines).

Figure 20. Cavity shapes at vertical section (X = 0 m) of each stage.

Figure 21. Cavity shapes at different cross sections of each stage.

Figure 22. Outlet brine concentration and cavern volume evolutions with leaching time.

Figure 23. Outlet brine concentration evolution with leaching time under different flow rates.

Figure 24. Cavern volume evolution with leaching time under different flow rates.

Figure 25. First-stage cavern shapes in vertical section (X = 0 m) under different flow rates.

Figure 26. Cavity’s final shapes in vertical section (X = 0 m) under different flow rates.

Figure 27. Outlet brine concentration evolution with leaching time under different injection concentrations.

Figure 28. Cavern volume evolution with leaching time under different injection concentrations.

Figure 29. First-stage cavern shapes in vertical section (X = 0 m) under different injection concentrations.

Figure 30. Cavity’s final shapes in vertical section (X = 0 m) under different injection concentrations.

Figure 31. Outlet brine concentration evolution with leaching time under different injection intervals.

Figure 32. Cavern volume evolution with leaching time under different injection intervals.

Figure 33. First stage cavern shapes in vertical section (X = 0 m) under different injection intervals.

Figure 34. Cavity’s final shapes in vertical section (X = 0 m) under different injection intervals.

Table 1

Basic parameters of the coupled model.

Parameters	Sign	Value	Units
Diffusion coefficient	D	1.6 × 10−9	m2/s
Molar mass of NaCl	M	58.5	g/mol
Density of NaCl	cs	2165	kg/m3
Saturated brine concentration	csat	318.4	g/L
Viscosity of saturated brine	μ0	0.001	Pa·s
Density of saturated brine	ρsat	1200	kg/m3
Dissolving rate of salt rock	Cdis	3.5 × 10−5	m/s
Depth of the formation	H	1500	m
T of the injection fluid	Tin	293.15	K
Initial T of the salt wall	Twall	333.15	K
Pressure of salt formation	P0	15	MPa
Inner tube diameters	din	177.8	mm
Outer tube diameters	dout	244.5	mm
Length of the solution cavity	Ls	50	m
Length of the reflux cavity	Lr	50	m
Injection flow rate	Q	200	m3/h
Leach time of every stage	T	50	day

Table 2

The technological parameters of the three groups of simulations.

Group	Borehole Length (m)	Step Time (Days)m	Flow Rate (m3/h)	Concentration (g/L)	Step Distance (m)
A	250	50	100, 150, 200, 250, 300	0	50
B	250	50	200	0, 50, 100, 150, 200	50
C	250	50	200	0	50, 60, 70, 80, 90, 100

References

1. Deng, J.Q.; Yang, Q.; Liu, Y.R. Time-dependent behavior and stability evaluation of gas storage caverns in salt rock based on deformation reinforcement theory. Tunneling Undergr. Space Technol.; 2014; 42, pp. 277-292. [DOI: https://dx.doi.org/10.1016/j.tust.2014.03.014]

2. Chen, X.S.; Li, Y.P.; Liu, W.; Ma, H.; Ma, J.; Shi, X.; Yang, C. Study on Sealing Failure of Wellbore in Bedded Salt Cavern Gas Storage. Rock Mech. Rock Eng.; 2019; 52, pp. 215-228. [DOI: https://dx.doi.org/10.1007/s00603-018-1571-5]

3. Yang, C.H.; Wang, T.T.; Li, Y.P.; Yang, H.; Li, J.; Xu, B.; Yang, Y.; Daemen, J.J.K. Feasibility analysis of using abandoned salt caverns for large-scale underground energy storage in China. Appl. Energy; 2015; 127, pp. 467-481. [DOI: https://dx.doi.org/10.1016/j.apenergy.2014.07.048]

4. Tao, C.; Xilin, S.; Jinlong, L. Simulation test on the horizontal-well leaching of salt-cavern gas storage. Oil Gas Storage Transp.; 2019; 38, pp. 1257-1264.

5. Atoyebi, T.M. The Preferred Natural Gas Conservation Option: Underground Storage of Natural Gas. Proceedings of the Nigeria Annual International Conference and Exhibition; Calabar, Nigeria, 31 July–7 August 2010; SPE-136984-MS.

6. Zhang, J.D.; Tan, Y.F.; Zhang, T.T.; Yu, K.; Wang, X.; Zhao, Q. Natural Gas Market and Underground Gas Storage Development in China. J. Energy Storage; 2020; 29, 101338. [DOI: https://dx.doi.org/10.1016/j.est.2020.101338]

7. Zhang, N.; Liu, W.; Zhang, Y.; Shan, P.; Shi, X. Microscopic pore structure of surrounding rock for underground strategic petroleum reserve (SPR) caverns in bedded rock salt. Energies; 2020; 13, 1565. [DOI: https://dx.doi.org/10.3390/en13071565]

8. Zhang, N.; Ma, L.J.; Wang, M.Y.; Zhang, Q.; Li, J.; Fan, P. Comprehensive risk evaluation of underground energy storage caverns in bedded rock salt. J. Loss Prev. Process Ind.; 2016; 45, pp. 264-276. [DOI: https://dx.doi.org/10.1016/j.jlp.2016.10.016]

9. Parkes, D.; Evans, D.J.; Williamson, P.; Williams, J.D.O. Estimating available salt volume for potential CAES development: A case study using the Northwich Halite of the Cheshire Basin. J. Energy Storage; 2018; 18, pp. 50-61. [DOI: https://dx.doi.org/10.1016/j.est.2018.04.019]

10. Peng, T.J.; Wan, J.F.; Liu, W.; Li, J.; Xia, Y.; Yuan, G.; Jurado, M.J.; Fu, P.; He, Y.; Liu, H. Choice of hydrogen energy storage in salt caverns and horizontal cavern construction technology. J. Energy Storage; 2023; 60, 106489. [DOI: https://dx.doi.org/10.1016/j.est.2022.106489]

11. Raju, M.; Khaitan, S.K. Modeling and simulation of compressed air storage in caverns: A case study of the Huntorf plant. Appl. Energy; 2012; 89, pp. 474-481. [DOI: https://dx.doi.org/10.1016/j.apenergy.2011.08.019]

12. Jean, P.M.; Li, Z.; Jacob, B. Dynamic modeling of compressed gas energy storage to complement renewable wind power intermittency. Int. J. Hydrog. Energy; 2013; 38, pp. 7867-7880.

13. Tackie-Otoo, B.N.; Haq, M.B. A comprehensive review on geo-storage of H₂ in salt caverns: Prospect and research advances. Fuel; 2024; 356, 129609. [DOI: https://dx.doi.org/10.1016/j.fuel.2023.129609]

14. Zhang, Z.X.; Liu, W.; Guo, Q.; Duan, X.; Li, Y.; Wang, T. Tightness evaluation and countermeasures for hydrogen storage salt cavern contains various lithological interlayers. J. Energy Storage; 2022; 50, 104454. [DOI: https://dx.doi.org/10.1016/j.est.2022.104454]

15. Tietze, V.; Stolten, D. Comparison of hydrogen and methane storage by means of a thermodynamic analysis. Int. J. Hydrog. Energy; 2015; 40, pp. 11530-11537. [DOI: https://dx.doi.org/10.1016/j.ijhydene.2015.04.154]

16. Mei, S.W.; Gong, M.Q.; Qin, G.L.; Tian, F.; Xue, X.; Li, R. Advanced Adiabatic Compressed Air Energy Storage System with Salt Cavern Air Storage and Its Application Prospects. Power Syst. Technol.; 2017; 41, pp. 3392-3399.

17. Ding, G.S.; Li, C.; Wang, J.M.; Xu, H.; Zheng, Y.; Wanyan, Q.; Zhao, Y. The status quo and technical development direction of underground gas storages in China. Nat. Gas Ind.; 2015; 35, pp. 107-112. [DOI: https://dx.doi.org/10.1016/j.ngib.2015.12.007]

18. Zhang, B.; Lv, B.; Wu, Y.H.; Cui, L.H.; Zhou, S.Q. Development and trend of salt-cavern gas storage in domestic and abroad. China Well Rock Salt; 2021; 52, pp. 21-24.

19. Ding, G.; Ding, Y.; Li, Y.; Tang, L.; Wu, Z.; Wanyan, Q.; Xu, H.; Wang, Y. Prospects of underground gas storage in China under the strategy of carbon neutrality. Oil Gas Storage Transp.; 2022; 41, pp. 1-9.

20. Yuan, G.J.; Wan, J.F.; Li, J.C.; Li, G.; Xia, Y.; Ban, F.; Zhang, H.; Jurado, M.J.; Peng, T.; Liu, W. Stability analysis of a typical two-well-horizontal saddle-shaped salt cavern. J. Energy Storage; 2021; 40, 102763. [DOI: https://dx.doi.org/10.1016/j.est.2021.102763]

21. Moditis, K.; Paidoussis, M.; Ratigan, J. Dynamics of a partially confined, discharging, cantilever pipe with reverse external flow. J. Fluids Struct.; 2016; 63, pp. 120-139. [DOI: https://dx.doi.org/10.1016/j.jfluidstructs.2016.03.002]

22. Wan, J.F.; Peng, T.J.; Yuan, G.J.; Ban, F.; Jurado, M.J.; Xia, Y. Influence of tubing/oil-blanket lifting on construction and geometries of two-well-horizontal salt caverns. Tunneling Undergr. Space Technol.; 2021; 108, 103688. [DOI: https://dx.doi.org/10.1016/j.tust.2020.103688]

23. Chang, X.N. Study on Subsurface Characteristics for Salt Deposits in China and Geological Evaluation on Salt Cavern Storage. Ph.D. Thesis; China University of Geosciences: Beijing, China, 2014; 5.

24. Xie, N.; Li, W.J. Thermal damage mechanism of salt rock in underground salt cavern gas storage. Acta Pet. Sin.; 2019; 40, pp. 357-369.

25. Wei, G.Q.; Zheng, Y.L.; Qiu, X.S.; Sun, J.; Shi, L.; Lai, X. Geological theory and application of underground gas storage in China. Acta Pet. Sin.; 2019; 40, pp. 1519-1530.

26. Yang, J.; Li, H.; Yang, C.H.; Li, Y.; Wang, T.; Shi, X.; Han, Y. Physical simulation of flow field and construction process of horizontal salt cavern for natural gas storage. J. Nat. Gas Sci. Eng.; 2020; 80, 103527. [DOI: https://dx.doi.org/10.1016/j.jngse.2020.103527]

27. Kunstman, A.S.; Urbanczyk, K.M. Modelling of Horizontal Cavern Leaching: Main Aspects and Perspectives. Proceedings of the SMRI Fall 1995 Meeting; San Antonio, TX, USA, 22–25 October 1995.

28. Liu, W.; Jiang, D.Y.; Chen, J.; Daemen, J.J.K.; Tang, K.; Wu, F. Comprehensive feasibility study of two-well-horizontal caverns for natural gas storage in thinly-bedded salt rocks in China. Energy; 2018; 143, pp. 1006-1019. [DOI: https://dx.doi.org/10.1016/j.energy.2017.10.126]

29. Li, Y.; Liu, W.; Yang, C.; Daemen, J.J.K. Experimental investigation of mechanical behavior of bedded rock salt containing inclined interlayer. Int. J. Rock Mech. Min. Sci.; 2014; 69, pp. 39-49. [DOI: https://dx.doi.org/10.1016/j.ijrmms.2014.03.006]

30. Li, Q.; Ning, Z.; Liu, J.; Xu, W.; Zhan, L.; Liu, J.; Chen, Y.; Shi, X.; Chen, X.; Li, J. Stability and economic evaluation of multi-step horizontal salt caverns with different step distances in bedded salt formations. J. Energy Storage; 2023; 57, 106192. [DOI: https://dx.doi.org/10.1016/j.est.2022.106192]

31. Urbańczyk, K. Leaching in thin-bedded deposits: Horizontal cavern modeling (Salt-cavern leaching process, Part V). Salt Rev.; 2022; 16, pp. 67-91.

32. Zeng, Q.S.; Wang, Z.M.; Wang, J.C.; Wanyan, Q.; Ding, G.; Li, K. Modeling horizontal salt cavern leaching in bedded salt formations. SPE J.; 2022; 27, pp. 1367-1378. [DOI: https://dx.doi.org/10.2118/209197-PA]

33. Thoms, R.L.; Gehle, R.M. Feasibility of controlled solution mining from horizontal wells. Proceedings of the SMRI Fall 1993 Technical Meeting at Louisiana; Lafayette, IN, USA, 24–28 October 1993.

34. Saberian, A. A preliminary model for horizontal well leaching. Proceedings of the SMRI Fall 1995 Technical Meeting; San Antonio, TX, USA, 22–25 October 1995.

35. Russo, A.J. HORSMIC solution mining code for horizontal caverns. Proceedings of the SMRI Fall 1995 Technical Meeting; San Antonio, TX, USA, 22–25 October 1995.

36. Kang, Y.P.; Jiao, Y.J.; Wang, J.F.; Ling, Z.; Wu, Z. Experiment on horizontal leaching characteristics of salt cavern gas storage. Oil Gas Storage Transp.; 2022; 41, pp. 1061-1068. (In Chinese)

37. Wan, J.F.; Peng, T.J.; Shen, R.C.; Jurado, M.J. Numerical model and program development of TWH salt cavern construction for UGS. J. Pet. Sci. Eng.; 2019; 179, pp. 930-940. [DOI: https://dx.doi.org/10.1016/j.petrol.2019.04.028]

38. Li, J.L.; Shi, X.L.; Zhang, S. Construction modeling and parameter optimization of multi-step horizontal energy storage salt caverns. Energy; 2020; 203, 117840. [DOI: https://dx.doi.org/10.1016/j.energy.2020.117840]

39. Li, J.L.; Yang, C.H.; Shi, X.L.; Xu, W.; Li, Y.; Daemen, J.J.K. Construction modeling and shape prediction of horizontal salt caverns for gas/oil storage in bedded salt. J. Pet. Sci. Eng.; 2020; 190, 107058. [DOI: https://dx.doi.org/10.1016/j.petrol.2020.107058]

40. Wang, J.C.; Wang, Z.M.; Zeng, Q.S.; Ding, G.; Li, K.; Wanyan, Q.; Wang, Y. Simulation of Flow Field of Solution Mining Salt Cavities for Underground Gas Storage. J. Energy Resour. Technol.; 2023; 145, 022001. [DOI: https://dx.doi.org/10.1115/1.4054913]

41. Wu, G.; Xiao, C.F.; Qiu, X.D. Study on the dissolution rate of rock salt. Chem. Min. Technol.; 1992; 21, pp. 18-22.

42. Xiao, C.F.; Yang, Y.K.; Wu, G.; Li, C.P.; Qiu, X.D. The dissolution characteristics of rock salt and mass transmission process. J. Chongqing Univ.; 1993; 16, pp. 51-57.

43. Wang, C.R. Experimental Study of the Influencing Factors on Dissolution Rate of Rock Salt; Chongqing University: Chongqing, China, 2009.

44. Souza, W.R.; Voss, C.I. Analysis of an anisotropic coastal aquifer system using variable-density flow and solute transport simulation. J. Hydrol.; 1987; 92, pp. 17-41. [DOI: https://dx.doi.org/10.1016/0022-1694(87)90087-4]

45. Shin, C.B.; Economou, D.J. Effect of transport and reaction on the shape evolution of cavities during wet chemical etching. J. Electrochem. Soc.; 1989; 136, 1997. [DOI: https://dx.doi.org/10.1149/1.2097124]

46. Igoshin, A.I.; Kazaryan, V.A.; Pozdnyakov, A.G.; Salokhin, V.I.; Vasjuta, Y.S. Creation of Tunnel Cavities in Rock Salt Beds of Small Thickness. Proceedings of the SMRI Spring 1997 Meeting; Krakow, Poland, 11–14 May 1997.