*Article* **Infrared Precursor Experiment to Predict Water Inrushes in Underground Spaces Using a Multiparameter Normalization**

**Kewang Cao 1,2, Furong Dong 1,2, Liqiang Ma 2,\*, Naseer Muhammad Khan 3,\*, Tariq Feroze <sup>3</sup> , Saad S. Alarifi <sup>4</sup> , Sajjad Hussain <sup>5</sup> and Muhammad Ali <sup>6</sup>**


**Abstract:** Rock failure is the root cause of geological disasters such as slope failure, civil tunnel collapse, and water inrush in roadways and mines. Accurate and effective monitoring of the loaded rock failure process can provide reliable precursor information for water inrushes in underground engineering structures such as in mines, civil tunnels, and subways. The water inrush may affect the safe and efficient execution of these engineering structures. Therefore, it is essential to predict the water inrush effectively. In this paper, the water inrush process of the roadway was simulated by laboratory experiments. The multiparameters such as strain energy field and infrared radiation temperature field were normalized based on the normalization algorithm of linear function transformation. On the basis of analyzing the variation characteristics of the original parameters, the evolution characteristics after the parameters normalization algorithm were studied, and the precursor of roadway water inrush was predicted comprehensively. The results show that the dissipation energy ratio, the infrared radiation variation coefficient (IRVC), the average infrared radiation temperature (AIRT), and the variance of successful minor infrared image temperature (VSMIT) are all suitable for the prediction of roadway water inrushes in the developing face of an excavation. The intermediate mutation of the IRVC can be used as an early precursor of roadway water inrush in the face of an excavation that is being developed. The inflection of the dissipation energy ratio from a declining amount to a level value and the mutation of VSMIT during rock failure can be used as the middle precursor of roadway water inrush. The mutation of AIRT and VSMIT after rock failure can be used as the precursor of roadway imminent water inrush. Combining with the early precursor and middle precursor of roadway water inrush, the graded warning of "early precursor–middle precursor–final precursor" of roadway water inrush can be obtained. The research results provide a theoretical basis for water inrush monitoring and early warning in the sustainable development of mine, tunnel, shaft, and foundation pit excavations.

**Keywords:** water inrush; precursor; strain energy; infrared radiation; normalized; sustainable mine

### **1. Introduction**

Due to the rapid development of underground spaces such as subways, tunnels, and caverns, large-scale geological disasters such as water inrushes, mud rushes, and rock collapses often occur in the construction stage [1–7]. In coal mining, water inrush accidents often occur, resulting in serious and often irreparable property losses and casualties, seriously affecting the normal production of coal mines [8,9]. It is therefore important to carry out relevant underground development face water inrush laboratory experiments to

**Citation:** Cao, K.; Dong, F.; Ma, L.; Khan, N.M.; Feroze, T.; S. Alarifi, S.; Hussain, S.; Ali, M. Infrared Precursor Experiment to Predict Water Inrushes in Underground Spaces Using a Multiparameter Normalization. *Sustainability* **2023**, *15*, 7570. https://doi.org/10.3390/ su15097570

Academic Editor: Jianjun Ma

Received: 16 February 2023 Revised: 26 April 2023 Accepted: 1 May 2023 Published: 5 May 2023

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

determine the variation characteristics of water inrush precursor information; this can not only improve the reliability of water inrush prediction in underground space development but also meet the needs of sustainable protection of mining water resources. Water resource protection is the theme of sustainable development of mining ecological environments and a key way to achieve green mining.

Mining and underground civil activities inevitably lead to the redistribution of underground rock stress and fracture damage of rock, which greatly change the permeability of surrounding rock. This can lead to water inrush mainly in the roof, through stress-induced natural fractures or geology features, causing safety accidents.

In subway construction, damage to the lining of shield tunnels may cause water inrush in the subway. This is caused by nonwatertight drilling of through-holes near the intersection of two subway tunnels. Under high hydraulic gradients, weak permeable areas form and extend, leading to soil water slurry explosions [10–13]. The infrared radiation information of rock changes during the process of stress redistribution and consequent fracture damage process [14–20]. Monitoring the infrared radiation released to the outside during the process of rock loading can predict the characteristics and process of rock deformation and failure. This provides reliable information for the establishment of rock failure precursors [21–29]. In recent years, many scholars have carried out considerable research using the infrared radiation characteristics of rock fracture and water seepage. Asakura et al. [30] studied the infrared radiation monitoring of water leakage in tunnel lining and proved the feasibility of monitoring water leakage by infrared technology. Liu et al. [31] studied the characteristics of infrared radiation in the process of concrete fracture and water seepage and found that the "initial increase followed by a decrease" in the curve of infrared radiation temperature was an abnormal precursor of the infrared radiation of concrete fracture and water seepage. Dou et al. [32] carried out the infrared radiation observation experiments of tunnel leakage, studied the infrared radiation variation characteristics, and wrote the MATLAB image processing program to extract the infrared image characteristics during the process of concrete leakage. Zhang [33] studied the characteristics of infrared radiation in the process of sandstone fracture and water inrush, in which the sudden decrease of infrared radiation temperature predicted the occurrence of roadway water inrush, and the accelerated rate of infrared radiation temperature could be used as a precursor for water inrushes.

To quantitatively analyze the characteristics of infrared radiation in the process of rock fracture and water seepage, Wu et al. [34] first proposed the average infrared radiation temperature (AIRT) index of rock surface. Liu [35] used the AIRT index to analyze the infrared radiation characteristics of dry and water saturated rocks during uniaxial loading, and found that water can promote the AIRT of a rock surface. However, different areas of the rock surface may simultaneously heat up and cool down during the unstable development crack development stage in the rock, resulting in no change in overall AITR index [36]. Therefore, Liu [37], Ma et al. [38], and Yang [39] proposed the infrared radiation variance (*IRV*), variance of successful minor infrared image temperature (*VSMIT*), and infrared radiation variation coefficient (*IRCV*) of rock surface. The results show that the above indexes can well reflect the differentiation characteristics of infrared radiation on the rock surface. For example, Ma and Zhang [38] studied the internal relationship between stress adjustment (due to excavating) and *VSMIT* index during the loading process of dry and water-saturated rocks, and found that *VSMIT* index can correlate well with the rock failure and water has an amplification effect on the mutation characteristics of *VSMIT* index.

Although scholars have carried out a lot of research work on the monitoring and warning of water inrush during underground development, the results to date have not proved reliable for detecting an imminent water inrush. This is mainly due to the localized complex geology and hydrology. Based on the analysis of the physical parameters such as stress, infrared radiation, and strain energy in the process of predicting roadway water inrushes, this paper normalizes each physical parameter, then comprehensively analyzes the multielement information evolution characteristics and their correlation with roadway water inrushes, and studies the comprehensive precursor characteristics of them. This study proposes an early warning precursor and monitoring the occurrence of roadway water inrush. The research findings will provide a theoretical basis for monitoring and early warning of water inrush in underground spaces for their safe and efficient development.

#### **2. Experimental Principle**

Stefan–Boltzmann's law states that any object above absolute zero will radiate electromagnetic waves to the outside world. Due to this, the radiation intensity and the temperature of the object satisfy the following formula [39]:

$$
\sigma^\* = \varepsilon \sigma T^4 \tag{1}
$$

where *J* ∗ represents the total energy radiated by the object per unit area, W·m−<sup>2</sup> ; ε represents the surface emissivity of the object, 0 < ε < 1; σ is the Stefan–Boltzmann constant, 5.670373 <sup>×</sup> <sup>10</sup>−<sup>8</sup> <sup>W</sup>·m−<sup>2</sup> ·K−<sup>4</sup> ; *T* is the absolute temperature of the surface, K.

*J*

This law explains how the infrared radiation on the surface of the object at room temperature is affected. Further, the relationship between the radiation intensity and the temperature of the object satisfies the fourth power. The force exerted on a solid causes changes in the distance between internal particles, resulting in thermodynamic changes and, thus, temperature changes. This phenomenon of temperature change caused by heat generated due to force can be referred to as the thermal–mechanical coupling effect [36]. Materials with different mechanical properties (such as elastic materials, elastoplastic materials, viscoelastic materials, etc.) and the same material have different thermal and mechanical coupling effects at different stress stages. These different thermodynamic coupling effects have different characteristics and laws due to their different microscopic mechanisms. In the elastic range, the object undergoes the process of tension or compression accompanied by the reversible conversion of heat. In the adiabatic environment, the sum of temperature and principal stress satisfies the linear relationship [40]:

$$
\Delta T = \frac{\alpha}{\rho \mathcal{C}\_{\sigma}} TS \tag{2}
$$

where ∆*T* and *T* are the change of object temperature and object temperature, respectively; α is the linear coefficient of thermal expansion; *C*<sup>σ</sup> is the specific heat coefficient under constant stress; *S* is the sum of the principal stresses.

The change of principal stress during uniaxial rock loading is only related to σ1. The fourth term represents the internal dissipated energy of the material, which is manifested in the thermoelastoplastic comprehensive effect at this stage. The work carried out by the external force is not all transformed into the internal thermal energy of the material but is mostly consumed in the process of internal microstructure change. Plastic deformation in the process of energy dissipation and thermal energy conversion is not reversible, and ∆*E* mainly includes the following three parts in the process of energy consumption [36]:

$$
\Delta E = \Delta E\_1 + \Delta E\_2 + \Delta E\_3 \tag{3}
$$

∆*E*<sup>1</sup> is the energy carried by the escape process of pore gas, in general, ∆*E*<sup>1</sup> < 0. ∆*E*<sup>2</sup> is the energy consumed by the expansion of pores, fractures, and joints and the generation of new fractures in the rock. The internal pores, fractures, and joint weak surfaces will first contract and close with the increase of stress. The pores will then collapse, and the primary fractures and joints will further expand, penetrate, and merge. Moreover, the new fractures will be generated by the increase in stress. The ∆*E*<sup>2</sup> is less than 0 due to the consumption of energy in this process. ∆*E*<sup>3</sup> is the energy generated by friction because there are friction behaviors among the pores, fissures, joints, and rock particles along all directions in the interior of the rock. Two factors influence the process of frictional heat generation: one is the positive pressure on the contact surface inside the rock, and the other is the friction coefficient. When the friction coefficient is fixed, the friction force is positively related to *3.2. Rock Samples* 

the normal stress on the contact surface. The larger the friction force, the more work will be carried out to overcome the friction force in the process of crack and particle sliding, resulting in the higher energy consumption. It is important to note that in this process the temperature of the contact surface due to friction heat production will increase, therefore ∆*E*<sup>3</sup> is greater than zero. **3. Experimental Design**  *3.1. Experimental Equipment*  The experimental loading equipment used an SANS electronic universal testing machine

#### **3. Experimental Design** system with a maximum vertical load of 1000 kN. The water pressure loading equipment

will increase, therefore Δ*E*3 is greater than zero.

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 4 of 18

fractures and joints will further expand, penetrate, and merge. Moreover, the new fractures will be generated by the increase in stress. The Δ*E*2 is less than 0 due to the consumption of energy in this process. Δ*E*3 is the energy generated by friction because there are friction behaviors among the pores, fissures, joints, and rock particles along all directions in the interior of the rock. Two factors influence the process of frictional heat generation: one is the positive pressure on the contact surface inside the rock, and the other is the friction coefficient. When the friction coefficient is fixed, the friction force is positively related to the normal stress on the contact surface. The larger the friction force, the more work will be carried out to overcome the friction force in the process of crack

#### *3.1. Experimental Equipment* adopts a Shanghai SB water pressure pump, the maximum working pressure is 10 MPa, and

The experimental loading equipment used an SANS electronic universal testing machine system with a maximum vertical load of 1000 kN. The water pressure loading equipment adopts a Shanghai SB water pressure pump, the maximum working pressure is 10 MPa, and the water pressure is set to 0.5 MPa. The infrared radiation detection device adopts the American FLIRA615 infrared thermal imager, whose thermal sensitivity is 0.025 ◦C, and the wavelength range is 7.5–14.0 m. The image acquisition rate was set at 25 frames/s. the water pressure is set to 0.5 MPa. The infrared radiation detection device adopts the American FLIRA615 infrared thermal imager, whose thermal sensitivity is 0.025 °C, and the wavelength range is 7.5–14.0 m. The image acquisition rate was set at 25 frames/s. The representative samples of sandstone collected from a coal mine in Shandong

#### *3.2. Rock Samples* province were used in the laboratory experimental process. All samples were obtained

The representative samples of sandstone collected from a coal mine in Shandong province were used in the laboratory experimental process. All samples were obtained from the same rock sample. The specimen design specification is a cuboid of 100 × 100 × 150 mm. The diameter and depth of the observation hole in the test block are 50 mm, and the diameter and depth of the water injection hole are 50 mm and 50 mm. A total of five specimens were prepared, represented by A1, A2, A3, A4, and A5. The actual measurement specifications of these specimens are shown in Table 1. The rock sample model is shown in Figure 1a, and the processed specimen is shown in Figure 1b. The water gushing from the tunnel mainly comes from the rich water in front of the tunnel face, so this special test piece shape is designed. The first section of the water injection hole of the rock sample is bonded together with the iron block of the fixed abrasive tool using strong adhesive, and then the fixed abrasive tool is reinforced by electric welding to resist the water pressure in the water injection hole after the water pump is running. The water injection pipe and the abrasive tool are tightened and fixed by screws to ensure that there is no water leakage on the water injection side during the experiment. The following figure shows the fixed mold of the rock sample. from the same rock sample. The specimen design specification is a cuboid of 100 × 100 × 150 mm. The diameter and depth of the observation hole in the test block are 50 mm, and the diameter and depth of the water injection hole are 50 mm and 50 mm. A total of five specimens were prepared, represented by A1, A2, A3, A4, and A5. The actual measurement specifications of these specimens are shown in Table 1. The rock sample model is shown in Figure 1a, and the processed specimen is shown in Figure 1b. The water gushing from the tunnel mainly comes from the rich water in front of the tunnel face, so this special test piece shape is designed. The first section of the water injection hole of the rock sample is bonded together with the iron block of the fixed abrasive tool using strong adhesive, and then the fixed abrasive tool is reinforced by electric welding to resist the water pressure in the water injection hole after the water pump is running. The water injection pipe and the abrasive tool are tightened and fixed by screws to ensure that there is no water leakage on the water injection side during the experiment. The following figure shows the fixed mold of the rock sample.

(**a**) rock sample structure (**b**) rock sample picture


**Table 1.** Actual measurement specifications of the specimen. **Table 1.** Actual measurement specifications of the specimen.

#### *3.3. Experiment Process 3.3. Experiment Process*

The rock specimen was loaded uniaxially with a closure rate of 0.1 mm/min. The data acquisition frequency of the testing machine was set as 10 times/s. The layout of the experimental equipment is shown in Figure 2. To facilitate the sorting and analysis of test data, the water pressure of the water pump was set as 0.5 MPa before the test, and the thermal imager was installed about 1 m away from the sample to observe the infrared temperature field changes on the sample surface. The experiment was started after the infrared radiation temperature on the rock surface remained stable. We synchronously calibrated the time of all test equipment. In addition, the start and end times of each equipment remained the same. Then, uniaxial loading was applied to the rock sample until inrush water appeared. The rock specimen was loaded uniaxially with a closure rate of 0.1 mm/min. The data acquisition frequency of the testing machine was set as 10 times/s. The layout of the experimental equipment is shown in Figure 2. To facilitate the sorting and analysis of test data, the water pressure of the water pump was set as 0.5 MPa before the test, and the thermal imager was installed about 1 m away from the sample to observe the infrared temperature field changes on the sample surface. The experiment was started after the infrared radiation temperature on the rock surface remained stable. We synchronously calibrated the time of all test equipment. In addition, the start and end times of each equipment remained the same. Then, uniaxial loading was applied to the rock sample until inrush water appeared.

**Figure 2.** Experimental layout. **Figure 2.** Experimental layout.

#### **4. Indicators 4. Indicators**

#### *4.1. Strain Energy 4.1. Strain Energy*

Assuming that a rock unit deforms under the action of external forces and the process occurs in a closed system, according to the first law of thermodynamics, the following can be obtained: Assuming that a rock unit deforms under the action of external forces and the process occurs in a closed system, according to the first law of thermodynamics, the following can be obtained:

$$
\mathcal{U}\mathcal{U} = \mathcal{U}^d + \mathcal{U}^e \tag{4}
$$

*U* = *Ud* +*Ue* (4) where *U* is the total strain energy, which is determined by the stress–strain curve and the area around the horizontal axis, *Ud* is the dissipated strain energy of the unit, and *Ue* is the where *U* is the total strain energy, which is determined by the stress–strain curve and the area around the horizontal axis, *U<sup>d</sup>* is the dissipated strain energy of the unit, and *U<sup>e</sup>* is the elastic strain energy released by the unit.

elastic strain energy released by the unit.

Figure 3 shows the relationship between *U<sup>d</sup>* and *U<sup>e</sup>* during uniaxial rock loading. The releasable elastic strain energy *U<sup>e</sup>* in the uniaxial loading process of rock can be rewritten as follows [41,42]: modulus 2 1 *<sup>f</sup> ( ) f( ) <sup>E</sup>* ε− ε <sup>=</sup> ε −ε (6)

1 <sup>2</sup>

*u*

*E*

In order to facilitate calculation, the elastic modulus *E*0 is generally used instead of *Eu*. In this paper, the average modulus is used to calculate the elastic modulus, and the stress–strain formula is assumed to be *σ* = *f* (*ε*). The following formula obtains the elastic

2

*e*

*U*

$$\mathcal{U}^{\varepsilon} = \frac{1}{2E\_{\mathfrak{u}}} \sigma^{2} \tag{5}$$

in the uniaxial loading process of rock can be rewritten

= σ (5)

during uniaxial rock loading. The

where *E<sup>u</sup>* is the unloading modulus of elasticity. the starting and ending points of the elastic phase, respectively.

where *Eu* is the unloading modulus of elasticity.

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 6 of 18

Figure 3 shows the relationship between *Ud* and *Ue*

releasable elastic strain energy *Ue*

as follows [41,42]:

**Figure 3.** Relationship between dissipated strain energy and elastic strain energy in the stress–strain curve. **Figure 3.** Relationship between dissipated strain energy and elastic strain energy in the stress–strain curve.

*4.2. Infrared Thermal Image*  Infrared thermal image is a series of object surface temperature distribution images output by an infrared thermal imager. The two-dimensional temperature matrix of frame *P* in the original infrared thermal image is [43]: In order to facilitate calculation, the elastic modulus *E*<sup>0</sup> is generally used instead of *Eu*. In this paper, the average modulus is used to calculate the elastic modulus, and the stress–strain formula is assumed to be σ = *f* (ε). The following formula obtains the elastic modulus

$$E = \frac{f(\varepsilon\_2) - f(\varepsilon\_1)}{\varepsilon\_2 - \varepsilon\_1} \tag{6}$$

where *P* is the frame number index of the infrared thermal image sequence; *x* and *y* represent the row and column numbers of the thermal imager temperature matrix, where *E* is the elastic modulus, and *f* (ε1) and *f* (ε2) are the stress values corresponding to the starting and ending points of the elastic phase, respectively.

#### respectively. *4.2. Infrared Thermal Image*

*4.3. VSMIT VSMIT* can reflect the dispersion degree of infrared radiation temperature value of Infrared thermal image is a series of object surface temperature distribution images output by an infrared thermal imager. The two-dimensional temperature matrix of frame *p* in the original infrared thermal image is [43]:

the entire rock sample surface, which is defined as follows [43]:

$$f\_p(\mathbf{x}, \mathbf{y}) \tag{7}$$
 
$$f\_{\mathbf{x}, \mathbf{y}, \mathbf{z}', \mathbf{z}, \mathbf{z}', \mathbf{z}'} = \mathbf{x} \tag{8}$$

2

1 1 *p p y x VSMIT f x, y AIRT M N* = = = − (8) where *p* is the frame number index of the infrared thermal image sequence; *x* and *y* represent the row and column numbers of the thermal imager temperature matrix, respectively.

#### *4.3. VSMIT*

*VSMIT* can reflect the dispersion degree of infrared radiation temperature value of the entire rock sample surface, which is defined as follows [43]:

$$VSMIT = \frac{1}{M} \frac{1}{N} \sum\_{y=1}^{N} \sum\_{x=1}^{M} \left[ f\_p(x, y) - AIRT\_p \right]^2 \tag{8}$$

Among them, *AIRT<sup>p</sup>* = <sup>1</sup> *M* 1 *N N* ∑ *y*=1 *M* ∑ *x*=1 *fp*(*x*, *y*); *M* and *N* are the maximum numbers of rows and columns for *x* and *y*, respectively.

### *4.4. IRVC*

The *IRVC* can measure the dispersion degree of the infrared radiation temperature field on the rock surface, which is defined as follows [39]:

$$IRVC = \sigma / AIRT \tag{9}$$

where *σ* is the standard deviation of infrared radiation temperature on the rock surface.

### **5. Experimental Results**

#### *5.1. Strain Energy*

The deformation and failure of rock is a process of energy input, elastic energy accumulation, energy dissipation, and energy release from the point of view of strain energy. Energy dissipation is mainly used for crack initiation and propagation, and energy release is the internal cause of the sudden failure of the rock mass. The elastic strain energy accumulated by the rock mass before excavation is the main source of the energy released by the ultimate failure of the surrounding rock mass, especially the deep hard brittle rock mass with good energy storage under the condition of high in situ stress. A large amount of elastic strain energy accumulated in the excavation is released instantly due to the excavation unloading effect, promoting the occurrence of rock mass failure and then connecting the water diversion fissure channel, finally resulting in a water inrush accident.

Figure 4 shows the strain energy evolution curve during the loading process of the rock sample. Due to space limitations, samples A1, A2, and A<sup>3</sup> were selected for analysis in this paper. As shown in Figure 4, the energy absorbed in the rock at the beginning of loading is mainly dissipated strain energy, because most of the strain energy is consumed by pore and microfracture compaction. The curves of elastic strain energy and dissipated strain energy diverge with the increase of stress, and the growth of elastic strain energy increases continuously, while the dissipated strain energy increases in a nearly straight line, which is used for the formation and expansion of plastic deformation and microcracks in rocks. The elastic strain energy drops sharply, while the dissipated strain energy increases sharply at the peak stress, which indicates that the internal microcrack propagation and penetration rate accelerated and the damage was aggravated. After the peak stress, the bearing capacity of the rock decreases rapidly and maintains a certain residual strength. The strain energy absorbed by the rock is transformed into dissipative strain energy in this period, which is used for the further development of rock fracture and shear deformation along the slip surface. Then, a macroscopic water diversion channel is formed in the rock, and water inrush eventually occurs in the roadway.

The dissipated strain energy ratio of rock refers to the proportion of dissipated strain energy to the total strain energy. The dissipated strain energy ratio curve shows a trend of decline before the peak stress, changes from a decline to a level near the peak stress when approaching the peak stress, and then begins to increase. The curve changes abruptly with the rock failure. After that, the dissipated strain energy ratio curve continues to increase until a water inrush occurs in the rock face. The analysis found that the strain energy dissipation ratio curve of samples experienced a "decline-level" process when the rock approached failure. This is due to the dissipated strain energy used for plastic strain and crack growth increasing with the rapid development of microcracks in the rock; although the elastic strain energy is still accumulating, the rock will reach the maximum energy storage limit. Therefore, the dissipated strain energy ratio of rock declines to the level, and then the elastic strain energy reaches the energy storage limit. The accumulated elastic strain energy is quickly released and causes rock failure. In summary, the turning point of dissipated energy ratio from decline to level can be used as the medium warning information of roadway water inrush.

summary, the turning point of dissipated energy ratio from decline to level can be used

as the medium warning information of roadway water inrush.

**Figure 4.** Evolution curves of strain energy during rock sample loading. **Figure 4.** Evolution curves of strain energy during rock sample loading.

#### *5.2. Infrared Radiation 5.2. Infrared Radiation 5.2. Infrared Radiation*

Figure 5 shows the photo of the water inrush instant of the roadway in the laboratory. The water diversion channel is formed after the rock in the roadway reaches failure, and then water inrush occurs in the mine. Hence, the hole is selected as the analysis area of infrared radiation data, as shown in Figure 6. The evolution characteristics of *AIRT*, *IRCV*, *VSMIT*, and infrared thermal image in rock holes during the rock loading failure and water inrush were analyzed, and the precursors of roadway water inrush were identified. Figure 5 shows the photo of the water inrush instant of the roadway in the laboratory. The water diversion channel is formed after the rock in the roadway reaches failure, and then water inrush occurs in the mine. Hence, the hole is selected as the analysis area of infrared radiation data, as shown in Figure 6. The evolution characteristics of *AIRT*, *IRCV*, *VSMIT,* and infrared thermal image in rock holes during the rock loading failure and water inrush were analyzed, and the precursors of roadway water inrush were identified. Figure 5 shows the photo of the water inrush instant of the roadway in the laboratory. The water diversion channel is formed after the rock in the roadway reaches failure, and then water inrush occurs in the mine. Hence, the hole is selected as the analysis area of infrared radiation data, as shown in Figure 6. The evolution characteristics of *AIRT*, *IRCV*, *VSMIT,* and infrared thermal image in rock holes during the rock loading failure and water inrush were analyzed, and the precursors of roadway water inrush were identified.

**Figure 5.** Water inrush in roadway. **Figure 5.** Water inrush in roadway. **Figure 5.** Water inrush in roadway.

**Figure 6.** Infrared radiation analysis region. **Figure 6.** Infrared radiation analysis region. **Figure 6.** Infrared radiation analysis region.

#### 5.2.1. AIRT 5.2.1. AIRT 5.2.1. AIRT

Figure 7 shows the time-varying curves of *AIRT*, *IRCV,* and *VSMIT* in the hole during the uniaxial rock sample loading. As shown in Figure 7, the *AIRT* in the hole of rock under loading showed a trend of stable fluctuation with the increase of stress during the process of rock uniaxial loading. The *AIRT* curves of rock samples A1 and A2 showed no obvious abnormal phenomena before rock failure, while the *AIRT* of the rock sample A3 gradually increased, with a temperature rise of about 0.1 °C. However, the *AIRT* of rock samples A1, A2, and A3 all dropped abruptly before water inrush, with a decrease temperature range of between 0.3~0.6 °C. This is due to the water seepage into the observation surface of the roadway absorbing part of the heat from the rock, resulting in a downward trend of *AIRT*. Therefore, the sudden drop of *AIRT* can be regarded as the precursor of water inrush. Two conditions must be satisfied to take a characteristic of infrared radiation index as a precursor: one is that all infrared radiation indexes of rock samples have this characteristic; the other is that this characteristic is easy to distinguish. The lead time of water inrush precursor of rock *AIRT* is 15~30 s before water inrush, hence *AIRT* decreased slightly and then increased. Figure 7 shows the time-varying curves of *AIRT*, *IRCV,* and *VSMIT* in the hole during the uniaxial rock sample loading. As shown in Figure 7, the *AIRT* in the hole of rock under loading showed a trend of stable fluctuation with the increase of stress during the process of rock uniaxial loading. The *AIRT* curves of rock samples A1 and A2 showed no obvious abnormal phenomena before rock failure, while the *AIRT* of the rock sample A3 gradually increased, with a temperature rise of about 0.1 °C. However, the *AIRT* of rock samples A1, A2, and A3 all dropped abruptly before water inrush, with a decrease temperature range of between 0.3~0.6 °C. This is due to the water seepage into the observation surface of the roadway absorbing part of the heat from the rock, resulting in a downward trend of *AIRT*. Therefore, the sudden drop of *AIRT* can be regarded as the precursor of water inrush. Two conditions must be satisfied to take a characteristic of infrared radiation index as a precursor: one is that all infrared radiation indexes of rock samples have this characteristic; the other is that this characteristic is easy to distinguish. The lead time of water inrush precursor of rock *AIRT* is 15~30 s before water inrush, hence *AIRT* decreased slightly and then increased. Figure 7 shows the time-varying curves of *AIRT*, *IRCV*, and *VSMIT* in the hole during the uniaxial rock sample loading. As shown in Figure 7, the *AIRT* in the hole of rock under loading showed a trend of stable fluctuation with the increase of stress during the process of rock uniaxial loading. The *AIRT* curves of rock samples A<sup>1</sup> and A<sup>2</sup> showed no obvious abnormal phenomena before rock failure, while the *AIRT* of the rock sample A<sup>3</sup> gradually increased, with a temperature rise of about 0.1 ◦C. However, the *AIRT* of rock samples A1, A2, and A<sup>3</sup> all dropped abruptly before water inrush, with a decrease temperature range of between 0.3~0.6 ◦C. This is due to the water seepage into the observation surface of the roadway absorbing part of the heat from the rock, resulting in a downward trend of *AIRT*. Therefore, the sudden drop of *AIRT* can be regarded as the precursor of water inrush. Two conditions must be satisfied to take a characteristic of infrared radiation index as a precursor: one is that all infrared radiation indexes of rock samples have this characteristic; the other is that this characteristic is easy to distinguish. The lead time of water inrush precursor of rock *AIRT* is 15~30 s before water inrush, hence *AIRT* decreased slightly and then increased.

**Figure 7.** Evolution curves of *AIRT*, *IRCV*, and *VSMIT* during rock sample loading. A–G: **Figure 7.** Evolution curves of *AIRT*, *IRCV*, and *VSMIT* during rock sample loading. A–G: different stages.

different stages.

#### 5.2.2. IRCV

The *IRCV* curve of rock sample A<sup>1</sup> shows a trend of steady fluctuation–decline–rise in the process of water inrush experiment under uniaxial loading, and the curve of *IRCV* mutates when loaded to 123 s, while the *IRCV* curves of rock samples A<sup>2</sup> and A<sup>3</sup> show a nearly horizontal trend, and both of them have a mutation in the middle of loading, with the occurrence moments of 94 s and 115 s, respectively. As shown in Figure 7, the *AIRT* value corresponding to the mutation of *IRCV* in the middle and late stages of loading decreases gradually, so the mutation of *IRCV* is caused by the mutation of infrared radiation standard deviation (*σ*). The authors propose that this is due to the rock having just entered the stage of unstable crack development, and microfracture events increase. The rock failure is dominated by microcracks induced by tensile failure. The tensile failure area corresponds to a drop in the rock surface temperature, while the shear failure area corresponds to a rise in the rock surface temperature, resulting in a gradual drop in *AIRT*. The increase of microfracture events leads to the occurrence of heating and cooling zones in different regions of the rock surface, and thus the standard deviation of corresponding infrared radiation temperature is suddenly changed. In conclusion, the curve of *IRCV* has suddenly changed in the middle and late stages of rock loading. Therefore, the *IRCV* mutation in the middle and late stages of rock loading can be regarded as the early precursor of water inrush.

As shown in Figure 7, when the *AIRT* of the rock drops suddenly before the water inrush (the precursor of water inrush), the *IRCV* curve of the rock increases or mutates gradually. Specifically, the *IRCV* curve of the rock sample A<sup>1</sup> increases gradually, and the rock samples A<sup>2</sup> and A<sup>3</sup> have a mutation. This feature, therefore, is not suitable as a precursor of roadway water inrush. If the water flows homogeneously into the roadway before water inrush, and the water has an amplification effect on the infrared radiation of rock [34], *AIRT* will drop abruptly, and the dispersion degree of infrared radiation temperature (*σ*) may also correspond to a sudden drop, which may cause no mutation in *IRCV*. If the water flows into the roadway inhomogeneously before water inrush, the *AIRT* will drop sharply and the differentiation of infrared radiation temperature (*σ*) will increase sharply, which will cause *IRCV* mutation. To sum up, *IRCV* curve mutation before water inrush is not universal, so it is not suitable as a precursor of roadway water inrush.

#### 5.2.3. VSMIT

The *VSMIT* curve of all the rock samples shows a general horizontal trend during the roadway water inrush test, and the *VSMIT* increases abruptly when the rock failure, with the mutation range, is 0.01~0.03. Due to the universality, synchronism, and significance of *VSMIT* mutation characteristics at the rock failure [26], the first mutation of *VSMIT* can be regarded as a midterm precursor of water inrush. With the process of rock loading, the stress gradually drops to the residual stress. When the *AIRT* of the rock declines abruptly (the precursor of water inrush), the *VSMIT* of all rock samples mutate for the second time, with an increased range of 0.02~0.05. This is due to the water seepage increasing the dispersion degree of infrared radiation temperature of the two adjacent frames, which also indicates that the sudden drop in *AIRT* will be accompanied by the sudden increase of the *VSMIT* index. Therefore, the second mutation of *VSMIT* can be used as a precursor of water inrush.

#### 5.2.4. Infrared Thermal Image

The above infrared radiation characteristics during the process of rock failure and water inrush were all obtained from the quantitative index analysis of infrared radiation, so only the time information of infrared radiation characteristics can be obtained. If the spatial information of infrared radiation during the process of rock failure and water inrush is needed, the infrared thermal image can be analyzed. To highlight the change of infrared radiation caused by the loading of the sample, and reduce the impact of the local radiation rate difference and environmental interference of the sample when processing the infrared

radiation experimental data, the thermal image obtained during the loading process is processed as a difference [31], that is, the first thermal image at the beginning of the loading is subtracted from each thermal image, and the change of radiation temperature field is analyzed by using the image after the difference.

As shown in Figure 7, the infrared thermal image of rock sample A<sup>1</sup> before failure shows the changing trend of bright (A)–dark (B)–bright (C)–dark (D), corresponding to the decline (AB)–rise (BC)–decline (CD) of the AIRT curve, and the temperature difference between the left and right sides of the sample is obvious at the rock failure (the temperature on the left side is high while that on the right side is low), and the infrared thermal image at the precursory point (F) of the roadway water inrush becomes dark as a whole. At the beginning of the loading stage (AB) of sample A2, the infrared thermal image of the rock sample changes from dark to bright, and there is no obvious abnormal change from point B to rock failure (E). The lower part of the rock shows the abnormal low-temperature area near the roadway water inrush precursory point (F), while the upper part of the infrared thermogram shows the abnormal high-temperature area when the roadway water inrush occurs, and the rest is the low-temperature area. At the beginning of loading stage (AB) of rock sample A3, there is no obvious change in the infrared thermogram, and the infrared thermogram becomes dark at the BC stage, while the abnormal hightemperature area appears at the lower part of the rock sample (D). With the increase of stress, the abnormal high-temperature area becomes an abnormal low-temperature area in the roadway water inrush precursor (E), and the high-temperature area appears on the left side. The high-temperature area extends upward with the loading process, and part of the original high-temperature area is eroded by water, so the right side of the rock sample presents a large area of low temperature.

#### **6. Multiparameter Normalization**

*6.1. Define*

If the value range of the sample data is [Min, Max], then the normalized expression of linear function is

$$y = (\text{x} - \text{Min}) / (\text{Max} - \text{Min}) \tag{10}$$

where *x* and *y* are the values before and after conversion, respectively; Max and Min are the maximum and minimum values of samples, respectively.

The linear function normalization has the following properties: (1) the sample size relation remains unchanged; (2) the relative distance of samples remains unchanged. The variation trend of each physical quantity obtained by linear function normalization is consistent with the original data curve, which can well reflect the key information such as fluctuation and mutation points in the original data curve [44].

#### *6.2. Analysis of Normalized Results*

Based on the analysis of strain energy and infrared radiation characteristics of roadway water inrush, it is found that dissipative energy ratio, *AIRT*, *VSMIT*, and *IRCV* are suitable as the main parameters to predict roadway water inrush. In the practical application of multiparameter joint monitoring of a roadway water inrush disaster, due to the differences in the range and dimension of each parameter, and the different physical parameters in different coordinate systems, it is not conducive to the rapid and intuitive identification of the sequence of the occurrence of the abrupt point of each physical parameter, and this also affects the analysis of the correlation between each physical parameter. Therefore, it is necessary to normalize each parameter in order to comprehensively compare and analyze the variation rule of each parameter in the same scale, which can provide a comprehensive basis for the early warning of roadway water inrush disasters.

Figure 8 shows the normalized curves of stress, dissipated energy ratio, *AIRT*, *VSMIT*, and *IRCV* with time collected in a laboratory experiment of sandstone roadway water inrush. Based on the comprehensive analysis of all rock samples, it can be found that *IRCV* mutation occurs at the stage of 0.55~0.65 TWI (TWI is the time of water inrush) in the

middle loading stage, which is the early precursor of water inrush in a sandstone roadway (as shown the EPWI in Figure 8). The dissipated energy ratio curve drops to the lowest point with the increase of stress, which is the first middle precursor of roadway water inrush (as shown by the FMPWI in Figure 8), and the early precursor is 20~33 s earlier than the middle precursor one. As the loading continues, rock sample failure occurs and accompanies the mutation of *VSMIT*, which is the second middle precursor of roadway water inrush (as shown by SMPWI in Figure 8). Thereafter, the stress drops rapidly and *AIRT* decreases abruptly at about 0.95~0.98 TWI after it decreases to residual strength. This is the first precursor of roadway imminent water inrush (as shown by FPIWI in Figure 8), corresponding to the sudden increase of *VSMIT*, which is the second precursor of imminent water inrush (as shown by SPIWI in Figure 8), indicating that roadway water inrush may occur at any time. that *IRCV* mutation occurs at the stage of 0.55~0.65 TWI (TWI is the time of water inrush) in the middle loading stage, which is the early precursor of water inrush in a sandstone roadway (as shown the EPWI in Figure 8). The dissipated energy ratio curve drops to the lowest point with the increase of stress, which is the first middle precursor of roadway water inrush (as shown by the FMPWI in Figure 8), and the early precursor is 20~33 s earlier than the middle precursor one. As the loading continues, rock sample failure occurs and accompanies the mutation of *VSMIT*, which is the second middle precursor of roadway water inrush (as shown by SMPWI in Figure 8). Thereafter, the stress drops rapidly and *AIRT* decreases abruptly at about 0.95~0.98 TWI after it decreases to residual strength. This is the first precursor of roadway imminent water inrush (as shown by FPIWI in Figure 8), corresponding to the sudden increase of *VSMIT*, which is the second precursor of imminent water inrush (as shown by SPIWI in Figure 8), indicating that roadway water inrush may occur at any time.

scale, which can provide a comprehensive basis for the early warning of roadway water

Figure 8 shows the normalized curves of stress, dissipated energy ratio, *AIRT*, *VSMIT*, and *IRCV* with time collected in a laboratory experiment of sandstone roadway water inrush. Based on the comprehensive analysis of all rock samples, it can be found

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 13 of 18

inrush disasters.

1.0 **Figure 8.** *Cont*.

0.0

**7. Discussion** 

0.2

0.4

0.6

Multiparameter normalization

0.8

AIRT

Dissipation energy ratio

IRCV VSMIT

0 50 100 150 200

(1) Due to the nonlinear process of roadway water inrush, the complexity and diversity of influencing factors, and the control of the accuracy of monitoring technology, the prediction of roadway water inrush with a single parameter has great limitations. The precursor of roadway water inrush cannot be accurately and effectively identified in practical application, which may lead to the problem of false alarms or missed alarms. The precursors of roadway water inrush can be identified more quickly by the normalized treatment of each physical parameter in the process of roadway water inrush. At the same time, the correlation between precursor information of roadway water inrush was obtained. Comprehensively considering

Time (s)

Sample A3 **Figure 8.** Multiparameter normalization curve in the process of roadway water inrush.

EPWI

Stress

Water inrush

FPIWI

FMPWI SMPWI

SPIWI

0 60 120 180

Time (s)

EPWI

Stress

VSMIT

IRCV

Dissipation energy ratio

AIRT

FPIWI

SMPWI

FMPWI

Water inrush

SPIWI

**Figure 8.** Multiparameter normalization curve in the process of roadway water inrush. **Figure 8.** Multiparameter normalization curve in the process of roadway water inrush.

#### **7. Discussion 7. Discussion**

0.0

0.2

0.4

0.6

Multiparameter normalization

0.8

1.0

	- (2) *AIRT* reflects the whole infrared radiation intensity of the rock surface, but there may be different heating and cooling zones in the process of rock loading and fracture due to which the *AIRT* remains unchanged. The *IRCV* of rocks reflects the dispersion degree of the original infrared radiation temperature, which has the advantage of avoiding the unit of measurement of data and neglecting the influence of numerical magnitude. Compared with *AIRT* index, *IRCV* can reflect the dispersion characteristics of infrared radiation caused by temperature rise and temperature drop areas. *VSMIT* reflects the dispersion degree of the difference in infrared radiation temperature between two adjacent frames. Compared with *IRCV* index, it eliminates the cumulative heating effect of loaded rock, and is easier to monitor the process of rock failure, instability, and seepage. From the sensitivity of *IRCV* index to the unstable crack development stage of rock, *IRCV* mutation was proposed as the early precursor of roadway water inrush. Based on the feature that *VSIMT* can monitor rock failure, the first mutation of *VSIMT* was proposed as the medium-term precursor of roadway water inrush. Established from the characteristic that *AIRT* and *VSIMT* are sensitive to water, the sudden drop of *AIRT* and the second mutation of *VSIMT* were proposed as the precursor of roadway imminent water inrush. In this paper, combined with the advantages of *AIRT*, *IRCV,* and *VSMIT*, the multiparameters precursory characteristics of roadway water inrush were determined.
	- (3) Infrared observation technology is a promising new method for monitoring rock samples, which has the advantages of noncontact and strong anti-interference, and can

be used for monitoring and warning the stability of bearing rock and surrounding rock of tunnels in underground engineering. Acoustic emission monitoring technology can detect the time and location of a microfracture in a rock mass. Therefore, acoustic emission monitoring (internal) and infrared radiation monitoring (external) should be combined in subsequent roadway water inrush experiments and underground engineering construction sites. Acoustic emission (AE) will be used to locate the waterconducting fissure passage in the rock, and the change rules and coupling effects of the stress field, infrared radiation temperature field, and seepage field before water inrush in underground engineering will be studied in order to build a multifield coupling model of "stress–temperature–seepage" of roadways based on infrared radiation and reveal the mechanism of water inrush in underground engineering.

### **8. Conclusions**

The forecast of highway water inrush with a single parameter has significant limits because of the nonlinear nature of the roadway water inrush process, the complexity and variety of contributing elements, and the control of the accuracy of monitoring equipment. The issue of false alerts or missed alarms may arise from the inability to precisely and efficiently identify the antecedent of highway water inrush in practical application. The normalized treatment of each physical parameter in the process of highway water inrush may help identify the antecedents of roadway water inrush more promptly. The association between the precursor data of the highway water inrush was discovered concurrently. The following conclusions were drawn:


**Author Contributions:** K.C.: conceptualization, writing—original draft, and experiments; F.D.: conceptualization, supervision, and project administration; L.M.: supervision and project administration; N.M.K.: software and data curation, writing—original draft; methodology; T.F.: visualization; S.S.A.: investigation and data curation; S.H.: Formal analysis; M.A.: experiments. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by Anhui Provincial Scientific Research Preparation Plan Project (2022AH050596). The authors also acknowledge the Researchers Supporting Project number (RSP2023R496), King Saud University, Riyadh, Saudi Arabia.

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

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding author upon reasonable request.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


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

**Jiayi Shen 1,2, Chenhao Sun <sup>1</sup> , Huajie Huang <sup>1</sup> , Jiawang Chen <sup>3</sup> and Chuangzhou Wu 1,2,\***


**Abstract:** Scale effects on the mechanical behavior of rock joints have been extensively studied in rocks and rock-like materials. However, limited attention has been paid to understanding scale effects on the shear strength of rock joints in relation to normal stress *σ*<sup>n</sup> applied to rock samples under direct shear tests. In this research, a two-dimensional particle flow code (PFC2D) is adopted to build a synthetic sandstone rock model with a standard joint roughness coefficient (JRC) profile. The manufactured rock model, which is adjusted by the experiment data and tested by the empirical Barton's shear strength criterion, is then used to research scale effects on the shear strength of rock joints caused by normal stresses. It is found that the failure type can be affected by JRC and *σ*n. Therefore, a scale effect index (SEI) that is equal to JRC plus two times *σ*<sup>n</sup> (MPa) is proposed to identify the types of shear failure. Overall, shearing off asperities is the main failure mechanism for rock samples with SEI > 14, which leads to negative scale effects. It is also found that the degree of scale effects on the shear strength of rock joints is more obvious at low normal stress conditions, where *σ*n < 2 MPa.

**Keywords:** rock joints; shear strength; scale effects; normal stress; JRC; PFC simulation

### **1. Introduction**

Rock joints play an important role in the estimation of the shear strength of rock masses [1–3]. Effective design of rock engineering projects, such as underground excavations and open pit slopes, requires precise estimation of the shear strength of rock joints [4–6]. However, it is well known that there is a scale effect on the shear strength of rock joints [7–9]. The main difficulty in determining how the shear strength of rock joints varies with scale is conducting expensive and time-consuming engineering scale in situ testing [10]. Although laboratory scale tests on a small jointed sample cannot generate the precise shear strength of rock joints, they can still reveal the mechanical behavior of jointed rock masses [11]. Therefore, laboratory tests are widely used by researchers to investigate how the shear strength of rock joints is affected by sample sizes. Table 1 presents a review of scale effects on the shear strength of rock joints, which shows conflicting results. The majority of results show that there is a negative scale effect on the shear strength, which means the shear strength decreases with the increase of joint sizes. Some results [12,13] show positive scale effects, which represent the shear strength increases when the joint size increases. While other results [13–15] show no scale effects.

Scale effects on the shear strength of rock joints could be explained in different ways. One explanation is that scale effects occur due to the contact area of joints changing with the increase in joint size [16]. Pratt et al. [14] and Yoshinaka et al. [17] attributed the decrease in shear strength to the smaller contact area of the sample where higher stress was concentrated on these contact surfaces. The other explanation is that the scale effect is associated with the change of undulations and asperities on a joint surface as joint

**Citation:** Shen, J.; Sun, C.; Huang, H.; Chen, J.; Wu, C. Scale Effects on Shear Strength of Rough Rock Joints Caused by Normal Stress Conditions. *Sustainability* **2023**, *15*, 7520. https:// doi.org/10.3390/su15097520

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

Received: 3 March 2023 Revised: 12 April 2023 Accepted: 26 April 2023 Published: 4 May 2023

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

size increases. A longer sample will result in higher undulation amplitude compared to a smaller sample [18]. Barton and Choubey [19] concluded that the shear behavior of larger rock samples is governed by larger and gentler asperities, while the steep and small asperities are the controlling mechanism in smaller rock samples. Giani et al. [20] stated that when the rock joint shear strength depends on the random distribution of asperity, it will produce a positive scale effect. If the rock joint shear strength depends on wavy and rough surfaces, then there is a negative scale effect. Therefore, more research is required to determine the exact nature of the scale effect on the shear strength of joints.


**Table 1.** Review of scale effects on the shear strength of rock joints [21].

"N" means negative scale effect; "P" means positive scale effect; "O" means no scale effect.

Based on the literature review, we noticed that the existing laboratory tests were carried out under various normal stress conditions ranging from 0.0245 to 10 MPa, as shown in Table 1. As we know, the failure mode of rock joints during the direct shear test will be affected by the application of normal stresses. When rock samples are under high normal stress conditions, the tips of asperities could be sheared off; therefore, the shear strength would be relatively higher compared to rock samples that are under low normal stress conditions where sliding is the controlling mechanism of rock failure.

Numerical simulations using the PFC are capable of simulating the asperity damage and degradation process during the shearing tests [28]. It has been proven that the shear strength results acquired from PFC modeling are typically comparable with experimental test results [29]. Therefore, PFC simulations as an alternative to physical testing can be used to reveal the fundamental mechanism of shear behavior of rock joints at various scales.

In this research, a synthetic rock model based on the two-dimensional particle flow code (PFC2D)-based synthetic sandstone rock model is used to study the influence of normal stress on scale effects on the shear strength of rock samples with standard joint roughness coefficient (JRC) profiles, and attempts to answer two questions: (1) Are scale effects on shear behavior affected by normal stresses? (2) What is the degree of scale effects affected by normal stresses?

In this paper, the synthetic rock model for numerical tests is introduced in Section 2. The verification of the synthetic rock model is shown in Section 3. Scale effect investigations are presented and discussed in Sections 4 and 5.

#### **2. Synthetic Rock Model for Numerical Tests**

*2.1. Synthetic Rock Model Based on PFC2D*

PFC2D is a discrete element program. The bonded particle model (BPM), a composite of rounded particles, simulates complete rock and does not require a continuum-scale constitutive model to depict the mechanical behavior of intact rock [30]. The parallel bond model, which can replicate the physical behavior of a substance similar to cement linking the two nearby particles, is one of the most fundamental and often used BPMs in the PFC2D, as illustrated in Figure 1. Bond breaking reduces stiffness because contact and bond stiffness both contribute to stiffness in a parallel bond model. While contact stiffness is active as long as particles are in contact, bond stiffness is instantly gone when a bond breaks [31]. Therefore, the parallel bond model is a more accurate bond model for materials that resemble rocks, since it allows for the possibility of bonds breaking in tension or shearing with a corresponding loss in stiffness.

**Figure 1.** Illustration of the parallel bond model.

On the other hand, by adding joints to a BPM assembly using the smooth joint model (SJM), jointed rock masses can be created. The BPM's original contact microscopic characteristics will be replaced with SJM properties with the names friction coefficient *µ*<sup>j</sup> , shear stiffness *k*sj, and normal stiffness *k*nj when the SJM is put into the BPM [32]. The synthetic rock model constructed by the BPM and SJM has the ability to simulate various mechanical responses of jointed rock masses including peak strength [31], scale effect [33], anisotropy [34], and cracking processes [30,33] in rocks and rock-like materials.

The numerical direct shear test used in this research is presented in Figure 2. The specimen (40 mm × 100 mm) is generated using the BPM. The rock joint is created using the SJM. During the direct shear test, the upper block receives the normal force in a vertical direction. The upper block is given a horizontal velocity of 0.03 m/s, while the lower block is held in place.

**Figure 2.** Numerical direct shear test set up.

#### *2.2. Calibration of Numerical Models*

In this study, the synthetic rock model was calibrated using laboratory data from Australia's Hawkesbury Sandstone [28]. Firstly, tests for uniaxial compression on rock samples (42 mm × 84 mm) with a loading rate of 0.02 m/s were performed to determine the BPM's parameters after a calibration process [31] to ensure that the mechanical properties of the synthetic rock model are close to laboratory data.

It should be mentioned that one of the key factors influencing the resilience of restricted materials to deformation and strength, such as rocks and cemented soil, is the loading rate [35–38]. In this research, we did not consider such loading rate effects on the mechanical properties of jointed rocks.

Using the calibrated BPM parameters indicated in Table 2 to perform the uniaxial compression test (Figure 3), the values of elastic modulus *E*, Poisson's ratio *v*, and UCS produced are comparable to experimental tests, as shown in Table 3.

**Table 2.** Micro-parameters of the BPM model.


**Figure 3.** Numerical uniaxial compressive test.

**Table 3.** Comparison of mechanical properties calculated from the numerical model and tested from laboratory.


Then, synthetic rock models (40 mm × 100 mm) with planar joints were constructed. The SJM has the following micro-parameters: friction coefficient *µ*<sup>j</sup> , shear stiffness *k*sj, and normal stiffness *k*nj. In this research, the values of *k*nj = 25 GPa, *k*sj = 13 GPa, and *µ*<sup>j</sup> = 0.75 were selected using the inverse-modeling calibration approach to ensure that the numerical rock model can give a similar response as that from laboratory tests with joint shear stiffness *K*<sup>s</sup> = 6.4 GPa/m, normal stiffness *K*<sup>n</sup> = 28.6 GPa/m, and joint friction angle *ϕ*<sup>b</sup> = 37.6◦ . The calibration procedure was as follows: (1) The normal deformability compression test was carried out to calibrate normal stiffness *k*nj. (2) The shear test was carried out to calibrate shear stiffness *k*sj under normal stress of 1 MPa condition. (3) Direct shear tests were undertaken and friction coefficient *µ*<sup>j</sup> was calibrated. Figures 2–4 present the final mechanical responses of the synthetic rock models after the final calibration.

**Figure 4.** The normal deformability test on the synthetic rock specimen with a planar joint.

Figure 4 shows the axial stress-displacement curves of the synthetic rock specimen (40 mm × 100 mm) under the normal deformability test with the loading rate of 0.02 m/s. The value of joint normal stiffness *K*<sup>n</sup> generated by the synthetic rock specimen is 28.6 GPa/m, which is close to the laboratory test results with *K*<sup>n</sup> = 28.8 GPa/m.

Figure 5 shows the shear stress-displacement curve of the synthetic rock model (40 mm × 100 mm) with a planar joint under a direct shear test (loading rate of 0.03 m/s) with the normal stress of 1 MPa. The value of joint shear stiffness *K*<sup>s</sup> generated by the synthetic rock specimen is 6.4 GPa/m, which is the same as laboratory test results with *K*<sup>s</sup> = 6.4 GPa/m.

**Figure 5.** The direct shear test on the synthetic rock specimen with a planar joint under normal stress of 1 MPa.

Figure 6 shows the failure envelope of the synthetic rock model (40 mm × 100 mm) with a planar joint under direct shear tests. The value of joint friction angle *ϕ*<sup>b</sup> generated by the synthetic rock specimen is 36.1◦ , which is close to laboratory test results with *ϕ*<sup>b</sup> = 37.6◦ .

*Sustainability* **2023**, *14*, x FOR PEER REVIEW 6 of 15

**Figure 6.** Shear strength of the synthetic rock samples under various normal stresses. **Figure 6.** Shear strength of the synthetic rock samples under various normal stresses.

#### **3. Validation of Synthetic Rock Models**

**3. Validation of Synthetic Rock Models**  To confirm the reliability of the synthetic rock model shown in Section 2, direct shear To confirm the reliability of the synthetic rock model shown in Section 2, direct shear tests on the synthetic rock models with 10 standard JRC profiles were performed and the shear strength values produced from numerical simulations were compared to those derived from Barton's empirical shear strength model.

### tests on the synthetic rock models with 10 standard JRC profiles were performed and the *3.1. Barton's Shear Strength Model*

shear strength values produced from numerical simulations were compared to those derived from Barton's empirical shear strength model. *3.1. Barton's Shear Strength Model*  One of the most widely adopted empirical strength criteria for estimating rock joint shear strength in rock engineering is the Barton's shear strength criterion. Based on the results of a large number of shearing tests on various rock joint profiles, Barton and his co-workers [19,39] proposed an empirical equation to estimate the shear strength of rock joints, as shown in Equation (1).

$$\tau = \sigma\_\mathbf{n} \tan \left[ \varphi\_\mathbf{b} + \text{[RC/g} \left( \frac{\text{JCS}}{\sigma\_\mathbf{n}} \right) \right] \tag{1}$$

We performed extensive numerical direct shear experiments on synthetic rock

failure envelopes generated by direct shear tests on synthetic rock models were compared to the empirical Barton's shear strength criterion (Equation (1)) with *φ*b = 36.1° and JCS = 27.4 MPa. Figure 7 compares the shear strength obtained from numerical simulations to Barton's model, indicating that the usage of synthetic rock models is

results of a large number of shearing tests on various rock joint profiles, Barton and his co-workers [19,39]proposed an empirical equation to estimate the shear strength of rock where *ϕ*<sup>b</sup> is the joint friction angle. JCS is the joint compression strength, which is equal to UCS of intact rock in this research. JRC stands for joint roughness coefficient and can be calculated using standard joint profiles.

#### joints, as shown in Equation (1). *3.2. Numerical Simulation Results*

=୬tan [ୠ + JRC lg (JCS ୬ )] (1) where *φ*b is the joint friction angle. JCS is the joint compression strength, which is equal to UCS of intact rock in this research. JRC stands for joint roughness coefficient and can be We performed extensive numerical direct shear experiments on synthetic rock models with varied JRC profiles in normal stress levels between 0.5 and 5 MPa. The failure envelopes generated by direct shear tests on synthetic rock models were compared to the empirical Barton's shear strength criterion (Equation (1)) with *ϕ*<sup>b</sup> = 36.1◦ and JCS = 27.4 MPa. Figure 7 compares the shear strength obtained from numerical simulations to Barton's model, indicating that the usage of synthetic rock models is capable of generating adequate shear strength of rock joints.

capable of generating adequate shear strength of rock joints.

calculated using standard joint profiles.

*3.2. Numerical Simulation Results* 

**Figure 7.** Comparison failure envelopes obtained from numerical simulations and the Barton's empirical model.

#### **4. Configuration of Rock Samples for Scale Effect Investigations**

Two methods are widely used for investigating scale effects on the shear strength of rock joints [15]. The first one is to divide a large rock joint into several smaller sizes of rock joints, as shown in Figure 8a, which presents an example of the division of the Barton's typical profile. The geometry characteristics or the values of JRC of smaller sizes of rock joints can be different from that of the original larger rock joint. The other method is the assembly of several repeated 100 mm profiles into larger rock joints many times the original profile length, as shown in Figure 8b. The joint roughness or the value of JRC of assembly samples is the same as that of the original joint surface.

**Figure 8.** Two types of configuration of samples for scale effect investigations (**a**) division of the Barton's JRC profile (**b**) assembly of repeated Barton's JRC profile.

In fact, the scale effect on the shear strength of rock joints includes two factors, which are the sample size itself and the geometrical characteristics of the joint surface. Rock samples generated by the division method have various sample sizes and geometry characteristics. However, the rock joints generated by the assembly method have various sample sizes but have the same geometry characteristics. It is well known that the geometry characteristics will affect the shear strength of rock joints [8]. Therefore, in this research, we adopt the assembly of a repeated model that has the same geometry characteristics and JRC values to research the influence of pure sample size on the shear strength of rock joints.

#### **5. Results and Discussion**

Once synthetic rock models were validated, a series of rock specimens with various JRC profiles and different sizes (40 mm × 100 mm, 80 mm × 200 mm, and 120 mm × 300 mm) were generated to study the effect of sample sizes on rock joint shear strength. The shear strength values of different sizes of rock samples under given normal stresses were calculated and are summarized in Figure 9.

In this research, the index *k* (see Equation (2)), which is the average slope of three points, was used to identify the types of scale effects.

$$k = \frac{N\sum x\_i y\_i - \sum x\_i \sum y\_i}{N\left(\sum x\_i^2\right) - \left(\sum x\_i\right)^2} \tag{2}$$

where *x<sup>i</sup>* is the joint length of the rock sample, *y<sup>i</sup>* the shear stress of the rock sample, and *N* is the number of the testing sample. *k* > 0 means the rock joint has a positive scale effect and *k* < 0 means the rock joint has a negative scale effect. The value of *k* can be calculated using three groups of data. For example, for rock samples with JRC = 2 under the normal stress *σ*<sup>n</sup> = 5 MPa, the shear strength of rock samples with joint lengths *l* = 100, 200, and 300 mm are 4.2, 4.4, and 4.6 MPa, respectively. Therefore, data (100, 4.2), (200, 4.4), and (300, 4.6) were put into Equation (2) to calculate the value of *k*. The result shows *k* = 0.4, which means the scale effect is positive. Table 4 shows comprehensive scale effect results of rock samples with various JRC profiles under different normal stress conditions. In Table 4, P means positive scale effect and N means negative scale effect.

The results presented in Table 4 are also plotted in Figure 10. It is found that the failure mode of rock joints during the direct shear test will be affected by the applying normal stresses and the joint roughness. When rock samples under high normal stress conditions, the tips of asperities with large joint roughness coefficient could be sheared off, therefore, the number of shear crack is relatively higher compared with rock samples with small joint roughness coefficients under low normal stress conditions where sliding is

the controlling mechanism of rock failure. For example, when a rock sample with JRC = 2 under the normal stress *σ*<sup>n</sup> = 0.5 MPa, the number of shear cracks is 10 and the scale effect is positive. However, when a rock sample with JRC = 20 under the normal stress *σ*<sup>n</sup> = 5 MPa, the number of shear cracks is 280 and the scale effect is negative. *Sustainability* **2023**, *14*, x FOR PEER REVIEW 9 of 15

**Figure 9.** Results of scale effects on the shear strength of rock samples under various normal stress conditions. **Figure 9.** Results of scale effects on the shear strength of rock samples under various normal stress conditions.

<sup>−</sup> <sup>=</sup>

In this research, the index *k* (see Equation (2)), which is the average slope of three

*N xy x y <sup>k</sup> Nx x*

where *xi* is the joint length of the rock sample, *yi* the shear stress of the rock sample, and *N*  is the number of the testing sample. *k* > 0 means the rock joint has a positive scale effect and *k* < 0 means the rock joint has a negative scale effect. The value of *k* can be calculated

( ) ( )<sup>2</sup> <sup>2</sup> *ii i i*

−

*i i*

(2)

points, was used to identify the types of scale effects.


**Table 4.** Results of scale effects of shear strength of rock joints.

**Figure 10.** Failure pattern and crack number of rock samples (40 mm × 100 mm) at peak shear strength.

Based on the results of Table 4, a scale effect index (*SEI*) which is equal to *JRC* plus two times of normal stress (MPa), as shown in Equation (3), was proposed to identify the types of scale effects.

$$SEI = fRC + \mathcal{D}\sigma\_n \tag{3}$$

The values of *SEI* for rock samples under different normal stress conditions are given in Table 4. For example, for rock samples with *JRC* = 2 under the normal stress *σ*<sup>n</sup> = 0.5 MPa, the value of *SEI* = 2 + 2 × 0.5 = 3. The number P3 in Table 4 means the rock sample with *SEI* = 3 has a positive scale effect. It was found that 20 out of 21 rock samples have negative scale effects when *SEI* > 14, and 29 out of 33 rock samples have positive scale effects when *SEI* < 14.

To find the possible reason why the use of SEI can identify the types of scale effects on shear strength, we monitored the crack number generated in the synthetic rock sample when the stress reaches the peak strength during the direct shear tests. When the parallel link between nearby particles in the PFC rock model is broken, a micro-tensile crack or micro-shear crack can occur. Figure 10 shows the failure pattern corresponding to the shear crack number of each rock sample when the shear stress reaches the peak strength. For example, the S = 10 represents a shear crack number of 10 and T = 7 represents a tension crack number of 7 for a sample with SEI = 3 (JRC = 2 and *σ*<sup>n</sup> = 0.5 MPa). The relations between SEI values and shear crack numbers of rock samples are also plotted in Figure 11. We can find that the number of shear cracks is low when SEI < 14. However, the number of shear cracks dramatically increases when the value of SEI is over 14, where the controlling failure mechanism transforms sliding to shearing off asperities. **Figure 10.** Failure pattern and crack number of rock samples (40 mm × 100 mm) at peak shear strength.

However, the number of shear cracks dramatically increases when the value of SEI is over 14, where the controlling failure mechanism transforms sliding to shearing off asperities.

*Sustainability* **2023**, *14*, x FOR PEER REVIEW 11 of 15

**Figure 11.** Relations between SEI and shear crack number of rock samples. **Figure 11.** Relations between SEI and shear crack number of rock samples.

Therefore, we can conclude that the results presented in Table 4 and Figure 11 show that the proposed SEI is capable of identifying types of scale effects. When SEI < 14, sliding over joints is the controlling mechanism of rock failure, which leads to positive scale effects; however, shearing off asperities could be the controlling mechanism of rock failure for rock samples with SEI > 14, which leads to negative scale effects. Therefore, we can conclude that the results presented in Table <sup>4</sup> and Figure <sup>11</sup> showthat the proposed SEI is capable of identifying types of scale effects. When SEI < 14, sliding over joints is the controlling mechanism of rock failure, which leads to positive scale effects; however, shearing off asperities could be the controlling mechanism of rock failure for rock samples with SEI > 14, which leads to negative scale effects. *Sustainability* **2023**, *14*, x FOR PEER REVIEW 12 of 15

> On the other hand, to further identify the degree of scale effects, the coefficient of variance (CV), which can calculate the value of Standard Deviation/Mean to quantify the random influence of a bunch of data, was further used as a scale effect magnitude index On the other hand, to further identify the degree of scale effects, the coefficient of variance (CV), which can calculate the value of Standard Deviation/Mean to quantify the random influence of a bunch of data, was further used as a scale effect magnitude index to quantify scale effects on shear strength of rock joints caused by normal stress conditions. The calculation results are presented in Figure 12. to quantify scale effects on shear strength of rock joints caused by normal stress conditions. The calculation results are presented in Figure 12.

**Figure 12.** The degree of scale effect caused by normal stresses. **Figure 12.** The degree of scale effect caused by normal stresses.

It is found that, in general, the values of CV decrease with the increase of normal stress. We also calculated the average value of CV for a given group of data. Figure 12 shows that the values of average CV decrease when the normal stress increases from 0.5 to 2.0 MPa, then, it tends to be stable with further increase of normal stresses, which means the degree of scale effects on shear strength of rock joints is more obvious at low normal stress conditions where *σ*n < 2 MPa. It is found that, in general, the values of CV decrease with the increase of normal stress. We also calculated the average value of CV for a given group of data. Figure 12 shows that the values of average CV decrease when the normal stress increases from 0.5 to 2.0 MPa, then, it tends to be stable with further increase of normal stresses, which means the degree of scale effects on shear strength of rock joints is more obvious at low normal stress conditions where *σ*<sup>n</sup> < 2 MPa.

Such a phenomenon can also be validated by the laboratory data published by Fardin [24], who carried out a laboratory study of the scale effect on the shear strength of

The CV values of samples under a specific normal stress condition were calculated and are shown in Figure 14. The value of CV is up to 0.4 when *σ*n = 1 MPa, then, it decreases sharply to 0.14 when *σ*n increases to 2 MPa. After that, there is a slight change in CV values with a further increase of *σ*n from 2 MPa to 10 MPa. Such change in CV values with

**Figure 13.** Laboratory data of scale effect on shear stress of rock joints under various normal stress

conditions (data from Fardin [24]).

0 50 100 150 200 250

Rock sample size (mm)

0

4

8

Shear stress (MPa)

12

16

20

normal stresses is similar to that of the numerical results in Figure 12.

10MPa 5MPa 2.5MPa 1MPa

Such a phenomenon can also be validated by the laboratory data published by Fardin [24], who carried out a laboratory study of the scale effect on the shear strength of concrete replicas with roughness joints. Laboratory test results are shown in Figure 13. The CV values of samples under a specific normal stress condition were calculated and are shown in Figure 14. The value of CV is up to 0.4 when *σ*<sup>n</sup> = 1 MPa, then, it decreases sharply to 0.14 when *σ*<sup>n</sup> increases to 2 MPa. After that, there is a slight change in CV values with a further increase of *σ*<sup>n</sup> from 2 MPa to 10 MPa. Such change in CV values with normal stresses is similar to that of the numerical results in Figure 12. Fardin [24], who carried out a laboratory study of the scale effect on the shear strength of concrete replicas with roughness joints. Laboratory test results are shown in Figure 13. The CV values of samples under a specific normal stress condition were calculated and are shown in Figure 14. The value of CV is up to 0.4 when *σ*n = 1 MPa, then, it decreases sharply to 0.14 when *σ*n increases to 2 MPa. After that, there is a slight change in CV values with a further increase of *σ*n from 2 MPa to 10 MPa. Such change in CV values with normal stresses is similar to that of the numerical results in Figure 12.

*Sustainability* **2023**, *14*, x FOR PEER REVIEW 12 of 15

conditions. The calculation results are presented in Figure 12.

**Figure 12.** The degree of scale effect caused by normal stresses.

01234567

Normal stress *σ*n (MPa)

to quantify scale effects on shear strength of rock joints caused by normal stress

JRC = 2 JRC = 4 JRC = 6 JRC = 8 JRC = 10 JRC = 12 JRC = 14 JRC = 16 JRC = 18 JRC = 20 Average

It is found that, in general, the values of CV decrease with the increase of normal stress. We also calculated the average value of CV for a given group of data. Figure 12 shows that the values of average CV decrease when the normal stress increases from 0.5 to 2.0 MPa, then, it tends to be stable with further increase of normal stresses, which means the degree of scale effects on shear strength of rock joints is more obvious at low

Such a phenomenon can also be validated by the laboratory data published by

normal stress conditions where *σ*n < 2 MPa.

0.00

0.05

0.10

CV 0.15

0.20

0.25

**Figure 14.** Laboratory investigation on the degree of scale effect caused by normal stresses. **Figure 14.** Laboratory investigation on the degree of scale effect caused by normal stresses.

#### **6. Conclusions 6. Conclusions**

available.

Synthetic rock models with standard JRC profiles were constructed in PFC2D to investigate scale effects on the shear strength of rock joints under various normal stress conditions. The capability of the synthetic rock model to simulate the shear behavior of rock joints was tested by comparing numerical simulations with the Barton's shear Synthetic rock models with standard JRC profiles were constructed in PFC2D to investigate scale effects on the shear strength of rock joints under various normal stress conditions. The capability of the synthetic rock model to simulate the shear behavior of rock joints was tested by comparing numerical simulations with the Barton's shear strength criterion.

strength criterion. Once synthetic rock models were validated, a series of rock specimens of different sizes (40 mm × 100 mm, 80 mm × 200 mm, and 120 mm × 300 mm) were generated to investigate the influence of sample sizes on rock joint shear strength under normal stress Once synthetic rock models were validated, a series of rock specimens of different sizes (40 mm × 100 mm, 80 mm × 200 mm, and 120 mm × 300 mm) were generated to investigate the influence of sample sizes on rock joint shear strength under normal stress ranges from 0.5 to 5 MPa.

ranges from 0.5 to 5 MPa. Numerical simulation results show that the types of scale effects could be affected by the JRC profiles and normal stresses. Therefore, a scale effect index (SEI) that is equal to JRC plus two times normal stress (MPa), as shown in Equation (3), was proposed to identify the types of scale effects. It is found that for the rock sample with SEI < 14, sliding over joints is the controlling mechanism of rock failure, which leads to positive scale effects. However, shearing off asperities could be the controlling mechanism of rock Numerical simulation results show that the types of scale effects could be affected by the JRC profiles and normal stresses. Therefore, a scale effect index (SEI) that is equal to JRC plus two times normal stress (MPa), as shown in Equation (3), was proposed to identify the types of scale effects. It is found that for the rock sample with SEI < 14, sliding over joints is the controlling mechanism of rock failure, which leads to positive scale effects. However, shearing off asperities could be the controlling mechanism of rock failure for rock samples with SEI > 14, which leads to negative scale effects.

We also further investigated the influence of normal stress on the degree of scale effects on the shear strength of rock joints. It is discovered that the degree of scale effect is

Finally, it should be noted that the finding of this research is based on the analysis of test data of Australia Hawkesbury Sandstone. Therefore, the finding of this research is open for further improvements as more shear strength data of various rock types become

**Author Contributions:** Conceptualization, J.S.; methodology, C.S. and H.H.; software, J.S. and H.H.; validation, J.S. and H.H.; formal analysis, J.S. and H.H.; investigation, J.S. and C.S.; resources, J.S. and C.W.; data curation, C.S. and H.H.; writing—original draft preparation, J.S. and H.H.; writing—review and editing, J.S., C.S., J.C. and C.W.; visualization, C.S. and H.H.; supervision, J.S.

**Funding:** This research was funded by the Zhoushan Science and Technology Bureau (2022C81001), the Key R&D of Zhejiang Province (2021C03183), and the Scientific Research Fund of Zhejiang University (XY2021011), the Sanya Yazhou Bay Science and Technology City (KYC-2020-01-001), the Finance Science and Technology Project of Hainan Province (ZDKJ202019).

**Data Availability Statement:** The data presented in this study are available on request from the

**Acknowledgments:** We express our sincere gratitude to the reviewers for their valuable comments

and C.W. All authors have read and agreed to the published version of the manuscript.

failure for rock samples with SEI > 14, which leads to negative scale effects.

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

**Informed Consent Statement:** Not applicable.

corresponding author.

and suggestions.

We also further investigated the influence of normal stress on the degree of scale effects on the shear strength of rock joints. It is discovered that the degree of scale effect is more obvious at low normal stresses conditions where *σ*<sup>n</sup> < 2 MPa.

Finally, it should be noted that the finding of this research is based on the analysis of test data of Australia Hawkesbury Sandstone. Therefore, the finding of this research is open for further improvements as more shear strength data of various rock types become available.

**Author Contributions:** Conceptualization, J.S.; methodology, C.S. and H.H.; software, J.S. and H.H.; validation, J.S. and H.H.; formal analysis, J.S. and H.H.; investigation, J.S. and C.S.; resources, J.S. and C.W.; data curation, C.S. and H.H.; writing—original draft preparation, J.S. and H.H.; writing—review and editing, J.S., C.S., J.C. and C.W.; visualization, C.S. and H.H.; supervision, J.S. and C.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Zhoushan Science and Technology Bureau (2022C81001), the Key R&D of Zhejiang Province (2021C03183), and the Scientific Research Fund of Zhejiang University (XY2021011), the Sanya Yazhou Bay Science and Technology City (KYC-2020-01-001), the Finance Science and Technology Project of Hainan Province (ZDKJ202019).

**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.

**Acknowledgments:** We express our sincere gratitude to the reviewers for their valuable comments and suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


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