Abstract
As a national development plan, ecological protection of the Yellow River Basin has attracted extensive social attention in recent years. Considering the influence of acid rain on the engineering characteristics of loess in this area, we investigated changes in the physical and mechanical characteristics of remolded loess under the combined action of acid rain and dry-wet cycles by means of triaxial tests, nuclear magnetic resonance spectroscopy, and scanning electron microscopy. The results are as follows: in the acidic environment, the stress-strain relationship of remolded loess undergoes stress hardening after dry-wet cycling. The cohesion and internal friction angle of remolded loess are negatively correlated with the number of cycles. From the multiscale analysis of the dry-wet cycle process under acid rain condition, the T2 spectrum of the test soil has three peaks at the micropore level. With the increasing number of cycles, the spectral area increases gradually, and the sample transitions from small pores to large- and medium-size pores. At the microscopic level, the clay mineral particles among soil particles decrease in size, the contact mode between soil particles develops from stable to unstable, the particles are gradually rounded, and the fractal dimension decreases. Chemical erosion and physical erosion are special features of this experiment. Physical erosion causes particle erosion and pore growth, while chemical erosion includes reactions by feldspar. Together, physical and chemical reactions aggravate the soil deterioration process. These research results have laid a good experimental foundation for the ecological protection of the Yellow River Basin.
1. Introduction
The ecological protection and high-quality development of the Yellow River Basin, a recent focus of national attention, is related to the national development plan. The Yellow River Basin is rich in coal, oil, natural gas, and nonferrous metal resources, including coal deposits that account for more than half of the national reserves. This region is an important source of energy, chemicals, and raw materials and a basic industrial base in China. However, the ecological environment of this region is fragile, river pollution is prominent, and soil erosion is serious [1]. Therefore, it is important to pay attention to the environmental and geotechnical problems in this area, especially because a large amount of fossil fuel combustion may bring acid rain.
The special rock and soil properties of the Yellow River Basin are responsible for its ecological fragility. Loess, a quaternary sediment, is widely distributed in the middle and upper reaches of the Loess Plateau [2]. Loess usually shows the characteristics of thick layer distribution mixed with paleosol interlayers [2], large intergranular pores [3], developed vertical joints [4], strong water sensitivity [5] etc., and loess is especially vulnerable to long-term water soil interactions caused by river erosion, rainfall infiltration, farmland irrigation, etc.
Research on soil-water interactions mainly focuses on the mechanical strength and pore distribution of soil under external conditions, such as seepage [6–8], water migration caused by freeze-thaw cycles [9, 10], and dry-wet cycles [11, 12]. Research on dry-wet cycles mainly focuses on the hydraulic behavior [13, 14], strength and deformation characteristics [15, 16], fracture evolution characteristics [17–19], and permeability strength characteristics [20] as well as the dynamic characteristics [21] and conductivity characteristics [22] of soil undergoing dry-wet cycles. The flow chart is shown in Figure 1. However, in recent years, the interactions between soil and water in an acid rain environment have become the focus of research.
An important means to study the interactions between acid rain and soil is the study of the leaching characteristics and release rules of phosphorus, calcium, magnesium, boron, and other major elements from samples [23, 24]. At the same time, the response characteristics of different soil ecosystems to acid rain, such as respiration [25, 26], are also the focus of attention. Many studies have been done on changes in engineering characteristics caused by the interaction between acid rain and soil; these studies investigate the shear strength, seepage characteristics, and pore characteristics of residual soil, loess, and expansive soil under the action of acid rain [10, 27, 28].
In summary, dry-wet cycling, a typical case of soil stress and deformation caused by soil-water interactions, has been studied extensively. However, solution properties, an important index of the test results, are not considered, especially as they relate to the macro- and microscale characteristics of soil under the combined action of acid rain erosion and dry-wet cycles. In view of this gap in the literature, this paper takes the environmental and geotechnical problems of the Yellow River Basin as background, studies the physical and mechanical changes in remolded loess for an acid rain environment and dry-wet cycles, and uses results obtained at the macro-, meso-, and microscale levels to explore damage laws of experimental rocks and soil. The research results are expected to provide a reference for ecological protection and prevention of soil erosion in the Yellow River Basin.
2. Test Process and Principle
2.1. On-Site Sampling
The city of Qingyang in Gansu Province, the sample collection site, is located in the hinterland of the Loess Plateau, has very serious soil erosion, and is an important energy base in China. According to long-term monitoring, a large amount of fossil fuel combustion causes acid rain in this area. Therefore, it is representative of conditions elsewhere in China to take samples in this area for testing. The test soil sample is Malan loess of the late Pleistocene epoch, with a yellowish-brown color and moderate soil moisture content. The sample was taken to the laboratory for basic physical and mechanical tests. The test process was in strict accordance with the standard for soil test methods (GB/T 50123-2019). The specific test parameters are shown in Table 1.
Table 1 shows that the water content and natural gravity of the loess samples are moderate, the water content of the soil is slightly larger in the plastic state, and the soil is hard. To better reflect the effects of soil erosion caused by acid rain, the mineral composition of the test soil samples should be tested, and the main components of the test soil samples are quartz, calcite, albite, anorthite, muscovite, kaolinite, and hydrotalcite.
2.2. Sample Preparation and Selection
Remodeled soils were selected as test samples, their density was set to 1.76 g/cm3, and their moisture content was 16%. During the process of sample preparation, the soil samples were rolled, dried, and sieved, and ordinary tap water was added to form soil samples with a moisture content of 16%. Then, they were wrapped in plastic wrap and allowed to stand for 24 hours to ensure even moisture diffusion. After that, the soil mixture was weighed to determine dry density and placed in a three-axis sample preparation instrument. Finally, the soil was compacted into 5 layers, and surfaces between layers were shaved with a utility knife. To ensure the success of sample preparation, a static pressure sample preparation method was adopted.
2.3. Preparation of Artificial Acid Rain Solution
The components of acid rain in China and the concentration of each component were fully considered. In formulating artificial acid rain, sulfuric acid, nitric acid, and hydrochloric acid were mixed with a molar ratio of sulfate, nitrate, and chloride ions of 5 : 1 : 1. This solution was diluted to obtain a simulated acid rain solution with pH = 1. To avoid the impact of impurities present in tap water on the test results and ensure that the simulated acid rain represented the natural situation, deionized water was used as the solvent during the preparation process.
2.4. Test Methods and Principles
Soil samples with the same initial damage were selected and divided into 5 groups for testing. The test conditions were 0, 1, 2, 3, and 4 dry and wet cycles, respectively. The test process is shown in Figure 2. On the basis of multiple experiments, it is believed that the soil sample changed from a saturated state to a state with a natural level of moisture after air-drying for 2 days. The soil sample reached a saturation state in 1 day under pressure-saturation conditions.
The mechanical test was a three-way synchronous isobaric consolidation test without drainage (CU). The strain control module and the autocontrol mode were used to maintain a loading rate of 0.02 kPa/min. The net confining pressures of the preliminary consolidation step were 100, 200, 300, and 400 kPa. Subsequently, the triaxial shear test was carried out at a loading rate of 0.02 kPa/min with a constant confining pressure and increased axial pressure. The point at which the axial deformation reached 14.5% was taken as the end of the test, and the stress-strain curve was drawn. A sample with an initial consolidation pressure of 400 kPa was selected for the microscopic test and was first scanned by NMR to obtain the T2 spectrum, which was used to determine the corresponding pore characteristics according to the relevant formula. After the porosity tests, the samples were separated and tested by SEM to analyze the influence of the microstructure on the mechanical properties. The specific experimental process is shown in Figure 3.
3. Test Results and Analysis
3.1. Stress-Strain Curve
Figure 4 shows the soil stress-strain curves obtained from triaxial compression tests with different numbers of cycles. The figure shows that all the curves contained the characteristics of stress hardening. As the strain increased, the stress generally showed an initial growth trend. That is, the deviator stress gradually increased with increasing strain, but the growth rate gradually slowed. The content of clay minerals in loess is considered to be relatively high. The cohesion in the soil first played a role, and then the internal friction angle determined the stress, which is the main reason for this phenomenon. In addition, the soil deformation and initial confining pressure were correlated. When the initial confining pressure of the soil was large, the pores between the soil particles were compressed, making the soil structure compact and resulting in small deformations due to the external force. Finally, the curves are roughly divided into two stages: the elastic stage and the elastoplastic stage. The relationship between stress and strain gradually transitions from linear to nonlinear.
(a)
(b)
(c)
(d)
(e)
These experimental results show that as dry-wet cycling progressed, the peak intensity of the triaxial compression curve gradually decreased, indicating that acid rain infiltration caused structural damage to the sample, and the mechanical index was reduced. To clarify the effect of the number of cycles on the cohesion and internal friction angle, the test results were plotted in a rectangular coordinate system according to the Mohr–Coulomb criterion, and the shear strength shown in Table 2 was calculated. Table 2 shows that as the number of cycles increased, the cohesive force gradually decreased, and the amplitude of the change gradually increased. Moreover, there was a negative correlation between the internal friction angle and overall number of cycles, but when the number of cycles was 2, there was an abrupt change.
3.2. Porosity Test Based on Nuclear Magnetic Resonance Technology
Rock-soil samples are frequently used as materials for porous media, and the pore structure often directly determines the mechanical strength of the sample. Therefore, it is necessary to test the pores of the rock and soil samples. In this study, a low-field nuclear magnetic resonance spectroscopy system was used for testing. The system mainly determines the relaxation time of the pore fluid based on the relaxation mechanism for fluids in porous materials, as shown in the following equation:where is the free relaxation time of the liquid, which is determined by the physical properties of the liquid (such as viscosity and chemical composition); is the diffusion coefficient; is the magnetic field strength; is the magnetic spin ratio; is the surface relaxation strength (the relaxation strength of the particle surface); and is the ratio of pore surface area to fluid volume. For spherical pores, the ratio of surface area to volume is r/3.
In the calculation, the free relaxation time is longer, and the diffusion relaxation time is smaller, so the above formula is simplified to obtain the relationship between the pore structure and the relaxation time, as shown in the following equation:
3.2.1. T2 Spectrum Distribution Law
The test results obtained by nuclear magnetic resonance spectroscopy were converted to a logarithmic coordinate system, and the nuclear magnetic resonance T2 spectra obtained after different numbers of cycles are shown in Figure 5. This figure shows that the T2 spectrum of the test soil basically has three peaks; the main peak is 75%∼85% of the total area of the T2 spectrum and the vast majority of the T2 spectrum curve. The main peak of the T2 spectrum is in the range of 0.98–23.8 ms, which provides a considerable amount of information about the pore distribution of the soil. The small peak accounts for the smallest proportion of the total spectrum area, generally between 0.08% and 0.19%. During the test, the pore evolution of the test soil sample showed a strong positive correlation with the number of cycles. When the number of cycles was large, the T2 data were mostly located on the upper and right sides of the spectrum, which indicated that the effect of acid rain erosion on the development of pores in the specimen was more important. Similarly, the pore area and number of cycles showed a strong positive correlation.
There is a good correlation between the NMR T2 spectrum and soil porosity. Therefore, it is very reasonable to use the spectral area to reflect the total amount of pores in the test soil. Taking the T2 spectrum area obtained after different numbers of cycles as an example, the spectral areas and number of cycles are shown in Table 3. Table 3 shows that, for the same immersion time, the greater the number of cycles, the greater the area of the T2 spectrum; specifically, when the number of cycles is 4, the performance is particularly obvious. At this point, the rate of change of the spectral area is the largest, and the erosion effect dominated by is more obvious.
3.2.2. Changes in Soil Pores
Because the core principle of NMR technology is based on the correlation between T2 relaxation time and the sizes of soil pores, the relationship between the two is established. Therefore, from many investigations and comprehensive consideration of the structural characteristics of the test soil samples, the pores of the soil sample are believed to be mainly spherical. Then, a surface relaxation strength of 3.0 μm/ms was selected and substituted into the above formula to obtain the relationship between the radius and relaxation time. Equation (3) is as follows:
The form of the T2 spectrum distribution was considered comprehensively, the sample with the initial consolidation pressure of 400 kPa was selected for testing, the percentage of each pore volume was calculated, and the curve indicating pore diameter change with the number of cycles was drawn (see Figure 6). Figure 6 shows that the pores in the sample were concentrated in the range of 0.14–131 μm after different numbers of cycles. However, the diameters of the pores corresponding to the main peak were in the approximate range of 1.0–46.5 μm, and the pore diameter corresponding to the peak of the relaxation time curve was also in this range. In addition, the test results show that, with the increasing number of dry and wet cycles, the pore size distribution was relatively broad, and there were more pores with larger sizes.
To directly reflect the changes in pore diameter using the nuclear magnetic resonance test results, the relaxation time was converted into the pore diameter by using equation (3). The test results show that the 1H atom information was not captured in the groups with diameters greater than 3,000, so the test results were divided into 8 groups of [0, 1), [1, 2), [2, 5], (5, 20], (20, 50], (50, 200], (200, 1 000], (1 000, 3 000], and (3 000, 6 000), soil pore sizes were plotted, and the percentages of the pore size groups were used to construct a histogram, as shown in Figure 7. Figure 7 shows that the pore radii of the soils are mainly concentrated in four groups, [1, 2), (5, 20], (20, 50], and (50, 200). Among these pores, the largest number was found in the range of (5, 20). Figure 7 shows that as the number of cycles increased, the proportion of pores in the [1, 2) group gradually decreased, and the proportions of pores in the (20, 50] and (50, 200) groups gradually increased, which indicates that the dry-wet cycling process caused the pore sizes of the sample to transition from small- to large- and medium-size pores under acid rain conditions.
3.3. Microstructure Based on Scanning Electron Microscopy
3.3.1. Qualitative Analysis of Scanning Electron Microscopy Results
The test samples were magnified 300, 500, 1,000, and 2,000 times to obtain scanning electron microscopy images, as shown in Figure 8. Comparing the images with different magnifications showed that there were pores of different shapes and sizes between the particles of the test soil samples; the pore areas were slightly smaller after dry-wet cycling, and some overhead pores were present. In addition, the contacts formed between particles mainly involved indirect contact and surface-to-surface contact, with less side-surface contact and little point-to-point contact. Only a large external force acting on soil particles destroys the contacts between particles. Therefore, the structure of the soil sample was very stable.
The above figure shows that the changes brought by the dry and wet cycle process mainly involve the structure and pores. In the soil sample structure, the contact mode between particles changed, and some surface-to-surface contacts gradually transitioned to edge-to-surface contacts and finally to edge-to-edge or point-to-surface contacts; that is, the contact form gradually changed from stable to unstable. Pore sizes gradually transitioned from small pores to macropores, but at the same time, some large pores decreased in size. This was mainly manifested by the continuous decrease in the sizes of particles of clay minerals located among the soil particles; some larger agglomerates gradually separated and decreased in size, but some fine particles reaggregated into blocks. The results of the coupling action between the two behaviors during the migration process are as follows. First, acid rain solution erosion caused the dissolution of some mineral components and formed some overhead pores; the water remained in the soil. Along the water migration path, water repeatedly scoured the original pores, entrained some small soil debris, and dissolved some clay minerals, which caused the pore sizes to increase continuously. The dissolution and loss of soluble minerals caused some of the larger clay mineral aggregates to gradually separate and disperse; at the same time, the dissolution and loss of some minerals aggravated the formation of overhead pores between soil particles, a partial vacuum zone appeared, and the contact areas between particles gradually decreased.
3.3.2. Quantitative Analysis of Scanning Electron Microscopy Results
SEM is an important method for studying the microstructures of rock and soil, and SEM images contain rich information. Based on images, there are two ways to study the microstructure of rock and soil. One way involves directly exploring the combination and distribution of particles and pores, which is not comprehensive. The other way involves using modern computer vision theory to extract more information about rock and soil mass microstructure, which constitutes a new research direction. However, the second approach needs to be based on image segmentation. This is because the quality of image segmentation directly affects the accuracy of subsequent experimental results. This paper uses the Otsu algorithm to ensure that the particles and pores are segmented accurately in the SEM images.
This algorithm is consistent with the basic idea of the traditional threshold segmentation method and uses a threshold to divide the pixels in the image into two categories. The difference between the two methods is that the Otsu algorithm uses the variance of the image to automatically determine the segmentation threshold. Variance is a measure of the uniformity of the gray distribution of an image. The greater the variance between the background and the target area, the greater the difference between the two parts of the image. When part of the background is incorrectly classified into the target, or part of the target is incorrectly classified into the background, the difference between the two parts decreases. Therefore, the segmentation that maximizes the variance between the classes has the lowest probability of misclassification.
The Otsu algorithm is relatively simple, and because it is readily adapted to light, it can avoid the effects of contrast and brightness. The Otsu algorithm considers that the image contains L gray levels. When the gray value is , the number of pixels is and the probability and mean value of the pixel are
According to the grayscale characteristics of the image, the threshold is used to divide the image into two categories: target and background . When and represent the thresholds , the probability of occurrence of and is as follows:where represents the variance between classes with a threshold value of T in the histogram, defined as
The optimal threshold is defined as the T value corresponding to the maximum variance between classes, namely,
In this paper, the Otsu algorithm is written using a Python program. It is assumed that the image pixels can be divided into two parts according to the threshold: the background [background] and the target [objects]. Then, the optimal threshold is calculated to effectively distinguish the two types of pixels.
(1) Data Extraction and Analysis. The parameters for the microscopic morphologies of the extracted particles were selected in ImageJ software. Because of its length, this paper focuses on the analysis of the changes in soil particle abundance and fractal dimension as a result of the effects of dry and wet cycles.(a)Soil grain abundance The abundance C is the ratio of the short axis to the long axis of soil particles, ranging from 0 to 1, and it is used to measure how close the particles are to the plane. The calculation formula is where and are the lengths of the short and long axes of the soil particle, respectively. Figure 9 shows the variation of soil particle abundance with the number of dry-wet cycles. The figure shows that the abundance values of the test samples are mostly in the range of 0.4–0.9, and there is no particle distribution in the range of 0.1, so it is considered that the roundness of the test soil sample is better. Overall, the initial particle content was low for the abundance value 0.2, but some particles with this abundance value appeared after the dry-wet cycles. Therefore, water migration divides some larger particles in the path to form smaller particles during the dry and wet cycles. In addition, although the curve shows that the variation of soil particle abundance values is relatively complex, the trend shows that with the number of dry-wet cycles, the soil particle abundance gradually increases, which also verifies that the dry-wet cycles constantly make the particles rounder, and the long strip-like particles constantly transform into flat and spherical-like particles; this is the mechanism for changes in grain size. However, when water divided the large particles, the small particles continued to adhere to form new particles, which is another reason why the abundance change is complex.(b)Fractal dimension of the shape distribution The fractal dimension is mainly used to measure the roughness of the soil particle unit using the box counting method and equivalent area perimeter, and it is calculated according to the following formula: where is the equivalent perimeter of the particle, is the equivalent area of the soil particle, refers to the fractal dimension of the soil, and is the fitting constant.
Figure 10 shows the variation curve of the fractal dimension of soil particles with number of dry-wet cycles. In the process of drying and wetting, water continuously washes and dissolves various minerals, which makes the original soil particles decompose and form particles of different sizes, resulting in the overall decreasing trend of the fractal dimension, especially after one drying and wetting cycle. However, the test results show that the fractal dimension of the test soil sample is larger as a whole. If a sample does not undergo the dry-wet cycle, the fractal dimension is 2.68, and the particle distribution is more uniform.
4. Discussion of Erosion Mechanism
The spatial structure characteristics of the rock and soil body are called the scale effect. Using the ranges >10−3 m, 10−6～10−3 m, and 10−9～10−6 m as standards, the damage to rock and soil can be divided into three types: macrodamage, mesodamage, and microdamage. One of the important directions in the development of geotechnical engineering involves comparing test results from multiple scales to reveal the mechanism of the experiment. As the main method for microscopic level detection, scanning electron microscopy can reveal images of soil under different magnifications. In the mesoscopic range, the T2 spectrum obtained by nuclear magnetic resonance spectroscopy indirectly reflects the changes in the pores of rocks and soil. The macromechanical properties are determined by the combined fine and microscopic soil structural characteristics. The test results obtained at the microscopic level showed that there were pores of different shapes and sizes between the reshaped loess grains. However, the pore area between the grains was relatively small, and there were overhead pores between the grains. In addition, the contact between grains mainly comprised indirect contact and surface-to-surface contact, with less side-to-surface contact and almost no point-to-point contact. At the mesoscopic scale, the distribution of the T2 spectrum showed three peaks, in which the area of the main peak accounted for the vast majority of the total area. Moreover, the result of curve inversion was good, and the pore distribution was relatively uniform.
As the dry-wet cycles progressed, the contacts between particles gradually transitioned from surface-to-surface contacts to edge-to-surface contacts and finally to edge-to-edge or point-to-surface contacts. After the sample was subjected to continuous water migration, the proportion of small pores gradually decreased, the proportions of large and medium pores gradually increased, the particle abundance value continued to increase, and the fractal dimension gradually decreased. At the mesoscale, the area of the main peak of the T2 spectrum still accounted for the vast majority of the total. As the number of cycles increased, the peak of the spectrum roughly increased, the spectral data continued to shift to the right, and the spectrum area gradually increased, showing a good correlation with the number of cycles. In the process of a dry-wet cycle, the internal structure and pore state of soil are changed by moistening, dissolving, scouring, and carrying processes. First, acid rain fills the pores of the sample and dissolves and erodes clay minerals such as agglomerates in the sample, which is one of the reasons for the change in soil cohesion. Second, under the joint action of dissolution and erosion, the contact form between soil particles changes to a certain extent, which is another main reason for the change in soil cohesion. Specifically, the decrease in clot area makes them directly desorb. The contact form of indirect contact is invalid. Thus, the stable contact form of the original surface-to-surface contacts develops into the unstable contact form of edge-to-edge contacts or point-to-surface to surface contacts, which weakens the interactions between molecules. The dissolution and erosion reactions interact and promote each other, and both become increasingly intense with the increase in the number of dry-wet cycles. In addition, for friction, after the dry-wet cycles, the protruding parts of soil particles are gradually smoothed, the particle shapes tend to be rounder, the surfaces are smoother, and the concentration of particle distribution becomes worse, which makes the sliding friction between particles and the soil erosion worse, resulting in the decrease of the internal friction angle with increasing number of cycles.
At the same time, it should be pointed out that the wet-dry process did not simply involve physical erosion, since this was accompanied by chemical reactions between calcite minerals and feldspar minerals under acid rain erosion conditions. Calcite contains a large amount of calcium carbonate, which reacts with to form overhead pores. At the same time, feldspar (including albite and potash feldspar) also undergoes chemical reactions under acid rain erosion conditions, which erodes some particles and increases the porosity of the sample. The two chemical reaction mechanisms are as follows:
The combined actions of physical erosion and chemical erosion caused an increase in the pore structure, a partial instability in the interactions between grains, and a gradual decrease in intergranular embedding and occlusal interactions. The deterioration caused by simultaneous sliding and occlusal friction reduced the mechanical strength of the soil as the number of cycles increased.
5. Conclusion
In this study, we researched remolded loess and investigated the influence of alternating dry and wet cycles on its mechanical properties and microstructural characteristics in an acid rain environment. The experimental results may be summarized as follows:(1)In an acidic environment, the deviatoric stress of remolded loess increases with increasing strain after dry-wet cycles, and the stress-strain relationship shows evidence of stress hardening. Under the double strength effect, cohesion is primary, followed by the internal friction angle, but the cohesion and internal friction angle are negatively correlated with the number of cycles.(2)Because the Otsu algorithm has a strong adaptive ability for light, it can avoid the influences of contrast and brightness, and it can be applied to the processing of soil SEM images.(3)The T2 spectrum of remolded loess basically has three peaks, the main peak accounts for 75%∼85% of the total area, and it has a great influence on the pore distribution of soil. When the number of dry-wet cycles increases, the spectrum area of the mesoscale level increases gradually; that is, the sample transitions from small pores to large- and medium-size pores. At the microscopic level, clay minerals and agglomerates decrease. This results in the weakening of the contact mode between soil particles and the cutting and rounding of large particles. The fractal dimension of the original soil particles is reduced due to continuous decomposition and grinding.(4)Compared with dry-wet cycles using neutral water, the dry-wet cycling test with acid rain has both certain commonalities and special characteristics. As long as the dry-wet cycle is a physical process, pore expansion, contact form weakening, and particle shape change occur. However, acid rain is often accompanied by certain chemical changes; that is, calcium carbonate in calcite and feldspar form overhead pores, which accelerates the deterioration of soil structure. The results show that the macromechanical strength of the samples decreases with increasing number of cycles due to the combined action of physical and chemical erosion.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors thank the workers, foremen, and safety coordinators of the main contractors for their participation. This work was supported by the National Natural Science Foundation of China (Grant nos. 42072319 and 41902299).