Coal is the most abundant and widely distributed conventional energy source in the world, (Miller 2005; Thomas 2013). China’s proven coal reserves amount to 143.2 billion tons, representing approximately 13.3% of the global total proven coal reserves (Amoco 2021). In 2021, China’s coal production reached 80.91 exajoules, accounting for 50.7% of the world’s total coal production (Amoco 2022). In the future, coal is expected to remain a significant component of China’s primary energy mix due to a continued rise in demand (Zhang et al. 2017; Yang et al. 2020; Li 2021). However, the coal mining industry is confronted with numerous challenges, particularly regarding resource exploitation. Firstly, as coal resources are being continuously extracted, China’s shallow coal reserves are approaching depletion. At present, the coal mining industry in China is shifting its focus from shallow to deep mining (Xie et al. 2012), at a rate of 8–25 m deeper per year, with more than 50 coal mines being mined at depths over 1000 m (Kang et al. 2023). With increase in mining depth, the coal mining is facing higher stresses (Brown and Hoek 1978; Liu 2011), higher temperature (Torres and Singh 2008; Duque et al. 2018), higher water pressure (Alam et al. 2014), and stronger mining-induced seismicity and disturbances (Ghosh and Sivakumar 2018; Keneti and Sainsbury 2018). The coal in such an environment exhibit more severe failure and poses challenges for prediction of instability, with release of greater magnitudes of energy, potentially resulting in human casualties and economic losses (Chase et al. 1994; Galvin 2016).

In terms of geographical regions, it is expected that coal resources in eastern China will be depleted soon, and the focus of coal mining is shifting to the western regions (Zhang et al. 2012; Ling et al. 2022). While the western regions of China possess over 70% of the country’s proven coal reserves, it holds merely 3.9% of the total national precipitation (Li 2021). The disproportionately distributed coal and water resources in the western regions of China exhibit an extreme imbalance (Cheng et al. 2016; Chi et al. 2022). Coal mining causes surface subsidence as a result of a significant amount of mine water and related poro-elastic effects in the strata. On average, about 2 tons of mine water are produced for every ton of coal mined in China (Gu et al. 2021). Traditional mining methods often do not incorporate efficient water management practices, leading to a substantial loss of water resources. Usually, mine water is not adequately captured, treated, and reused. Instead, it is discharged as wastewater, causing both water losses and environment contamination. This is a disaster for the ecological system of the western regions in China, where coal and water are distributed disproportionately (Lv et al. 2019), seriously disrupting local water circulation. This unsustainable water management highlights the need for improved strategies in mining industry to maximize water reuse. Effective utilization of mine water is an important approach to protect the ecological stability of arid areas. Based on decades of engineering practices and theoretical research, Gu (2015) proposed a new framework for mine-water storage. The framework proposed is centered around the principles of water transfer, storage, and utilization. It entails channeling the mine water into an underground water reservoir known as Coal Mine Underground Water Reservoir (CMUWR). The CMUWR is constructed using artificial dams and coal pillars, which serve as storage structures for the water, see Fig. 1a. The aim is to effectively capture and store the mine water, allowing for its subsequent use in various applications and minimizing water wastage in the mining process. Currently, this technology has been used in more than 30 coal mines in China, with a maximum water storage capacity of 31 million cubic meters (Gu et al. 2021), achieving efficient recycling of mine water in the arid areas. With applications of CMUWR, the safety and stability of coal pillars has become a key issue. Previous studies have mainly focused on the stability evaluation of coal pillars during the mining process (Zhang et al. 2018; Sinha and Walton 2019; Wang et al. 2019b; Xia et al. 2021), while the analysis of the safety and stability of coal pillars during the operation of an underground water storage is not yet sufficiently revealed. Induced stresses in coal pillars in CMUWR are due to various sources, including static or dynamic stresses from overlying strata, rock mass collapse, water pressure, mining activities, and seismic events in adjacent areas, see Fig. 1b. Additionally, the coal pillars experience periodic immersion due to fluctuating water levels within the CMUWR. Understanding the behavior of coal pillars under these conditions is essential for ensuring their safety and stability and more research is needed to investigate the effects of these factors on the integrity and bearing capacity of coal pillars within a CMUWR system.

a An overview of the CUMWR b Possible directions and types of loads acted on coal pillars in CMUWR c The designed laboratory tests based on the actual load pattern of coal pillar dams

Currently, the design of underground coal pillars mainly revolves around optimizing structural dimensions (Reed et al. 2017; Kumar et al. 2019; Mark and Gauna 2021; Tuncay et al. 2021), as well as investigating the effects of water immersion (Poulsen et al. 2014; Wang et al. 2019a), wet-dry cycles (Yu et al. 2022), and loading patterns (Lu et al. 2020) on the strengths of coal. Experimental studies on the influence of water content confirm that water immersion significantly reduces the mechanical properties of some rocks and coal, such as the uniaxial compressive strength (UCS), tensile strength, and elastic modulus (Hawkins and McConnell 1992; Erguler and Ulusay 2009; Wild et al. 2015; Vishal et al. 2015; Kim and Changani 2016; Yao et al. 2016; Maruvanchery and Kim 2019). Additionally, the different distributions of water within the sample can result in differences in rock strength, an increase in water content can exacerbate crack propagation (Dusseault and Fordham 1993). Through compression-shear tests, it is found that the failure surface of saturated coal samples differs from that of dry coal samples, as the fracture propagation deviates from the shear plane (Yao et al. 2020; Chen et al. 2023).

The aforementioned studies indicate that current investigations on coal failure behavior stem mainly from lab experiment. The investigation of the failure mechanism at the micro-scale is still a subject of ongoing research. Further, research on the stability of coal pillars under complex stress paths is not yet systematically performed. Coal pillars, as the critical structures in CMUWR, are subjected not only to the vertical stresses caused by overlying strata and hydrostatic pressure, but also shear stresses resulting from mining-induced seismic events. Therefore, it is crucial to comprehensively consider the mechanical responses and failure modes of coal pillars under complex stress paths and immersion effects to reveal their underlying failure mechanisms. By taking into account the combined influences of these factors, a deeper understanding of how coal pillars behave and fail in such conditions can be gained. This article employes a particle-based discrete element method (DEM) to perform numerical simulations, providing a microscopic perspective of failure mechanisms of coal pillars under various water immersion and complex stress paths, which include the monotonic shear load (MSL) and cyclic shear load (CSL) under constant normal load (CNL) and dynamic-cyclic normal load (DNL). The primary objective is to reproduce major experimental results and gain insights into their failure mechanisms through numerical analysis. Furthermore, the study utilized a self-designed AE algorithm to analyze the microseismical characteristics of coal under dynamic stress paths.

2.1 Configuration of the numerical model

Shear tests were performed on coal samples that had undergone different levels of water soaking, while subjecting them to normal stress. The resulting experimental datasets were utilized to calibrate numerical simulation. The particle based simulations are often carried out to replicate the direct shear behaviors of rocks under constant load (Bewick et al. 2014a, b), the results show that the mineral grain and grain boundary strength both have dominant effects on rupture pattern. For cyclic shear tests, the loading regimes commonly employed in laboratory tests are either stress- or strain-controlled. In this study, a stress-controlled loading mode was adopted. Figure 2 shows the testing device used in lab test and the configuration of the numerical model. The blue particles (clumps) in Fig. 2b represent the shear box, which consist of an upper and lower loading jaw with a gap (approx. 0.625 mm) in the middle. The upper loading jaw can only move vertically, while the lower jaw can only move horizontally. The loading jaw consists of 658 overlapped pebbles with the same diameter, ensuring the flatness of the contact surface and preventing stress concentration caused by uneven contacts. The coal sample is represented by over 7000 particles. The contact between particles inside the coal sample model is Linear Parallel Bond Model (LPBM) (Potyondy and Cundall 2004) and the contact between the loading jaws and the coal sample is modelled via the Linear Model (LM). The detailed introductions for LPBM and LM are provided by Itasca Consulting Group Inc. (2014). When the local stress between particles exceeds the assigned bond strength, the connection between particles breaks, resulting in a linear contact, and strength parameters of the bond group in LBPM are longer available. Therefore, the local damage is characterized by bond failure between cemented particles.

A general overview of the experimental-set up and numerical model a shear box, the horizontal movement of upper jaw is fixed and normal load is transferred via the upper jaw, the lower jaw can move horizontally b numerical model configuration with boundary conditions (i.e., the upper jaw can only move vertically and the lower jaw can only move horizontally). c an overview of the DJZ 500 shear box. Details were given in (Dang et al. 2022)

2.2 Experimental datasets

An overview of the servo-controlled testing device (DJZ-500) is shown in Fig. 2c. The detailed technical parameters of DJZ-500 are as follows (Dang et al. 2022):

• Max normal force 500 kN, Max shear force ± 500 kN

• Frame stiffness 5 GN/m

• Dynamic stress/strain frequency up to 10 Hz, data acquisition frequency Max 10 kHz.

• Maximum shear displacement of ± 40 mm

• Shear velocity 0.01 to 240 mm/min

• Accuracy of load measurement < 0.1 kN,

• Accuracy of displacement measurement < 0.1 μm.

The coal specimens used in this study are cubic with (i.e., 100 mm × 100 mm × 100 mm). Six sets of stress–strain data were used for numerical model calibration, including two sets of monotonic shear load (MSL) tests and four sets of cyclic shear load (CSL) tests. Petrophysical and geometric properties of the six coal samples were measured prior to tests (see Table 1). Coal samples M1 to M4 were in the naturally dried state, while coal samples M7 to M8 underwent a treatment of 30-day water immersion. M1 and M2 were subjected to MSL with constant normal load (CNL) and dynamic-cyclic normal loads (DNL), respectively. The obtained shear stresses at failure were used for designing the stress levels in the subsequent cyclic shear tests. The detailed loading schemes are shown in Table 2. Dry samples M3 and M4 and wet samples M7 and M8 were subjected to CNL and CSL tests were performed under frequencies of 0.1 and 0.2 Hz (see Table 3). Figure 3 shows the stress paths and shear strengths of the six coal samples. Figures 3a,b show the paths of MSL (dry samples M1 and M2). Figure 3a shows CNL, while Fig. 3b shows DNL. Figures 3c–f show the CSL tests under CNL with different shear frequencies on dry and wet samples. Figure 3g suggests that different normal load boundary conditions (CNL or DNL) lead to differences in shear strength. Coal under MSL exhibit larger strength in CNL compared to DNL. The difference in shear strength under different frequencies is more significant. At a loading frequency of 0.2 Hz, the shear strengths reached 9.7 MPa (M4) and 8 MPa (M8). For a lower frequency (0.1 Hz), a shear strength of 6 MPa (M3) and 4.9 MPa (M7) were found. Furthermore, for both shear frequencies (0.1 and 0.2 Hz), the fatigue life of dry samples (0.1 Hz M3: 56 cycles; 0.2 Hz M4: 83 cycles) was consistently higher than for wet samples (0.1 Hz M7: 41 cycles; 0.2 Hz M8: 72 cycles).

Stress paths for the 6 coal samples and statistics of shear strengths: Monotonic shear test: a M1 b M2; Cyclic shear test: c M3 d M7 e M4 f M8 g Statistics of the shear strengths at failure

3.1 Assignment of microscopic parameters

Due to the differences in stiffness and strength between coal samples, each model representing an individual test requires calibration. The ratios of shear to normal stiffness, including the linear group k_s/k_n and the parallel bond group k_s′/k_n′, need to be determined from the normal and shear stress–strain curves. The ratios of shear to normal stiffness were different under MSL and CSL tests. In MSL, the ratio of shear to normal stiffness is 1.5, while the ratio is 2 in case of the CSL simulation. Subsequently to stiffness parameters calibration, a calibration was performed to determine the values of the strength parameters, including tensile strength σ_t, cohesion c and friction angle ϕ. The parameters are compiled in Table 4 and the process will not be detailed here.

3.2 Monotonic shear test

3.2.1 Verification of stress–strain curves

The shape of stress–strain curves were calibrated by adjusting the microscopic parameters to match the general shape and stiffness of the experimental stress–strain curves. First, the calibration of the MSL test was conducted (i.e., M1 and M2). Figure 4 illustrates the stress paths in the numerical simulation versus timestep, which are consistent with the stress paths obtained from experiments (Figs. 3a, b). Figure 5 compares the normal and shear stress–strain curves of selected experiments (red solid lines) and numerical simulation (black solid lines). In numerical simulations, when the samples are subjected to both normal and shear stresses, the fit of the numerical data with the measured stress–strain curves is moderate as compared to the case of being subjected to only normal stress. This may be attributed to the relatively soft nature of the coal. During the shear process, the upper shear box compacts the generated coal fragments, preventing the occurrence of significant normal dilation. In the PFC model, the sample is composed of rigid spheres that do not deform or break, presenting a significant disparity between numerical simulation and experiments. This feature causes a substantial normal dilation in the upper part of the specimen in numerical simulation, induced by rigid spheres protruding from the shear surface. Nevertheless, the simulated stiffness and strength shown in Fig. 5 sufficiently replicate the experimental results.

Numerical simulation of stress paths for two monotonic shear tests, including normal and shear stresses: a M1 b M2, the five indicated time intervals ①–⑤ are selected to present simulation results in the following sections

Normal and shear stress–strain curves for two MLS tests (M1 and M2), the M1 was subject to the CNL and M2 was subject to DNL load, the results include the lab test and numerical simulation a Normal stress–strain curve of M1 b Shear stress–strain curve of M1 c Normal stress–strain curve of M2 d Shear stress–strain curve of M2

3.2.2 AE energy and b-value

In the numerical model, the particles are bonded via LPBM. The occurrence of an AE event is characterized by the breakage of a bond between the cemented particles. When the bond breaks, the energy stored in the bond is released and detected and localized in real time. The Gutenberg-Richter relationship is widely used to characterize the intensity of rock damage based on AE monitoring (Iturrioz et al. 2014; Schütz and Konietzky 2016; Liu et al. 2020), which was originally used for earthquake characterization, and its mathematical expression is as follow (Gutenberg and Richter 1944):

where N is the number of earthquakes in the range of M to (M + ∆M), M is the magnitude of an earthquake, a and b are fitting parameters, where b is the parameter characterizing the relationship between the magnitude and the number of earthquakes. The AE event produced by rock fracturing is usually considered as a micro-earthquake, and thus the change in b-value can be used to characterize the intensity of rock damage and the energy magnitude. Figure 6 shows the spatial distribution of AE events for the two MSL tests (M1 and M2). Figure 6 illustrates the development of AE events during the experiment at five loading stages. The sizes of the circles are proportional to the magnitudes of released kinematic energy. Rows a and b in Fig. 6 show AE events at the selected loading time intervals, while rows c and d show all AE events. For both samples M1 and M2, initial microcracks emanate around the gap between the upper and lower jaw of the shear box. As the shear strain increases, the patterns of crack propagation exhibit different characteristics. In the case of DNL (M2), the cracks mainly develop along the horizontal shear plane, forming a macroscopic horizontal crack (Fig. 6d). In the case of CNL (M1), the coalescences of microcracks deviate from the horizontal plane and a cluster of wing cracks form. Furthermore, a significant difference in the magnitudes of AE energy is observed. For M1 the magnitude of AE energy is larger, indicating more intensive damage processes. Figure 7 shows the evolution of the total AE events for the two MSL tests (M1 and M2). In region B of M1, the green areas correspond to the testing phase in which the number of AE event increases, indicating the development of wing cracks on both sides of the sample. In region C of M1, the green area indicates a sharp rise in the number of cracks prior to shear failure. Wing cracks continued to propagate, eventually forming a macroscopic shear plane by coalescence of an individual crack, crossing the entire specimen laterally. In the case of M2, the number of AE events in region B show several times a sudden rise. In region C, the increase in AE events occurs exclusively at low normal stresses, indicated by the green areas. With increasing normal stresses during cyclic loading, the number of AE events decayed, suggesting a normal stress threshold that controls failure processes and thus AE activity. During unloading, AE events are detected again when the normal stress approaches the above-mentioned threshold stress. This behavior is likely attributed to the fact that for low normal stress, the shear stress exceeds the shear resistance within the coal matrices. As the normal stress increases, the applied shear stress becomes insufficient to overcome the enhanced frictional resistance (Dang et al. 2019, 2021).

Spatial distribution of AE events for two MSL coal samples, M1 (CNL) and M2 (DNL), ①–⑤ indicate 5 time intervals shown in Fig. 4. a Newly generated AE events for M1 b Newly generated AE events for M2 c Total AE events for M1 d Total AE events for M2, the sizes of the circles are proportional to the dimensionless magnitudes of AE energy (see the legend at the top)

Evolution of normal stress and cumulative AE events vs. timesteps for two MSL tests a Overall evolution of AE for M1 b Zoomed-in region B for M1 c Zoomed-in region C for M1 d Overall evolution of AE for M2 e Zoomed-in region B for M2 f Zoomed-in region C for M2

Figure 8 illustrates the distribution of numerical b-value for the two MSL tests at the five time intervals documented in Fig. 4. All magnitudes of the modelled AE energy were scaled by a constant factor in order to achieve realistic magnitudes of micro-seismic events. The sample sizes for calculating the b-value were 70 AE events. Figure 8a shows the evolution of b-values obtained from the M1 specimen, indicating a significant decrease (from ④ to ⑤) in b-value preceding the failure. Figure 8b show a relatively moderate variation in b-value for M2, with a decrease prior to failure (from ③ to ⑤).

Distribution of b-values of two MSL tests at five time intervals documented in Fig. 4a b-value at five time intervals for M1 b b-values at time intervals for M2

3.2.3 Stress distribution during monotonic shear

The localized stress tensor between particles was calculated by volume average of summing up contact reactions according to Eqs. (2)–(5) (Peters et al. 2005; Zhang and Konietzky 2022). Where V is the volume of particles, N^c is the number of contacts, f_i^c the ith component of the force acting at the contact, and r_j^c the jth component of the radius vector from the centre of particle to the point of contact, σ_ij are the components of localized stress tensor. σ₁/σ₃ are respectively the maximum and minimum principal stresses. Here we adopt a tension-positive and compression-negative convention, σ₃ is hence the most compressive principal stress. The angle of the σ₃ with respect to vertical (positive y-axis) was calculated using Eq. (5). The calculation of the maximum shear stress (τ_max) is shown in Eq. (6). The distribution of the σ₃ for M1 and M2 are shown in Figs. 9a, b. Due to significant shear displacements at failure, only stress distributions within an area covering 80% of the central specimen surface (dimension of 8 cm × 8 cm) are shown. The five figures correspond to the five time intervals are documented in Fig. 4.

Stress distributions during two MSL tests, the results are based on the five time intervals documented in Fig. 4a σ₃ distribution of M1 b σ₃ distribution of M2 c τ_max distribution of M1 d τ_max distribution of M2. The stresses are the corresponding components of localized stress tensor calculated by volume average of summing up particle contacts, the magnitudes of the stress can refer to the legend at top

The stress distribution in Fig. 9a shows that under CNL, when the applied shear load is low, the concentration of σ₃ initially occurs around the gap between the upper and lower jaws of the shear box. As the shear load is gradually increased, the regions of stress concentration grow at an angle of 60°–80° with respect to the vertical. When failure occurs, the angle approaches approximately 45°, and low-stress regions are observed in the upper left and lower right corners. For DNL conditions (Fig. 9b), the concentration zone of σ₃ appears at an angle that is normal stress dependent, and approach angles significantly less than 50° with respect to vertical prior to failure. The patterns of τ_max distributions for M1 (Fig. 9c) and M2 (Fig. 9d) are respectively similar to the distributions of their σ₃ in location and shape. However, the width of the zone is larger under CNL compared to DNL. This region coincides with the formation of a series of wing cracks.

Figure 10 shows the relations of magnitudes and angles between σ₃ and τ_max. The ratio of |σ₃| to τ_max is shown as colored dots and the values can be determined using the legend at the top. The right y-axis shows the angle (θ) of σ₃ in respect to vertical calculated from Eq. (5). M1 and M2 exhibit the following common patterns (see Fig. 10):

Evolution of relations between the ratio of |σ₃| to τ_max for two MSL tests, the results are based on the five time intervals in Fig. 4: a Evolution of M1 under CNL b Evolution of M2 under DNL, the right axis indicates the angle of σ₃ with to vertical and corresponds to the values of x shaped black symbols, the colorful dots correspond to the bottom and left axis where the values for σ₃ and τ_max can be read, the ratio of |σ₃| to τ_max is determined using the legend at the top c Explanation of the ratio of |σ₃| to τ_max and its implications based on M-C criteria implemented in LPBM model

1. (1)

   The σ₃ concentration zone tends to develop towards ± 45°.

2. (2)

   For low |σ₃|, the ratio of |σ₃| to τ_max, spans a broad range, and its maximum value exceeds 8. As |σ₃| increases, the ratio drops toward 2 and displays a significant convergence.

In LPBM, the failure criterion for the bond between two particles is controlled by assuming a Mohr–Coulomb (M-C) envelope, where τ_max is equal to 1/2 |σ₁–σ₃|. Therefore, the ratio of |σ₃| to τ_max is 2, when σ₁ = 0. This indicates no tensile state for the maximum principal stress and the system is in a critical state. When the ratio of |σ₃| to τ_max is lower than 2, σ₁ is tensile with the Mohr circle moves left. This implies that the probability of failure, either by tensile or shear fracture significantly increases. For higher ratio of |σ₃| to τ_max (e.g., 8), the probability for failure is lower. Hence, a decrease in ratio of |σ₃| to τ_max increases the probability of bond failure. When this ratio is less than 2, failure can occur in both tensile and shear modes. It also explains why as stress increases, the ratio decreases and approaches 2.

3.3 Cyclic shear tests

3.3.1 Stress–strain curves and bond diameter multiplier

The micro parameters assigned to the numerical assembly are influenced by the decay in bond diameter, which originates from the parallel-bonded stress corrosion (PSC) model proposed by (Potyondy 2007) to simulate the rocks under static fatigue with constant creep load. The PSC model represents the rock as a bonded granular material and then removes the bond based on different rates under the external load. The reduction in bond diameter is facilitated by changing the bond radius multiplier as shown below:

where D, R, λ are respectively the bond diameter, radius and diameter multiplier, the values before and after reduction in bond diameter are denoted by non-superscript and superscript symbols. If the reducing rate of bond diameter is constant, it can be expressed as Eq. (8), ∆t is the duration and γ is a fixed reducing rate. The interrelations between the bond force, displacement, rotation, and effective stiffness are shown in matrix form in Eq. (9), where A, I, and J are the area, moment of inertia, and polar moment of inertia of the bond cross-section. F^pn and F^ps are normal and shear components of the parallel bond force. M_t and M_b are twisting and bending rotation moments. U_n and U_s are local normal and shear displacements of two particles. θ_t and θ_b are local twisting and bending rotation between two particles. According to (Potyondy 2007), when the diameter of a bond is reduced, the effective stiffness is weakened because of the reduction in A, I and J, whereas the stiffness per unit area, k_n′ and k_s′ are not affected.

The updated stiffness matrix after damage induced from bond diameter reduction is shown in Eq. (10). The local translational and rotational strains and rotation moments after reduction of bond diameter will be enlarged when exposed to a fixed loading rate, see Eq. (11).

In our previous research on rock fatigue using modified PSC models, the nonlinear evolution path of bond multiplier can properly simulate the stress–strain curves under both, constant or multi-level cyclic loading. The path is directly linked to the magnitude of stress increments, number of cycles, and sudden strain increases. The four samples (M3, M4, M7, M8) subjected to CSL tests exhibit distinct stress–strain curves prior to failure, hence requiring calibration of each dataset. Figure 11 shows the applied evolution paths for bond diameter multipliers and the stress paths over timesteps for the four samples. Figure 12 summarizes the evolution paths of bond diameter multipliers. All the multipliers are set to 1 to indicate a pristine state before loading. For M4 and M8, with a loading frequency of 0.2 Hz, the decrease in bond diameter multipliers is relatively gentle, remaining above 0.8 at failure. However, for M3 and M7, with a frequency of 0.1 Hz, exhibits a more intensive drop in bond diameter multiplier by reaching a much lower value of 0.12 at failure. According to the testing results, the frequency has a notable impact on strengths of coal sample under CSL. A lower frequency can result in a shorter fatigue life and premature failure, therefore the evolution for the bond diameter multiplier under the lower frequency (M3 and M7) display a sharper decline. The results also suggest (1) for a given frequency the strength remains almost equal for wet and dry coal (2) frequency is dominating the damage evolution, see Fig. 12.

Evolution of bond diameter multipliers and stress paths versus timesteps for four selected samples a M3 0.1 Hz b M4 0.1 Hz c M7 0.2 Hz d M8 0.2 Hz. the five moments ①–⑤ indicated are the selected time intervals to show the simulation results

Comparisons of the evolution paths of the bond diameter multipliers versus the cyclic loading stage (CLS) for the four selected samples, the data for the 5th CLS of M7 is not presented

Figure 13 show normal (Figs. 13a–d) and shear (Figs. 13e–h) stress–strain curves respectively obtained from lab tests (yellow dashed lines) and numerical simulations (black solid lines). Experimental results show nonlinear, concave upward stress/strain curve at low stress which corresponds to closure of pores and voids. This phase was not considered for numerical back-calculations Figs. 13a–d show that the elastic moduli obtained from simulation are consistent with those obtained from lab testing. Since the non-linear part of the experimental stress strain curve was not considered in the calibration, numerical results were shifted in Fig. 13 along the X-axis to fit the linear elastic portion of the experimental stress strain curve (referred to as strain shift). At peak strength, the normal strains of Figs. 13b, d suggest that failure is accompanied by dilation associated with abrupt change from a bulk compression to extension.

Comparison of normal and shear stress–strain curves between lab tests and numerical simulations based on 4 samples under CNL and CSL, where the black curves are for simulations, the gray arrows indicate the strain shift, and the yellow curves are for laboratory tests: Normal stress–strain curves: a Dry sample M3 b Dry sample M4 c Wet sample M7 d Wet sample M8; Shear stress–strain curves: e Dry sample M3 f Dry sample M4 g Wet sample M7 h Wet sample M8

The shear stress–strain simulation results are shown in Figs. 13e–h. In the first two cyclic loading stages, the shear hysteresis loops indicate a high stiffness, and little to no plastic deformation. For samples M3 and M7, oscillations in the stress strain curve at low shear stresses are almost vertical (i.e., the shear strain is almost zero). This is likely due to the fact that the absolute magnitude of the shear stress is low compared to the applied normal load. For two samples it was even observed at low shear stress that the lateral strain induced by the applied normal load exceeds the shear strain induced by the applied shear stress, see Figs. 13e, g. Thus, the above described behavior of the samples at low shear stresses is not used for numerical model calibration. Figure 13 further suggests from hysteresis loops that damage accumulates with increasing number of cycles. However, for one sample M3, the damage increases substantially during load cycles close to failure.

3.3.2 AE characteristics and crack patterns

Figure 14 presents the spatial distribution of AE events for four selected samples under CNL-CSL. Rows a-d show the generated AE events for five selected time intervals (Fig. 11), while rows e–h show the total AE events. Similar to the MSL tests documented in Fig. 6, initial cracks occur at the gap between the upper and lower jaws of shear box. As the cyclic shear test is continued, cracks continue to propagate, and coalesce toward the center of the sample, eventually penetrating the whole sample. M3 and M4 ultimately exhibit a horizontal shear zone, accompanied by a limited number of wing cracks. The amplitudes of AE energy in M3 were generally lower, while the corresponding AE energy for M4 was significantly higher. This is likely due to the lower shear stress at failure and thus the lower quantity of stored energy. Cracks in M8 continuously propagate along horizontal plane, ultimately forming a macroscopic horizontal crack without a wing crack. For M7, initial cracks initially form mainly at the contact between the samples and the shear box. With increasing shear stress, new cracks appear in the center of the sample eventually forming a sub-horizontal shear band that further increases in the width with increasing shear stress. Figure 13c and Table 4 show that M7 has a much lower normal elastic modulus and parallel bond stiffness compared to the other samples. Thus, during MSL, interparticle dislocations are more pronounced resulting in more local stress concentrations and broader distribution of microcracks. This may explain why in M7 a shear zone rather than a clear shear plane is ultimately forming.

The spatial distribution of AE events for the model of four selected samples under cyclic shear loading, ①–⑤ are the selected five time intervals in Fig. 11, distribution of newly generated AE events: a M3 b M4 c M7 d M8. Distribution of total generated AE events: e M3 f M4 g M7 h M8

Figure 15 shows the distribution of b-values for the four selected samples. At time interval ①, when the sample starts to fail, the b-value is relatively high, indicating that AE activities generated at this time interval are mainly small-energy events. It is observed that the b-values of the four samples drop prior to shear failure (from ④ to ⑤), indicating an increase in large-magnitude AE events. This significant decrease in the b-value can be considered as precursor signal for coal failure. Compared to other samples (M4, M7, and M8), the b-value of the coal sample in Fig. 15a (M3) is obviously higher at all five time intervals, indicating a predominance of small-energy AE events. This is consistent with the lower magnitudes of AE energy in Fig. 14 for M3. The normalized seismic magnitude of M3 in Fig. 15a is around 4, while the normalized magnitudes of the other three samples are larger than 10, with M4 in Fig. 15b reaching a maximum magnitude of around 21.

Distribution of numerical b-values of four selected samples at five time intervals a The distributions of b-value at time intervals for M3 b The distributions of b-value at five time intervals for M4 c The distributions of b-value at five time intervals for M7 d The distributions of b-value at five time intervals for M8

We monitored the fracturing processes of two coal samples under CNL-CSL by recording videos (M9 and M10, which have been exposed to the same stress paths as M3, M7 and M4, M8). The traces of the cracks are shown in Fig. 16. Their macroscopic failure surfaces exhibit approximately horizontal planes without obvious wing cracks. At failure, some branching fractures formed corresponding to the distribution of located AE events shown in Fig. 14. The crack pattern observed under CNL-CSL differs significantly from the pattern shown in Fig. 6 for sample M1 under CNL-MSL. The results suggest that for CNL-CSL conditions, cracks occur predominantly along a horizontal plane, while for CNL-MSL conditions a complex pattern of horizontal shear planes and wing cracks form. Our simulation results are consistent with experimental observations.

Sketches for the evolution of cracks of two samples M9 and M10, the red lines are the traces of macrocracks

3.3.3 Spatial distribution of σ ₃ and τ _max

Figure 17 presents the spatial distribution of σ₃ at the five time intervals for four selected CNL-CSL tests. Stress distributions is shown within an 8 cm × 8 cm area covering 80% of the specimen surface. σ₃ concentrations occur at the middle of both shear box side walls. With the continued application of cyclic shear, the numerical values of σ₃ progressively increase near the center of the specimen, forming a distinct high-stress region. This high-stress region gradually expands towards the center, converges, and eventually interconnects. The band of stress concentrations in the central part of the sample is inclined with an angle of 65°–75° from vertical. The magnitude of the inclination is similar to M2 (60°–80°) reported in Fig. 9b and is much larger than that for M1 (< 50°) in Fig. 9a, which is subject to CNL-MSL. This indicates that both DNL and CSL can significantly increase the angle between the stress concentration regions and the vertical direction at failure. This may also explain why the failure surfaces under CSL conditions are rarely associated with wing cracks formation. In Fig. 17, the fifth image in each row represents the stress distribution in the post-peak phase. The shear displacement increases significantly in this phase and the previous distribution of σ₃ is disrupted. Specifically, the high-stress concentrations no longer penetrate through the central part of the specimen. This is caused by loose contacts between particles on the shear plane. The dislocations of particles interrupt the stress transfer and explain the disappearance of high-stress regions during the post-peak phase. However, Fig. 17a, also indicates two remaining stress concentrations, likely due to interlocked particles. The distribution of τ_max is presented in Fig. 18. The general pattern is similar to the pattern of σ₃. However, the magnitude for τ_max is lower than σ₃. Figure 19 shows the locations of AE events for three selected samples (M1, M4, M8) together with the distribution of σ₃ and τ_max. The occurrence of AE events with large AE energy are typically found within high-stress regions. The distribution of stress concentrations is influenced by factors such as particle grading, stress path and fabrics.

σ₃ distribution for four selected coal samples under CNL-CSL: a σ₃ distribution for M3 b σ₃ distribution for M4 c σ₃ distribution for M7 d σ₃ distribution for M8, the magnitude of stress can refer to the legend at top

τ_max distribution for four selected coal samples under CNL-CSL: a τ_max distribution for M3 b τ_max distribution for M4 c τ_max distribution for M7 d τ_max distribution for M8, the magnitude of τ_max can refer to the legend at top

Comparisons of σ₃, τ_max and AE events distribution a M1 σ₃ b M4 σ₃ c M8 σ₃ d M1 τ_max e M4 τ_max f M8 τ_max

Figure 20 illustrates the relations between σ₃ and τ_max for the four selected CNL-CSL samples. Similar to the distribution patterns of the two monotonically sheared samples presented in Fig. 10, it is observed that for low |σ₃|, the ratio of |σ₃| to τ_max covers a range between 0 and 8, indicating a stable state (Fig. 10c). As |σ₃| increases, the ratio approaches 2, which implies a critical state ₁. Further, the angle between the σ₃ with respect to vertical approaches ± 45° as |σ₃| increases. In the post peak phase, the angles are within a range of 0°–90° in general, and approximately 0° for M8.

Evolution of relations between the ratio of |σ₃| to τ_max for four CNL-CSL tests, the results are based on the five time intervals in Fig. 11: a evolution of M3 b evolution of M4 c evolution of M7, d evolution of M8, the right axis indicates the angle of σ₃ with to vertical and corresponds to the values of x shaped black symbols, the colorful dots correspond to the bottom and left axis where the values for σ₃ and τ_max can be read, the ratio of |σ₃| to τ_max can refer to the legend at top

As a sample approaches failure, decreasing b-values are often observed (Rivière et al. 2018; Dublanchet 2020; Ritz et al. 2022), suggesting that as the stress or strain increase, larger magnitude AE events become relatively more frequent. This change in the b-value is often associated with increased stress concentrations, crack propagation, and formation of localized failure zones, but depends on the material properties, loading regime and stress patterns. b-value evolution of six selected samples together with the number of AE events and the cumulative AE energy is shown in Fig. 21. Figure 21a shows the variation of b-value with the increase of AE events, and Fig. 21b shows the variation of the b-value with the cumulative AE energy. Specifically, the distributions of b-value can be categorized into two types: (1) An almost monotonic decrease (M4, M8); (2) An increase followed by a decrease shortly preceding failure (M1, M2, M3, M7). Numerical simulations suggest for both, monotonic or cyclic shear cases, a decrease in b-value prior to sample failure. This is consistent with the pattern of AE in rock mechanical experiments, indicating that the energy released through bond failure can be properly used to characterize magnitudes and to calculate the numerical b-value. Whether the b-value exhibits a monotonic decay is directly related to the evolution of bond diameter multiplier. For M4 and M8, the decay of bond diameter is much milder compared to M3 and M7, see Fig. 12, resulting in a monotonic decrease in b-value. This indicates a stable crack propagation. On the other hand, for sharp decrease in bond diameter (M3 and M7), local stiffness deterioration occurs, leading to disturbances in localized AE energy, which is reflected in the non-monotonic decay of b-value. By combining the above patterns, we can predict the strength of the sample using the evolution of b-value in numerical simulations, which is consistent with the experimental b-value evolution in rock tests.

Evolution of b-value of six selected coal samples with respect to the number of AE events and cumulative AE energy magnitude a Variation of b-value with the number of AE events (cracks) b Variation of b-value with AE energy. The solid lines are used to indicate data tendency

Due to the complexity of loads acting on CMUWR, the stability of coal dams is influenced by multiple factors including normal and shear load directions and water immersion. Our research numerically investigated the mechanical properties of coal samples under dynamic stress paths and dry–wet conditions. The results indicate that the patterns of normal and shear load significantly influence the distribution of minimum principal stress, maximum shear stress and the morphological characteristics of shear bands. This suggests that specific loading conditions, such as dynamic loads during rock strata fracturing, constant normal loads from strata weight, or dynamic shear loads caused by seismic waves, will result in different degrees of damage to the CMUWR. Based on the numerical simulation conducted in this study the following conclusion can be drawn:

1. (1)

   By employing a particle based numerical model, we replicate the stress–strain characteristics of coal samples subjected to complex dynamic stress paths in both, the normal and shear loading directions. Distinct damage evolution paths were established in numerical models for dry and water-immersed coal samples under two frequencies (0.2 and 0.1 Hz). Low frequency causes more server damage and thus a shorter fatigue life. Under the identical frequency, water-immersion results in a shorter fatigue life.

2. (2)

   Our algorithm successfully captures the magnitude of released energy at bond breakage, allowing for the determination of b-value and visualization of AE events over time including their energy magnitudes. The observed decrease in b-value prior to failure may have predictive significance for coal pillars. The simulation reveals the distribution of σ₃ and τ_max under monotonic and cyclic shear conditions. It is found that cyclic stress can reduce the probability of wing crack formation and also increase the angle between the σ₃ in respect to vertical at failure.

3. (3)

   When the ratio of |σ₃| to τ_max is high, the probability of failure is lower. When the ratio is less than 2, failure can occur in both tensile and shear modes. The significant number of bond failures is associated with the left shift of Mohr circle.

4. (4)

   Future research may consider the mechanical responses of coal under normal and shear stresses with different stress rates. It is recommended to introduce the random stress spectrum or cyclic loading with different loading and unloading rates. This will provide a more realistic characterization of in-situ loading conditions and stress patterns.