Stability analysis of artificial dam in coal mine underground water reservoir based on the hydro-mechanical damage model

Abstract A numerical hydro-mechanical damage model is proposed in this paper to investigate the stability and failure characteristics of artificial dam in coal mine underground water reservoir. The localized damage model for the coupled pore pressure and rock failure of a stressed rock is developed, and a local damage constitutive law is employed to describe its stress-strain relationship and fracture propagation. The Mohr-Coulomb criterion and maximum tension criterion are used to Judge the damage and failure of rock. Experimental data for rock samples are used to validate the proposed numerical model, and it accurately replicates the stress-strain curves and failure pattern compared with the experimental results. We then explore the influence of water pressures, cutting depth and dam strengths of artificial dam in coal mine underground water reservoir, and the simulation results show that the dam strength is the main factor affecting the strength of artificial dam. The modeling described herein is expected to assist in the management and optimization of underground water reservoirs.


Introduction
The western region of China (Shanxi, Shaanxi, Mongolia, and Gansu) is the primary coal producing area, accounting for almost two-thirds of the country's proven reserves and production, but only 3.9% of the country's water resources (Gu 2015;Zhang et al. 2021).Therefore, the technological system of coal mine underground water reservoir is ingeniously proposed by Shenhua Group to circumvent the issues of outside evaporation loss, surface water treatment plant construction, and operating cost.Underground water reservoirs have evolved into a crucial technical approach for water resource protection in the process of coal mining in western China (Gu et al. 2016).
In contrast to surface dams, underground water reservoir dams are subjected not only to osmotic water pressure but also to high abutment pressure resulting from initial rock stress and mining disturbances (Yao et al. 2020).Under the action of hydromechanical coupling, the deterioration and instability of underground reservoir dam are the result of the failure and deterioration of the dam mechanical properties (Zhang et al. 2021).On the one hand, the dam is damaged and broken by the repeated action of in-situ rock stress and high abutment pressure caused by coal mining; on the other hand, groundwater accelerates the expansion and evolution of fractures in surrounding rock via pore pressure and water physicochemical interactions, thereby aggravating the damage and deterioration of the mechanical properties of surrounding rock (Wang et al. 2020;Andr es et al. 2017).
Numerous researchers have conducted the seepage characteristics test of rock using an advanced electro-hydraulic servo-control rigid testing machine and obtained the permeability variation law of the entire process of rock deformation and failure under specific conditions.Some scholars improved the conventional test equipment (Souley et al. 2001;Wang and Park 2002;Aker et al. 2014), established a complete rock seepage and creep coupling test system, and conducted experimental research on the creep and permeability characteristics of rock under the coupling effect of stress-seepage.In order to study the seepage characteristics of rock which could not be completed by experimental equipment, the hydro-mechanical coupling model and numerical simulation of rock have been explored by many scholars (Xu and Yang 2016;Xu et al. 2013;Xu et al. 2012;Yang et al. 2014;Zhu et al. 2014).However, there are few studies on the properties of underground water reservoir rock mass, taking into account the coupling effect of seepage and stress (Zhang et al. 2015;Hu et al. 2016;Zhang et al. 2016;Hu et al. 2017;Hu et al. 2017;Yang et al. 2018;Zhang et al. 2022;Luo et al. 2021;Scheperboer et al. 2022;Krzaczek et al. 2020;Zhao et al. 2022;Hu et al. 2017;Zhao et al. 2022).
In view of the importance of understanding rock mechanical behavior in underground reservoirs (Bi et al. 2020;Xu et al. 2018;Zhou et al. 2020;Xu and Yang 2016), a hydro-mechanical creep model is presented in this paper to describe the deformation and failure characteristics of rock samples in underground water reservoirs.The hydro-mechanical coupled time-dependent model is first formulated and then validated using experimental data.The results of the brittle creep simulations with varying confining and pore pressures are then discussed in further details.Then, the influence of water pressures and dam strengths of artificial dam are explored through a simplified 3D damage numerical model in coal mine underground water reservoir.

Formulation of constitutive model
This section proposes a hydro-mechanical creep mode for the coupled pore pressure and rock failure problems to characterize the deformation and fracture of rock samples.The coupled effect of rock deformation and water flow must be taken into consideration while developing a model that can be used to investigate deformation and failure under varying constant deviatoric stress and pore pressure conditions.

Hydro-mechanical model
Rock is a kind of porous medium material, and it has obvious nonlinear characteristics.However, for the sake of simplification, assuming that rocks are elastic materials and ignoring its plastic deformation in this paper.The seepage characteristic has a great influence on its deformation and liquid flow, Darcy's law as a well-developed and practical constitutive relation is used to describe the fluid flow in a porous medium as follows: where P is pore fluid pressure, z is the vertical coordinate, k is the coefficient of permeability, q l is the fluid density (kg/m 3 ), and g is the acceleration due to surface gravity (m/s 2 ).Substituting Eq. ( 1) into the fluid conservation equation yields (Biot 1956): where K S is the effective bulk modulus, e V is volumetric strain, and K 0 ðK 0 ¼ 2Gð1 þ lÞ=3ð1 À 2lÞÞ is the drained bulk modulus, u is the porosity, and b l denotes the bulk modulus of fluid.
Initially rock as a kind of porous medium material is assumed elastic, and a generalized Hooke's law is defined as its constitutive relationship (Zhu et al. 2014;Xu et al. 2006).Therefore, the ultimate governing equation of a deformed porous medium can be derived from the static stress equilibrium equation, the Cauchy strain tensor equation, and the constitutive stress-strain relation.It can be expressed in displacement form as follows: where the constants G and k are called the Lam e constants, u i is the components of displacement, P is the pore pressure, and the parameter a ( 1) is the Biot's coefficient, which depends on the compressibility of the constituents and can be defined as a ¼ 1 À K 0 =K s :

Damage model
The mechanical properties of rock are very complex, with heterogeneity, nonlinearity and non-continuity.In order to better reflect the damage and failure process of rock, a Weibull distribution is adopted in this study is used to assign a random distribution of the mechanical parameters to reflect the material heterogeneity at the mesoscale.The Weibull distribution function is expressed in the following form: where v represents the homogeneity index, the larger the value of homogeneity index, the more homogeneous the rock; u 0 denotes the average element scale parameter, and u a represents the scale parameter of an individual element.
When the stress of rock subjected to confining pressure and pore pressure at each loading step exceeds a damage threshold, the damage in tension or shear is judged according to the maximum tensile stress criterion and Mohr-Coulomb criteria (Zhu et al. 2014;Xu et al. 2006;Bi et al. 2020;Xu et al. 2018), respectively, as expressed by: (5) where r 1 and r 3 are the maximum and minimum principal stresses, respectively, r t0 is the uniaxial tensile strength, r c0 is the uniaxial compressive strength, u is the internal friction angle, and F 1 and F 2 are damage threshold functions.
If the element is determined to damage in tensile or shear, the damage variable D can be calculated as: where e t0 and e c0 are the maximum uniaxial tensile strain and compressive strain at the elastic limit, respectively; r is a constant with a value of 2.0 in this paper; e 1 is the major principal strain at each loading step, e ¼ À ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi , where e 1 , e 2 and e 3 are the three principal strains, hi stands for the positive part of a scalar, and hxi ¼ ðx þ jxjÞ=2: The elastic damage constitutive law (Xu et al. 2006) for an element under uniaxial compression and tension conditions is illustrated in Figure 1.
If the element is damaged, the elastic modulus of the element will be reduced monotonically as damage evolves according to the linear elastic damage principle.It describes the stress-strain relationship of an element subjected to an instantaneous loading as (Zhou et al. 2020): where E and E 0 are the damaged and undamaged elastic modulus, respectively, and D is the isotropic damage variable.Regardless of the type of damage, the damage variable, D, is always between 0 and 1 and is irreversible once damage occurred.Negative numbers indicate tensile damage, whereas positive numbers represent shear damage, such that the two types of damage, i.e. tensile and shear damage modes, can be distinguished in the post-processing figures.
This paper documents an iterative numerical procedure that is based on the finite element method.A flow chart of the model is shown in Figure 2 to clarify the implementation of the numerical model.In each loading step, stresses and strains for each element are calculated according to Eq. (1-3) by COMSOL Multiphysics (Xu and Yang 2016).Then, the maximum tensile stress criterion and Mohr-Coulomb criteria are used to judge whether an element is damaged, the tensile stress criterion (Eq.( 5)) is always used whether the element is damaged in tension, and the Mohr-Coulomb criterion (Eq.( 6)) is used to judge whether the element is damaged in shear.When an element is damaged due to loading, the damage variable is calculated according to Eq. ( 7) and a degradation of the damaged elements occurs based on the elastic damage constitutive law (Eq.( 8)).With ongoing loading damage evolution occurs until complete failure and simultaneous stop of the calculation procedure.

Validation of the model
In this paper, the floor rock mass of an underground water reservoir serves as the primary research object.As the floor rock mass of an underground water reservoir is predominantly composed of sandstone, we begin by validating our model using previously published data (Okubo et al. 2010) on sandstone.A series of biaxial compression tests at varying confining pressures of 5, 10, and 20 MPa are numerically conducted to extract the mesoscale input parameters from the macroscopic physicalmechanical properties of sandstone.The sandstone utilized in the numerical tests has dimensions of 100 mm in height and 50 mm in width.The Young's modulus and Poisson's ratio in biaxial compression tests are 20 GPa and 0.3, respectively.A series of predetermined displacement increments of 0.001 mm/step is applied to both ends of the numerical samples at constant confining pressures of 5, 10, and 20 MPa and a   pore pressure of 3 MPa.The mechanical properties of the numerical samples are listed in Table 1.
Figure 3 depicts the stress-strain curves and failure pattern predicted by the numerical model, as well as the experimental stress-strain curves and failure mode observed in laboratory tests.As demonstrated in Figure 3, the simulated stress-strain curves and failure pattern are in reasonable agreement with the experimental curves.Moreover, the ultimate failure patterns of numerical biaxial compression tests reveal that as confining pressure increases, compressive failure in the macroscopic crack area increases while tensile failure decreases.

Effect of water pressures on artificial dam
In order to study the influence of water pressures, cutting depth and dam strengths on the stability of artificial dam of underground reservoir, A simplified 3D numerical damage model is established based on the conditions of underground reservoir in Ningdong mining area.Using the established 3D numerical damage model, several calculation schemes are formulated for the water storage height and dam body strength, etc.The model has a strike length of 260 m, a strike width of 175 m and a height of 115 m, the numerical model of underground reservoir and artificial dam are shown in Figure 4.The tetrahedral mesh was used for the numerical simulation, and the numerical model had a total of 26,000 elements.Fixed boundary conditions are applied at the bottom of the model, the other four sides except the bottom and upper sides can move freely in the vertical direction, but cannot move in horizontal direction, and gravity is applied in the vertical direction.According to the geological data of Ningdong mining area, the mechanical parameters of each rock layer are listed in Table 2.
In order to study the influence of the stability of artificial dam body under different water pressures, eight numerical simulation schemes are designed based on the conditions of the underground reservoir in Ningdong mining area.The artificial dam is a plate, the strength of artificial dam is 25 MPa, dam thickness is 1 m, side cutting depth of the artificial dam is 0.5 m, roof cutting depth is 0.3 m, floor cutting depth is 0.2 m, head pressures are 5 m, 10 m, 15 m, 20 m, 25 m, 30 m, 35 m, 40 m.
Figure 5 shows the maximum principal stress on upstream surfaces and downstream surfaces of artificial dam body at different head pressures.As illustrated, the maximum tensile stress on upstream surfaces of artificial dam appears at the cutting edge around the dam body, and the maximum is at the cutting edge of the roof.In the downstream surfaces, the maximum principal stress appears in the bottom area, and the minimum is at the cutting edge of the bottom.As can be seen from Figure 5, the stress distribution under different head pressures is basically the same, but the stresses will increase with the increase of water pressures.Figure 6 shows the damage distribution of artificial dam on upstream surfaces and downstream surfaces at different head pressures.As can be seen from Figure 6, damage on the upstream surface first appears at 5 m head pressure, but no damage on the downstream surface, so the artificial dam body do not appear penetrating damage.With the increase of head pressure, the damage of downstream surface first appears at 25 m head pressure at  the bottom cut of artificial dam.It can be seen from the curve of the maximum displacement of artificial dam with the head pressures, the maximum displacement of the dam increases linearly with the water pressures, and the displacement of the dam increases about 6 Ã 10 À6 m for every 1 m increase in the head pressure.

Effect of cutting depth on artificial dam
In this section, the control variable method is used to change the cutting depth of the side groove, the top groove and the bottom groove respectively to study their influence on the stability of the dam body.The default strength of artificial dam is 25 MPa, dam thickness is 1 m.In order to study the influence of cutting depth on the stability of artificial dam, in the first set of simulations, the default size of floor cutting depth is 0.2 m and roof cutting depth is 0.3 m, side cutting depth of the artificial dam are 0.2 m, 0.4 m, 0.6 m, 0.8 m, respectively.In the second set of simulations, the default size of side cutting depth is 0.5 m, floor cutting depth and roof cutting depth of the artificial dam are 0.2 m, 0.4 m, 0.6 m, 0.8 m, respectively.Figure 7 depicts the distribution of damage on upstream surfaces and downstream surfaces of artificial dam with different side cutting depths.As can be seen from Figure 7, obvious damage occurs on the upstream surfaces and downstream surfaces as the side cutting depths are 0.2 m and 0.4 m.When the side cutting depth increases to 0.6 m and 0.8 m, the damage of upstream and downstream surfaces is within the acceptable range.In particular, there is only slight damage on the downstream surface, indicating that the stress and water pressure have not completely damaged the artificial dam.Therefore, the simulation results show that when the side cutting depth is 0.6 m, it can meet the safety requirements of the underground reservoir at 25 m water pressure.
Figure 8a and b shows the distribution of damage on upstream surfaces and downstream surfaces of artificial dam with different roof and floor cutting depths.As shown in Figure 8a and b, the damage on downstream surface can be ignored when the depth of roof and floor cutting depths are 0.6 m.As can be seen from the curve of maximum displacement of artificial dam with the cutting depth, the maximum displacement of artificial dam gradually slows down as the cutting depth increases.Taking different side cutting depths as an example, the change rates of the curve are 1.6%, 0.68%, 0.29% and 0.20% respectively, and the maximum displacement value also decreases from 1.06 Â 10 À4 m to 1.016 Â 10 À4 m.It can be seen from the simulation results that increasing the cutting depth has certain influence on reducing the displacement of artificial dam, and it can be concluded from the damage distribution and displacement curve that the cutting depth should be no less than 0.6 m.

Effect of dam strength on artificial dam
The artificial dam is mostly made in concrete, and the type of concrete used in the artificial dam determines the strength of the dam.Therefore, this section studies the effect of different concrete on the strength of artificial dam body through numerical simulation.In the numerical simulation, the head pressure is the limit water pressure 25 m, the cutting depth is 0.6 m, the strength of the artificial dam is 15, 20, 25, 30 and 35 MPa, respectively.According to Code for Design of Concrete Structures (GB50010-2010), the Poisson ratio of concrete is set to 0.2, the ratio of compressive to tensile strength is 10, Young's modulus is 2.2, 2.55, 2.8, 3.0 and 3.15 GPa, respectively.
As can be seen from Figure 9a, with the increase of the strength of the artificial dam, the peak value of the maximum principal stress of the artificial dam increases,  but the change rate decreases.Under different strength of artificial dam, the peak value of the maximum principal stress on upstream surface of artificial dam is distributed at the intersection of roof cutting and roadway (Figure 9b), and the peak value of the maximum principal stress on downstream surface of artificial dam is at the intersection of floor cutting and roadway.Figure 10 presents the damage distribution of artificial dam on upstream and downstream surface at different dam strength.When the dam strength is 15 MPa, severe damage occurs in the upstream and downstream surfaces of the artificial dam body.When the dam strength is 20 MPa, damage occurs at the intersection of roof cutting and roadway of the artificial dam, and slightly damage occurs at the intersection of floor cutting and roadway.There is no damage either on the upstream surface and downstream surface of the artificial dam when the dam strength is 25 and 30 MPa.Therefore, under the current engineering conditions, the strength of the artificial dam should not be less than 25 MPa.

Conclusions
A numerical hydro-mechanical creep model is proposed in this paper to investigate the stability and failure characteristics of brittle rocks in an underground water reservoir.A series of biaxial compression tests are carried out to validate the proposed hydro-mechanical creep model.The simulation results demonstrate that the proposed numerical model is capable of simulating rock failure.In order to study the influence of water pressures and dam strengths on the stability of artificial dam of underground reservoir, A simplified 3D numerical damage model is established based on the conditions of underground reservoir in Ningdong mining area.According to the engineering geological conditions of Ningdong mining area and the specific parameters of the artificial dam, the maximum water pressure that the artificial dam body can withstand is 25 m.The cutting depth has certain influence on reducing the displacement of artificial dam, and it should be no less than 0.6 m.The strength of artificial dam has great influence on the stability of artificial dam body, and the strength of the artificial dam should not be less than 25 MPa.The results of these numerical simulations indicate that the proposed numerical model is capable of investigating the stability and failure characteristics of rock in underground water reservoirs.

D
Damage variable E, E 0 Young's moduli of damaged material and undamaged material F 1 , F 2 Tensile and shear damage threshold functions r t0 , r c0 Uniaxial tensile strength and uniaxial compressive strength e t0 , e c0 The maximum tensile and compressive strain r 1 , r 2 The maximum and minimum principal stress f i Body forces G Shear modulus k Lame's constant k The coefficient of permeability K 0 Bulk modulus of the porous medium K s Effective bulk modulus e V Volumetric strain

Figure 1 .
Figure 1.Elastic damage constitutive law for an element under uniaxial compression and tension.

Figure 2 .
Figure 2. Flow chart of the coupled hydro-mechanical model.

Figure 4 .Figure 5 .
Figure 4.The simplified three-dimensional model of underground reservoir in Ningdong mining area and.

Figure 6 .
Figure 6.Damage distribution diagram and displacement of artificial dam on upstream surfaces and downstream surfaces at different head pressures.(a) head pressure 15 m; (b) head pressure 20 m; (c) head pressure 25 m; (d) head pressure 30 m; (e) head pressure 35 m; (f) head pressure 40 m.

Figure 7 .
Figure 7. Damage distribution of artificial dam on upstream and downstream surface with different side cutting depth.

Figure 8 .
Figure 8. Damage distribution and displacement of artificial dam on upstream surfaces and downstream surfaces with different roof and floor cutting depth.

Figure 9 .
Figure 9. Results of maximum principal stress under different dam strength.(a) Curve of maximum principal peak stress at different dam strength.(b) Maximum principal stress on upstream surface at the dam strength of 25 MPa.(c) Maximum principal stress on downstream surface at the dam strength of 25 MPa.

Figure 10 .
Figure 10.Damage distribution of artificial dam on upstream and downstream surface at different dam strength.

Table 1 .
Mechanical properties of numerical samples at the mesoscopic level.