#### Abstract

There is a coupling relationship between surrounding rock stress, deformation, and fracture evolution, especially in the microdynamics of the crust caused by mining activities and earthquakes. Previous research has investigated many cases regarding the coseismal water level responses and proposed a method to calculate the aquifer parameters by tidal analysis. However, to date, measurement of the degree of rock damage in the field has not been reported. Quantifying the fracture characteristics is essential for accurate evaluation of rock stability. This study has analyzed the relationship between the seismograms and hydroseismograms in response to the *Mw* 7.8 Solomon Islands earthquake and the *M*w 7.8 Kaikōura earthquake, both events occurring in 2016. The calculated and measured changes in water level in the X10 well were fitted in order to study the relationships among the volumetric strain, the deviatoric strain, and the oscillations in the pore pressure. Then, we further estimate the degree of rock damage and the hydraulic characteristics of the aquifer. The results showed that the values for the rock damage parameter, 0.662 < *α*_{D} < 0.754, and the Skempton coefficient, −0.100 < *A* < 0.026, estimated for the Solomon Islands earthquake signified higher damage and dilatancy in the X10 well. Also, the respective values for the parameters, 0.293 < *α*_{D} < 0.363 and 0.226 < *A* < 0.251, calculated for the Kaikōura earthquake signified a lower degree of rock damage. It is concluded that the changes in the pore pressure were influenced by both the volumetric strain and the deviatoric strain. The degree of rock damage and the hydraulic properties of the aquifer estimated from the water level fluctuations in the wells which were induced by the seismic waves represent the actual aquifer characteristics.

#### 1. Introduction

Understanding the mechanical and hydrological properties of rock is important for engineering, including mining, sources evaluation, landslide stability assessment, infrastructure stability, reservoir geomechanics, and so on [1–4]. Fractures significantly control the properties of rock, and thus quantifying fracture characteristics is essential for accurate evaluation of rock stability. Usually, an empirical damage-mechanism model is introduced and is widely used in geophysics and engineering for irreversible deformation [5]. The rock damage is described by 0 < *α*_{D} < 1 which is influenced by aquifer properties and stress state [6]. High damage value (*α*_{D} > 0.8) means that the rock is highly fractured and its mechanical properties are low while its permeability is high. Low damage value (*α*_{D} < 0.2) means stiff and strong rock with low permeability. However, measuring the rock damage in the field is still impossible.

Water level influenced by barometric pressure, Earth tides, and earthquakes indicates that aquifer is sensitive to the crustal strain which provides a way to probe the material characteristics [7–12]. It is generally recognized that the water level in aquifers responds to volumetric strain [13, 14]. Moreover, the water level responses in aquifer system caused by earthquakes may be explained by the dynamic interactions between the seismic waves and the aquifers [15, 16]. The Rayleigh and waves are the seismic waves that have been documented most widely as being responsible for the expansion and contraction of aquifers [17–19]. However, recent studies have shown that the water level of well also responds to *S* and Love waves because of the anisotropic poroelastic characteristics [15, 20]. Laboratory studies on the effect of shear stress on pore pressure demonstrated that the deviatoric strain significantly influenced the pore pressure, and there was a strong correlation as the degree of rock damage increased [21, 22]. Recently, Shalev et al. [23] showed that the water level in wells responding to deviatoric strain depended on the rock damage, and the deviatoric strain of the damaged rock led to high amplitudes in the water level. For rocks without damage, there was no response from the water level in the well [23]. Therefore, it should, in principle, be possible to determine quantitatively the degree of damage to the rock from the responses of the water level to the deviatoric strain. In this study, the relationship between the seismograms and the hydroseismograms in response to the *M*w 7.8 Solomon Islands earthquake (2016-12-9) and the *M*w 7.8 Kaikōura earthquake (2016-11-13) will be investigated. The X10 well in Xinjiang is selected for study because it is located at the intersection of two faults with extensive rock fractures. The rock damage and the hydraulic properties of the aquifer in the well will be estimated from the fluctuations in the water level induced by the seismic waves.

#### 2. Observational Background and Hydrogeological Setting

The groundwater well (X10) was established to the south of Urumqi City in Xinjiang of China by the China Earthquake Administration in September 1980 to monitor the possible reservoir-induced seismicity. It is located at the intersection of Liushugou-Hongyanchi fault (43.70°N, 87.62°E; with the elevation of 1056 m above sea level) (Figure 1). The depth of the open well is 28 m. The aquifer is located in the depth between 16 and 24 m, consisting of sandstone and muddy siltstone of Permian Hongyanchi group. It has a static water of approximately 1 m below ground surface. Extensive fractures could be found, and a screen section was installed.

#### 3. Seismic Data

The water level in the well was recorded using an SWY-II type digital meter from October 2016 to the present with an accuracy of ±0.05% F.S. and a precision of ±0.02% F.S. in the range of 0–50 m [16]. The sampling interval was 1 s, and the information recording capability was a significant improvement. Furthermore, an SLJ-100 type three-component seismograph was used for collecting the radial, transverse, and vertical components of the seismograms following the earthquakes. The hydroseismograms and the seismograms following the *M*w 7.8 Solomon Islands earthquake and the *M*w 7.8 Kaikōura earthquake are shown in Figure 2. The epicenter of the Solomon Islands earthquake was 9,543 km southeast of the X10 well while the epicenter of the Kaikōura earthquake was 12,798 km southeast of the well.

The X10 well is connected indirectly to the aquifer through the fault. The water flow in the pores of the aquifer may flow indirectly to the well bore along the fracture channel of the fracture zone. Taking into account the arrival time of the P, S, Love, and Rayleigh waves, the data of 1 : 50–3:05 (GMT +8 : 00) (Solomon Islands earthquake) and 18 : 43–19 : 58 (GMT +8 : 00) (Kaikōura earthquake) were selected for this study. Spline interpolation was used to fill in the missing data, and the data were smoothed by a high-pass Butterworth filter with a bottom corner frequency of 0.01 Hz to remove the long period trends. The times before the arrival of the seismic waves—1:50 : 00–2:00 : 42 on 2016-12-9 (Solomon Islands earthquake) and 18 : 43 : 00–19 : 16 : 28 on 2016-11-13 (Kaikōura earthquake)—were regarded as the noninterference states. The periods of 2 : 00 : 43–2:10 : 42 or 2 : 10 : 43–2:20 : 00 for the Solomon Islands earthquake corresponding to P, S, Love, and Rayleigh waves, respectively, were regarded as the interference state, and the same was considered for the periods 19 : 16 : 29–19 : 27 : 37 or 19 : 27 : 38–19 : 38 : 04 in the case of the Kaikōura earthquake. When the long period waves of Love and Rayleigh arrived, the largest amplitudes for the level of the well water occurred [24]. Compared with the P and S period, the water level response was more distinct during the Love wave period and Rayleigh wave period. The vertical and radial seismograms were larger than those in the transverse direction. The shapes of radial, transverse, and vertical seismograms were similar to those of water level changes. Because of the imbalance of the total water head between the wellbore and the aquifer caused by air pressure, the response of the water level lagged behind that of the seismic waves. This phenomenon was attributed to the wellbore storage effect, which could cause attenuation of the water level, the water pressure, the phase shift, and the aquifer permeability [7, 8, 13, 17, 18, 25, 26]. Due to wellbore storage effect, the measured changes in the water level are not accurate. In order to maintain equilibration in the aquifer-borehole, when estimating the rock failure by changes in the water level, it is necessary to correct for the wellbore storage effect. Such a correction method has been reported by Liu et al. [18]. Also, the attenuation factor of an open well can be calculated from equations [4, 17, 18]. The attenuation factor depends mainly on the transmissivity, the storage coefficient, the aquifer thickness, and the effective height of the water column in the well-aquifer system. In this study, the attenuation factor of the X10 well was 2.049 as calculated for the Solomon Islands earthquake while the value was 8.967 for the Kaikōura earthquake.

#### 4. Methods

##### 4.1. Estimation of the Degree of Rock Damage by Seismic Analysis

The amplitudes of the oscillations in the water level and the changes in the pore pressure are affected by the wavelength and frequency of the specific seismic waves. Poroelastic theory has been used extensively to study the oscillations in the water level in response to earthquakes, barometric pressure, and Earth tides, especially with respect to the seismicity in oil-rich areas [27–29]. Several studies have successfully applied poroelastic theory in explaining the relationship between the well-aquifer system and seismic waves [16, 19, 23]. Poroelastic theory attested to the fact that the strain was affected by the both pore pressure and stress in an isotropic porous elastic medium [30]. According to Henkel and Wade [31], Holtz et al. [32], and Skempton [33], changes in pore pressure were also affected by the deviatoric strain [31–33]. It was reported that based on Biot’s theory, changes in the pore pressure could not be caused by deviatoric stress [34]. However, Skempton [33] allowed for this possibility under undrained conditions. A two-dimensional stress state was employed to describe the relationship between the deviatoric strain, the volumetric strain, and the pore pressure [23]:where △*P*_{f} represents the change in pore pressure, *K*_{u} is the undrained bulk modulus [35], *ε*_{v} is the volumetric strain, *N* is the coupling coefficient of the shear strain, *N* *=* −4 *μB*(*A-*1/3), *μ* is the shear modulus, and *ε*_{d} is the deviatoric strain where *ε*_{d}=(*ε*_{1} − *ε*_{3})/2. *A* and *B* are referred to as Skempton coefficients which relate to rock compressibility. *B* depends on the rock lithology in the range from 0 to 1. *A* may be expressed as follows:

*A* is related to *μ* and an additional modulus *γ* which couples the shear and volumetric stress of the nonlinear materials. *μ* = *μ*_{0} + *α*_{D}*ζ*_{0}*γ*, and *γ* = *α*_{D}, where *μ*_{0,}, and *ζ*_{0} are the elastic moduli. *ζ*_{0} is a critical value of the strain invariant ratio at the onset of damage accumulation and is connected to the internal friction angle of the friction law for rocks. *ζ*_{0} = −1, which is based on the typical value for sandstones [8]. The damage parameter *α*_{D} represents the elastic characteristics of a damaged solid.

Substituting *N* = -4 *μB*(*A*-1/3) into (2) gives the value for *N*:where the tidal coefficient *BK*_{u} is a constant. It is characterized by rock elasticity. The *α*_{D} value is equal to 1 for a material completely damaged and equal to 0 for a damage-free material.

The relationship between the volumetric strain, the deviatoric strain, and the pore pressure can be established by both tidal analysis and seismic analysis. Then, the rock damage may be calculated.

##### 4.2. Estimation of Parameters by Tidal Analysis

In a previous study by Doan and Cornet [36], the relationship between the phase shift and the poroelastic properties was used to determine *B* and *K*_{u}. In addition, the permeability and transmissivity may be estimated by tidal analysis [7, 37].

In a well, when the poroelastic medium of the rock changes due to external force, the water level will fluctuate. The relationship between volumetric strain and pore pressure may be expressed as follows:where *P* is the pore pressure, is the strain, *BK*_{u} is characterized by the rock elasticity, and *B* depends on the rock lithology (the value is in the range of 0 to 1).

Substituting (5) into (4), *BK*u may be expressed as follows:where is the tidal factor, which represents the fluctuations in the amplitude ratio of the water level to volumetric strain, *ρ* is the water density, *ρ* = 1 g/cm^{3}, is the acceleration due to gravity, *h* is the water level, and *ε* is the volumetric strain.

In an ideal, confined, and isotropic well-aquifer system, the parameters can also be estimated by tidal analysis [7, 38]:where is the tidal factor, △*h* is the change in water level, △*ε* is the change in the volumetric strain, and *S*_{s} is the specific storage. may be determined from tidal analysis. According to (7), the tidal factor and the specific storage obey a reciprocal relationship. Moreover, the relationship between *S*_{s} and *n* is as follows:where *n* is the porosity, *B*_{e} is the barometric efficiency calculated by a regression deconvolution method or Clark’s method [16, 25], and *β* is the water compressibility, *β* = 4.6 × 10^{−10 }m^{2}/N. Lithology dilatancy is an important characteristic for granules especially for dense sand. An increase in *n* represents an increase in the expansion of the lithology [39].

##### 4.3. Estimation of Parameters by Seismic Analysis

###### 4.3.1. Strain in Surface Waves

The actual state of rock failure during the deformation process is difficult to estimate. Hence, this research is meaningful for studying the formation of old ruptures and the generation of new ruptures. Love and Rayleigh waves propagate along the surface of the Earth, and they can change in their propagation direction and attenuate with depth. Also, stress can be estimated from the seismic wave velocities:where *u*_{R} is the radial displacement, is the radial seismic wave velocity, *u*_{T} is the transverse displacement, and is the transverse seismic wave velocity. *V*_{R} (3.5 km/s) is the velocity of the Rayleigh wave, and *V*_{L} (3 km/s) is the velocity of the Love wave, which depends on the particular phase [23, 40].

Combining (9), (10), and (1), the following formula is obtained:

When only the effect of the volumetric strain on the pore pressure is considered [23]:

When only the effect of the deviatoric strain on the pore pressure is considered:

###### 4.3.2. Strain in Body Waves

The volumetric strain and the deviatoric strain calculated for body waves are as follows:where *u*_{z} is the vertical displacement and is the vertical velocity^{23}. *V*_{p} (2.5 km/s) is the velocity of the wave and *V*_{s} (1.5 km/s) is the velocity of the *S* wave, which are dependent on the particular phase [23, 40]. Both phases are the acknowledged seismic waves which can lead to aquifer expansion and contraction [41].

Substituting (14) and (15) into (1) results in the following equation:

#### 5. Results

##### 5.1. Tidal Analysis

In this study, the Solomon Islands earthquake and the Kaikōura earthquake were selected. The time series selected was divided into 10-day interval and a 5-day sliding window. The tidal factors obtained by BAYTAP-G in each interval did not change with time, which indicated that the aquifer characteristics were stable. The tidal factor estimated for the Solomon Islands earthquake was 4.167 × 10^{5} m, and that for the Kaikōura earthquake was 4.831 × 10^{5} m.

In the Solomon Islands earthquake, *BK*_{u} was 4.167 GPa as calculated by (6). Considering the negative reciprocal relationship of the specific storage and the tidal factor, a value of *S*_{s} *=* 2.40 × 10^{−6} m^{−1} was obtained. Meanwhile, *n* = 0.27 was estimated. For the Kaikōura earthquake, *BK*_{u} = 4.831 GPa, *S*_{s} = 2.07 × 10^{−6} m^{−1}, and *n* = 0.22, as specified in Table 1.

##### 5.2. Seismic Analysis

###### 5.2.1. The Love and Rayleigh Phases

In this study, the time series of the Love and Rayleigh phases during the Solomon Islands earthquake (2 : 10 : 43–2:20 : 00) and the Kaikōura earthquake (19 : 27 : 38–19 : 38 : 04) were selected. Using the seismic velocities, the measured water level changes and the calculated water level changes were fitted. The results are listed in Table 2. When assuming only volumetric strain, , *BK*_{u} = 3.616 GPa with a correlation coefficient of *R*^{2} = 0.675 for the Solomon Islands earthquake (Figure 3(a)) and *BK*_{u} = 12.364 GPa with a correlation coefficient of *R*^{2} = 0.001 for the Kaikōura earthquake (Figure 3(b)). When assuming only deviatoric strain, *N* = 17.363 GPa with a low correlation coefficient value of *R*^{2}_{=} 0.258 for the Solomon Islands earthquake (Figure 3(c)) and *N* = 2.529 GPa with *R*^{2} = 0.557 for the Kaikōura earthquake (Figure 3(d)). The relatively smaller *R*^{2} value for the former indicated that only the volumetric or the deviatoric strain in the X10 well had small effects on the changes in the pore pressure. When considering both the volumetric strain and the deviatoric strain, , *BK*_{u} = 3.229 GPa, and *N* = 8.437 GPa with *R*^{2} = 0.728 (Solomon Islands earthquake) (Figure 3(e)) and , *BK*_{u} = 3.717 GPa, and *N* = 2.529 GPa with *R*^{2} = 0.560 (Kaikōura earthquake) (Figure 3(f)). The best correlation coefficient between the calculated and measured water pressure showed that the changes in the pore pressure were influenced by both the volumetric strain and the deviatoric strain.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

###### 5.2.2. The P and S Phases

When considering both the volumetric strain and the deviatoric strain of the *P* and *S* phases, the results were *BK*_{u} = 0.397 GPa and *N* = 7.384 GPa with *R*^{2} = 0.554 for the Solomon Islands earthquake (Figure 3(g)) and *BK*_{u} = 3.197 GPa and *N* = 2.238 GPa with *R*^{2} = 0.201 for the Kaikōura earthquake (Figure 3(h)). The amplitude during the *P* and *S* periods was small, and many phases were contained in the same packet, so the estimated strain by the *P* and *S* waveforms was not well resolved. Surface waves can cause significant expansion of the aquifer compared to that of the body waves. Therefore, the pressure fluctuations were caused mainly by surface waves, and the parameters should be calculated with reference to the Love and Rayleigh phases.

#### 6. Discussion

##### 6.1. Estimation of Parameters

The obtained *BK*_{u} values for the two earthquakes by tidal analysis were similar (4.167 GPa for the Solomon Islands earthquake and 4.831 GPa for the Kaikōura earthquake) while the respective values were 3.229 GPa (Solomon Islands earthquake) and 3.717 GPa (Kaikōura earthquake) by seismic analysis. Zhang and Huang [42] showed that the *BK*_{u} values for sandstone measured in the laboratory ranged from 0.039 to 26.40 GPa. Also, Liao et al. [43] reported the measurement of *BK*_{u} values at different effective pressures in overseas trials. These independent studies demonstrated that the results of this study were consistent with real-world data. In addition, according to the experimental data for Berea Sandstone under a mean effective stress of 10 MPa, *BK*_{u} varied from 8.5 to 15.5 GPa [30]. *BK*_{u} for the X10 well obtained in this study was very close to the range for sandstone in the experimental examples. Therefore, the rock damage and the aquifer properties estimated from the water level fluctuations in wells induced by seismic waves represent the actual aquifer characteristics. Considering the parameter inversion results and the different *BK*_{u} values from previous studies, it is concluded that there is good agreement between the seismic analysis and the tidal analysis, and this is in accordance with the experimental results for sandstone. The best correlation coefficient between the calculated and measured water pressure and the water level fluctuations indicated that the pore pressure was influenced by both the volumetric strain and the deviatoric strain.

##### 6.2. Estimation of Rock Damage

A linear poroelastic material without dilatancy has *A* = 1/3 and *N* = 0 and an undrained Poisson’s ratio *V*_{u} = 0.25. It is known that stress expansions exist in fractured rock [39]. For porous sandstone without cohesion under undrained conditions, with stress increasing, the sensitivity of pore pressure to deviatoric strain increases [21]. With accumulation of the rock damage, the coupling between the pore pressure and strain becomes increasingly apparent, and the shear stress coupling coefficient increases. This situation tends to appear in active tectonic areas which have high levels of stress. Thus, in combining the tectonic and geological conditions, the high stress coefficient values mean that the rocks around the well are close to damage [23]. In this study, *N* was calculated by (3). Based on the response of the pore pressure to stress, the damage parameter *α*_{D} and the Skempton coefficient *A* can be estimated by (2) (Figure 4).

**(a)**

**(b)**

From the results of seismic analysis, the shear stress coupling coefficient *N* = 8.437 GPa corresponded to *α*_{D} values of 0.662 to 0.754, and *A* values for this damage state were between −0.100 and 0.026 in the X10 well as estimated for the Solomon Islands earthquake (Figure 4(a)); also, *N* = 2.529 GPa which corresponded to respective values of 0.293<*α*_{D} < 0.363 and 0.226<*A*<0.251 as estimated for the Kaikōura earthquake (Figure 4(b)). The high shear stress on the faults signified that the strain was large enough to lead to remarkable water level oscillations.

##### 6.3. The Degree of Rock Damage

Shalev et al. [23] showed that the response of the water level to deviatoric strain depended on the rock damage and that the deviatoric strain of the damaged rock led to high amplitudes in the water level. In the undamaged state, there is no response [23]. In analytical studies, a coupling of *A* = 0.2 related to *α*_{D} = 0.45 near Parkfield was considered to be slightly below *A* = 1/3 which indicated that there was no obvious dilatancy in the linear poroelastic medium [44]. In contrast, in the Gome 1 well established on the Dead Sea Transform, *BK*_{u} values of 14 GPa and 11 GPa were calculated by tidal analysis and seismic analysis, respectively, for the condition *V*_{u} = 0.38–0.4. The final result was *α*_{D} = 0.84–0.9 [23]. According to an experiment whereby a deviatoric stress of 65 MPa with a confining pressure of 17 MPa was imposed to sandstones with *n* = 21% in Berea and a deviatoric stress of 68 MPa with a confining pressure of 17 MPa was imposed to sandstones with *n* = 22% in Navajo, the conclusion was that an increasing *A* meant an increase in the original rupture and the production of a new rupture[21].

In this research, it was found that stress would cause the porous elastic media parameters to change which signified that the rock was damaged. The results for *α*_{D} and *A* suggested different degrees of rock damage in the X10 well during the two earthquake periods: 0.662 < *α*_{D} < 0.754 and −0.100 < *A* < 0.026 values estimated for the Solomon Islands earthquake (2016-12-9) were indicative of higher damage and dilatancy in the linear poroelastic medium of the X10 well.

Sun et al. [16] proposed a methodology that exploited spectral analysis based on the response of the water level fluctuations to Rayleigh waves to calculate the hydraulic parameters. Using this method, the optimal hydraulic conductivity *K* was 324 m/d for the Solomon Islands earthquake and 69 m/d for the Kaikōura earthquake. The higher the degree of rock damage is, the higher would be the hydraulic conductivity. Laboratory tests have demonstrated that the loading-induced increase in the volumetric strain would cause the permeability to expand [45–52]. However, according to (7) and (8), a specific relationship between the rock damage and the specific storage or porosity was not found [16]. Low permeability rock may have very high porosity (>99%) [20]. Based on the parameters for the Solomon Islands earthquake, the results of 0.293 < *α*_{D} < 0.363, 0.226<*A*<0.251, *K* = 69 m/d, and *n* = 0.22 for the X10 well calculated for the Kaikōura earthquake (2016-11-13) had a lower degree of rock damage. On the day before the Solomon Islands earthquake, there was a *M*w 6.2 earthquake (Hutubi earthquake, 2016-12-08) in Changji county of Xinjiang Province, China. The epicenter was located 113 km west of the X10 well, and there were numerous aftershocks. This earthquake destroyed the structure of the strata and altered the hydrogeological conditions of the X10 well, especially the stress-strain conditions, which led directly to higher rock damage and higher hydraulic permeability. It is concluded that using the relationship between the seismic waves and the water waves to invert the degree of rock damage is a feasible proposition. For the same aquifer conditions, there is a connection between the degree of rock damage and the aquifer properties. The higher the degree of rock damage is, the higher will be the hydraulic conductivity, and thus a method to estimate the degree of rock damage and the properties of the aquifer from water level fluctuations in wells induced by seismic waves has been realized.

#### 7. Conclusions

A method for estimating the degree of rock damage and the hydraulic properties of the aquifer in the X10 well from the fluctuations in the water level induced by seismic waves from the *M*w 7.8 Solomon Islands earthquake (2016-12-09) and the *M*w 7.8 Kaikōura earthquake (2016-11-13) has been developed. Based on the curve fitting data, the best correlation coefficient between the calculated and measured water pressures showed that the changes in the pore pressure were influenced by both the volumetric strain and the deviatoric strain. By comparing the results with the tidal response, it was shown that the results represent the actual aquifer characteristics. In this work, the damage parameter 0.662 < *α*_{D} < 0.754 for the Solomon Islands earthquake was related to the shear strain coupling coefficient *N* = 8.437 GPa and the Skempton coefficient −0.100 < *A* < 0.026; in the case of the Kaikōura earthquake, the respective values were 0.293 < *α*_{D} < 0.363, *N* = 2.529 GPa, and 0.226 < *A* < 0.251. It was shown to be feasible to invert the value for the degree of rock damage by using the relationship between the seismic waves and the water waves. For the same aquifer conditions, there is a connection between the degree of rock damage and the hydraulic properties of the aquifer, namely, the higher the degree of rock damage is, the higher will be the hydraulic conductivity of the aquifer. Thus, characterization studies on the hydraulic properties of the aquifer may be used to study the degree of rock damage.

#### Data Availability

Requests for access to these data should be made to the corresponding author.

#### Disclosure

The original water level and seismic wave data supporting the work presented in this paper are from the monitoring centre of the Earthquake Agency of the Xinjiang Uygur Autonomous Region. These data are essential for this research.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Authors’ Contributions

Di Zhao was responsible for methodology, formal analysis, writing, review and editing, and visualization. Yifan Zeng was responsible for writing and formal analysis. Xiaolong Sun was responsible for formal analysis and visualization. Aoshuang Mei was responsible for review and editing.

#### Acknowledgments

This research was supported by the National Natural Science Foundation of China (42072284, 42027801, and 41877186 to Yifan Zeng; 41972253 to Xiaolong Sun) and the Major Science and Technology Projects of Inner Mongolia Autonomous Region (2020ZD0020-4).