Mechanism and Sensitivity Analysis of Collapse in Large Section Mountain Neighborhood Tunnels

Collapse is one of the most frequent geological disasters in mountain tunnel engineering, posing major safety concerns for underground structures and construction crews. According to a catastrophic collapse in the shallow buried area of the Huangjiakuang tunnel in Weihai City, Shandong Province, the contributing factors, that is, the surrounding rock property, influence of neighboring tunnel construction, and tunnel over-excavation are systematically investigated. The tunnel collapse, breaks through the ground surface, is inverted using MIDAS GTS NX. A tunnel deformation analysis model is created using the quantitative methods of grey relation analysis and entropy weight methods based on crown settlement and lateral wall horizontal contraction during tunnel excavation. The surrounding rock property, the distance between the left and right tunnel faces, and the over-excavation height have a significant effect on tunnel deformation, which are quantitatively interpreted using the tunnel deformation analysis model. This study is valuable for the multi-factor analysis of tunnel deformation and determining the main contributing factors to the collapse quantitatively.


INTRODUCTION
The continual development of China's transportation industry has resulted in more flawless infrastructure construction. Tunnels are becoming an increasingly important component of transportation networks, both in terms of number and total length. According to data from the People's Republic of China's Ministry of Transport, China added 14,032 new motorway tunnels between 2010 and 2020, totaling more than 169,000 km. Meanwhile, tunnel accidents are becoming more frequent. There were 97 tunnel accidents between 2002 and 2018, with tunnel collapse accounting for 62.89 percent of mountain tunnel geological disasters, and 150 workers or management losing their lives . The collapse accident prompted the project's extension, resulting in countless property losses and casualties, making it one of the most damaging disasters in tunnel engineering history.
Consequently, many scholars have studied the failure mechanism of the tunnel surrounding rock in detail, and the main research directions are theoretical analysis, model experiment, and numerical simulation. Among the theoretical methods, the classical pressure theory was the first to be used in the research on the surrounding rock's stability. Based on the Hoek-Brown failure criterion, Fraldi and Guarracino, (2009) ;Fraldi and Guarracino (2010) obtained the exact solution of tunnel collapse prediction in the field of plastic theory by using the classical tools of the variational method. Yang and Huang (2011) used the nonlinear Hoek-Brown failure criterion to carry out limit analysis on the collapse of the shallow circular tunnel under support pressure; obtained the collapse mechanism of the shallow circular tunnel. By utilizing the nonlinear Hoek-Brown failure criterion, Zhang et al. (2014) deduced the analytical solution of the collapse block shape curve by replacing the elementary function with the total potential energy in the form of functional and introduced the criterion to judge whether the roof collapse occurred in the deep tunnel. Li and Yang (2017) modified Fraldi and Guarracino (2009) ;Fraldi and Guarracino (2010) tunnel collapse solution by introducing a simplified method based on the variable separation velocity of the yield surface. Ding et al. (2019) proposed a new three-dimensional collapse mechanism, which shows that the working face stability of a circular tunnel in pure clay soil is related to the cover diameter ratio and provides a scheme for soft soil tunnel collapsing stability. By utilizing the upper bound theorem of limit analysis and the nonlinear Hoek-Brown yield criterion, Lyu and Zeng (2019) obtained the range of potential collapse blocks and the total gravity equations of rectangular and circular tunnels in shallow tunnels under the inclined strata using the variational principle. These theoretical investigations serve as a theoretical foundation for the application of other research approaches.
Concurrently, many articles studied tunnel collapse by applying the modeling experiment (Barla, 2008;Yang et al., 2019;Liu et al., 2020;Hao et al., 2021;Wang et al., 2021;Lan et al., 2022;Li et al., 2022). The simulation experiment allows for a more intuitive understanding of the structure's failure process and a detailed investigation of its failure mechanism. Li et al. (2014) developed an artificial speckle to monitor the change in the strain field of the specimen, and based on the developed degradation constitutive model, the failure mechanism and deformation law of tunnel surrounding rock are studied. Cui et al. (2021) applied monotonic and cyclic loading tests in the laboratory and determined that the carbonate fault breccia on the sliding layer has a high liquefaction potential. Lei et al. (2015) systematically studied the variation and distribution of structural stress and surrounding the rock stress and the failure mechanism of the lining. Zhang et al. (2016) utilized a variety of methods to measure the convergence of the large deformation, the displacement, and damage mode of the critical points of the sidewall and obtained the change process and regular curve of the displacement of the surrounding rock and the critical issues of the sidewall in the excavation stage. Xiao et al. (2016) conducted a geomechanical model test to investigate the deformation, stress, and potential damage of the Baihetan hydropower tunnel. After the triaxial compression test, Xia et al. (2022) established standard specimens by 3D printing technology and obtained the microscopic failure mechanism of the columnar jointed rock mass. Zhang et al. (2020) quantitatively analyzed the peak strength, effective elastic modulus, and related failure characteristics of fractured marble using the discrete element method and numerical experiments. Liu et al. (2019), Liu et al. (2022) studied the progressive failure process of surrounding rock using field microseismic (MS) monitoring technology and provided a method for predicting collapse. Zhou et al. (2021a) studied the stability of the crown pillar from induced microseismicity and that it is helpful for the characterization of seepage channel development. Sx et al. (2020) established the tunnel model of gully terrain, analyzed the relationship between tunnel construction collapse and longitudinal terrain, and provided theoretical support for constructing shallow tunnels under gully terrain. However, model testing takes a significant amount of time and money. This strategy frequently involves a significant amount of time and effort to achieve the desired results. It is hard to conclude universally applicable conclusions.
Over the past few decades, global technological levels have multiplied. Numerical simulation methods, such as the discrete element method (DEM), finite element method (FEM), and DEM-FEM, are widely used to study the failure of the tunnel surrounding rock. The variation laws of stress and strain for these structures are obtained during the entire process of inversion construction (Jia and Tang, 2008;Kong et al., 2016;Li et al., 2020a). Zhang et al. (2011) calibrated and predicted the excavation damage zone of the tunnel section using continuous fast Lagrangian analysis (FLAC) and particle flow code (PFC). Huang et al. (2013) applied experiments and numerical simulations to study the influence of weak interlayers on the failure mode of the tunnel surrounding rock. Zampieri et al. (2018) used ABAQUS software to conduct quantitative analysis combined with construction conditions, obtained the deformation law of the tunnel without support, and confirmed the causes of motorway tunnel collapse. Li et al. (2019) used finite element software to study the stress, strain, and potential damage area of tunnel construction in hilly terrain with slope stratification. Yan et al. (2021) used FLUNT software to study the collapse accident caused by water inrush in the Qingdao subway and proposed remedial measures using advanced pipe roof grouting reinforcement. These studies help to understand the mechanism of tunnel collapse and improve the stability of tunnel construction. Nonetheless, it is impossible to evaluate the causes of tunnel collapse caused by multiple factors. The influence of different factors on tunnel collapse must be accurately interpreted, and the leading causes of tunnel collapse cannot be established.
Meanwhile, various statistical methods or others have been used (Ma et al., 2020;He and Kusiak, 2017;Li et al., 2021a;Li et al., 2021b). In these tools, statistical methods account for the majority. In disaster prevention and mitigation of geotechnical engineering, grey relation analysis can determine the relation degree between disaster factors according to the similarity of change curves, so this method has been widely used (Deng, 2014). It can quantitatively analyze the development of dynamic processes, complete the comparison of geometric relations of time series in the system, and attain the grey relation degree between the reference sequence and the comparison sequence. Zhang et al. (2019) used this method to evaluate the relationship between soil elastic modulus, stress release rate, and ground settlement during tunnel construction and established an equation for predicting ground settlement. Qiu et al. (2021) modified this theory and determined the relation degree of the water inrush risk level by taking fault zone properties, fault structure characteristics, and other factors as criteria. Based on 24 collapse accidents, Kim et al. (2022) proposed a tunnel collapse assessment model based on a quantitative classification system and significant influencing factors of tunnel collapse risk. These research studies give the impact of different factors on accidents from a statistical perspective but do not consider the weight of various factors on landslides. Alternatively, judging weights by an expert review is not objective enough.
The entropy weight method is an objective weighting method. The information entropy is used to estimate the entropy weight of each index according to the variation degree. Then, the weight of each index is corrected by entropy weight. The entropy weight method overcomes the intense subjectivity of the expert evaluation method and equal distribution method and determines the index weight according to the measured data. A more reasonable and reliable index weight can be attained (Zhou et al., 2021b). Thus, the combination of grey relation analysis and the entropy weight method can more effectively determine influencing factors. Currently, multiple pieces of literature have studied this problem (Zhang and Yang, 2018;Li et al., 2020b;Zheng et al., 2021). There are many studies on tunnel collapse. Nonetheless, in many studies, some scholars attribute the cause of tunnel collapse to a particular factor. Alternatively, many factors have been considered but there are fewer related studies (Wang, 2021). The method to determine the leading cause of the collapse was not introduced in the multifactor accident. The research system to analyze the failure under multi-factor was not established. Therefore, it is essential to establish an effective collapse analysis method for engineering accident investigation. This article studied the collapse accident which occurred at a large section of the mountain neighborhood tunnel's shallow buried area in Weihai City, Shandong Province. The properties of the surrounding rock, the distance between the left and right tunnel faces, and the over-excavation height are analyzed in detail, and the collapse accident was inversed with the MIDAS GTS NX. On this basis, the three factors are analyzed by controlling variables. Based on the crown deformation and horizontal convergence, the grey relation analysis-entropy weight method is used to calculate the relation degree. The main factors influencing tunnel deformation among the three factors are determined. This research proposes a new analysis approach for the investigation of multi-factor collapse events by performing a quantitative analysis of numerous components impacting tunnel deformation. This study is valuable for the multi-factor analysis of tunnel deformation and for determining the main contributing factors to the collapse quantitatively.

Project Overview
As illustrated in Figure 1A, the engineering location is the Huangjiakuang tunnels in Weihai City, Shandong Province, China. The tunnels are located southwest of the village of Huangjiakuang. These tunnels are twin-lined six-lane roadways with a design speed of 60 km per hour, divided into left and right tunnels. The left line tunnel has an import and export mileage of ZK5 + 540 and ZK6 + 000, whereas the right line tunnel has an import and export mileage of YK5 + 533.5 and YK6 + 113.5. The tunnels span 580 m in length and reach a maximum depth of 47 m. The project is a large section neighborhood tunnel in the mountain region, measuring 134.89 m 2 in the section area. The tunnels have a maximum width of 15.5 m, while the middle rock pillar has a width of 10.09 m. The whole tunnels' terrain is steep in the north and flat in the south, with a slope of 20°-45°. YK5 + 540 to YK6 + 120 is depicted in Figure 1B. According to the survey study, gneiss covers the tunnels in large part. The buried depth varies significantly, from 4 to 45 m. The foundation soil of the Huangjiakuang tunnels is separated into three strata from top to bottom: plain fill soil, strongly weathered granitic gneiss, and weakly weathered granitic gneiss. China divides the strength of the surrounding rock of the highway tunnel into six grades according to the rock mass basic quality (BQ). It is a comprehensive indicator which includes three aspects. The surrounding rock classification basis is obtained according to Eq. 1. The strength of the surrounding rock decreases from Ⅰ to Ⅵ.
where BQ is the rock mass basic quality, R c is the uniaxial compression strength, andK v is the intactness index of rock mass.
According to the geological description, the tunnel ground is primarily composed of three types of soil. At the shallow buried section of the tunnel, the rock mass basic quality is 207 and its uniaxial compression strength is 20 MPa. On the basis of the classification standard of the surrounding rock in China, this kind of geology is judged as grade Ⅴ surrounding rock.

Tunnel Collapse
On 2 October 2018, the right tunnel to the area denoted in Figure 1B was constructed. The working face (YK5 + 680) was abnormal. Soil and stones began to fall above the left area of the tunnel. Many soils and fragmentary rocks started to fall, ultimately becoming a collapse accident directly to the surface. A ground settlement area of 7 m long and 5 m in depth was formed after the collapse, and both tunnels were compelled to halt. Figure 2 shows the photograph of the collapsed roof fall wreck. It can be seen from the pictures that the soil in the collapse area is predominantly soil. After investigation, the engineering data show that the collapsed roof area is grade V surrounding rock. The tunnel construction method is the bench-cut method, and the cyclic footage is 4 m. Only the upper step which is 16 m in length was carried out. The soil above the tunnel face collapsed after constructing the YK5 + 678 area. So, the size of the excavation section expanded. The highest collapsed soil was 1.8 m above the tunnel vault contour.

Analysis of the Disaster Causes
After the collapse accident, the tunnel construction was compelled to stop, resulting in delays in the construction period and economic and property losses. The construction conditions and the mechanical parameters of the soil at the collapse position are collected. The collapse occurred on October 2 nd , and there was no substantial rainfall at the Huangjiakuang tunnels before the disaster. In September, the rainfall at this site is 36.0 mm, while the annual rainfall is 596.8 mm. It only makes up a small fraction of the total. According to the survey report, the shallow buried section has minimal groundwater resources. The water outlet formed during the excavation process is primarily linear and drip seepage. As a result, the reasons for hydrological and climatic conditions are not thoroughly investigated. It primarily investigates the nature of the surrounding rock, excavation events, and structural aspects of tunnel building.

Surrounding Rock Characteristics
The collapse location of tunnels pertains to the shallow buried section (YK5 + 648.5-703.5). The surrounding rock of grade V covers this part of the right tunnel. The geological exploration report shows that tunnels are mainly composed of broken soil and gravel containing sand and gravel. Tunnel collapse (YK5 + 680) is  Rock classification Ⅲ ≤0 Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 904655 6 strongly weathered by granitic gneiss, rock fragmentation, and joint development. The rock mass basic quality (BQ) is 207 and its uniaxial compression strength (Rc) is 20 MPa. The tunnel face condition is shown in Figure 3A.
Ground-penetrating radar (GPR) is one of the most popular electromagnetic methods in tunnel ahead prospecting, which can propose information about the subsurface with high resolution . Before the collapse, the ground-penetrating radar was used to detect the excavation face. The use frequency of ground penetrating radar is 100 MHz, and the maximum depth that can be identified is 30 m. Two lines are laid on the working face, line 1 from left to right and line 2 from right to left. Radar image processing is shown in Figure 3B. YK5 + 678-YK5 + 708 regional geological radar image is shown as Figure 3B. The geological report shows this place is grade V surrounding rock. The starting point of the radar image is 2 m before the collapse position, and the ending point is 28 m after the accident location. The collapsed position is in the red area and is shown in Figure 3B. It can be seen that the pulse signal in the red region is predominantly medium-low frequency from Figure 3B. The waveform is disordered and discontinuous, and the phase axis is cracked. The amplitude has multiple oscillations, and the pulse frequency is different and irregular. The reflection pulse of the broken rock is distributed in a band, and the reflection intensity of medium-and low-frequency signals is high. The geological radar image of YK5 + 843.5-YK5 + 873.5 area is also illustrated in Figure 3B. The geological report demonstrates that grade III surrounding rock covered this area. Figure 3B shows that the blue area 2 m in front of the radar is grade III surrounding rock. The reflected pulse signal in this area is predominantly medium and high frequency, which changes slightly. The in-phase axis is steady and uniform, and the signal amplitude is small. The properties of the surrounding rock are better than those shown in Figure 3B.
To see more intuitively the influence of different surrounding rocks on tunnel deformation, the crown settlement monitoring data of grade III (BQ = 371) and V (BQ = 207) surrounding rocks are collected, as shown in Figure 3C. The tunnel crown settlement can be divided into three stages according to the settlement curve (Gao et al., 2018). In the first stage, the crown settlement of the grade V surrounding rock is faster than that of the grade III surrounding rock. The settlement of the grade V surrounding rock rises swiftly to 10 mm, and the grade III surrounding rock is just 6 mm. The second stage takes 40 days for the third-level surrounding rock to enter a stable state and 60 days for the grade Ⅴ surrounding rock to reach a steady condition. Currently, the settlements are 13 and 7 mm, respectively. In the third stage, the crown settlement of the three-level surrounding rock and five-level surrounding rock attained a stable state. Currently, the settlement of grade V surrounding rock was steady at 8 mm, and that for grade V surrounding rock was 14 mm. The maximum settlement and times to reach a stable state for grade V surrounding rock are greater than those of the grade III.

Excavation Disturbance of the Left Tunnel
The tunnels' section area is 134.89 m 2 , and the distance between the left and right tunnels is 10.09 m. The clear distance between the left and right tunnel is less than the tunnel's diameter. The distance between the two tunnels is shown in Figure 4. The construction process involves excavating the left tunnel first; then, the right tunnel is excavated. The two tunnels are constructed using the bench-cut method. The original equilibrium state in the soil is lost as a result of excavation disturbance. The stress in the rock mass is redistributed until a new equilibrium is found. The excavation disturbance has induced a change in the size, direction, and character of the original stress in the rock. The tunnel is excavated by blasting, which causes damage accumulation to the shared rock in the neighborhood tunnel (Cao et al., 2020). Therefore, the disturbance to the right tunnel after constructing the other tunnel is considered in this study. The underground cavern excavation will lead to the redistribution of soil stress in the range of 6R 0 . The R 0 is the diameter of a single circular tunnel excavation. (Han et al., 2020). In 'Technical Guidelines for Construction of Highway Tunnels' (U.S. Department of Transportation Federal Highway Administration, 2009), the mutual influence degree of two neighborhood tunnels is divided, as summarized in Table 1. In the Table, the parameter b is the width of the excavation section of a single tunnel. Under the condition of grade V surrounding rock, the thickness of middle rock pillar ≤0.75 b is seriously affected. The rock width in the Huangjiakuang tunnel is 10.09 m, and the tunnel span is 15.5 m. The tunnel's height is 9.77 m. So, the tunnels belong to the category of serious impact. The location of the right tunnel is in the excavation influence area of the left tunnel. Considering that the tunnel section is a noncircular tunnel, the equivalent radius of the tunnel section is calculated using the equivalent method (Peng et al., 2010). The equation of the equivalent method is.
where b is the tunnel span, h is the tunnel height, and the equivalent radius R 0 = 7.94 m is calculated. In this study, b is 15.5 m and h is 9.77 m, according to the equation: where σ r is the minor principal stress, σ 0 is the field stress, and γ is the distance of the left tunnel edge.
where σ θ is the maximum principal stress.
Only the variation of the maximum principal stress is discussed here. Thus, the influence area of the right tunnel is 1.10-1.20 σ 0 .

Tunnel Over-Excavation
The tunnel collapse occurs in the shallow buried area with the surrounding rocks of grade V. The upper and lower steps method is used for excavation in these areas. After excavation, the soil is hauled out of the tunnel by machinery, and then the anchor rod and shotcrete are set up. The primary support is C25 concrete. The C25 concrete is defined as concrete with a compressive strength of at least 25 MPa after 28 days of curing in laboratory settings. The thickness of the primary support is 30 mm. Figure 5A shows the construction diagram of the tunnel bench-cut method, and only the upper bench excavation is carried out during the collapse.
The right line of the tunnel was constructed to the YK5 + 677 area after the workers carried out the blasting operation, the rock fell over the right side of the tunnel, and the maximum rock particle size was 0.5 m, forming an over-excavation accident. The over-excavation position is located at the top right of the tunnel contour line vault, and the maximum height of over-excavation is about 1.8 m from the tunnel contour line. After over-excavation, the cross-section area of the tunnel expanded, and the cavity was formed above the right of the tunnel. The over-excavation area diagram is shown in Figure 5B. The over-excavation will change the stresses around the tunnel face and affect the deformation of the ground surface (Xu et al., 2021).  The aforementioned three factors are obtained from the investigation of construction conditions and the search of relevant geological survey reports. Whether it is related to the collapse requires further verification. Consequently, the FEM is used to verify our conjecture.

Inversion of the Tunnel Collapse
In this study, based on the FEM, the construction simulation of tunnel excavation is used to explain the deformation of the tunnels. Based on Figures 2, 4, a three-dimensional finite element simulation model of the tunnels in the collapse section is established, as shown in Figure 6. The tunnel size in the model is identified as the actual tunnel measurement. The model's dimensions in X and Y directions are 100 and 210 m, respectively, and the distance from the roof to the bottom is 50 m. The width of the tunnel is 15.5 m, and the height is 9.77 m. The clear distance between the two tunnels is 10.09 m. In the model, the left and right boundary constraints X and Y direction displacement, the bottom constraints X-, Y-, and Z-direction displacement. The building approach of tunnels is the benchcut method, divided into two steps. First, 16 m long upper steps were excavated, and the primary support construction was implemented next. The surrounding rock adopts a threedimensional solid element, and the effect of steel arch and steel mesh is converted into primary support using the equivalent stiffness method. The primary support adopts a two-dimensional plate element with a 0.30 m thickness. In actual tunnel construction, the secondary lining support was not carried out when the right tunnel collapsed. So, the secondary lining support was not considered in the tunnel simulation process, just the primary support was analyzed. The specific parameters in the model are summarized in Table 2.
The Z-direction displacement diagrams of the right tunnel are extracted before and after the collapse, which are shown in Figure 7. At the collapse position, the displacement diagram is different before and after the failure. In Figure 7A, the largest displacement occurs at the tunnel portal, with a maximum displacement of 11mm; in Figure 7B, it occurs at the collapse location, with a maximum displacement of 117 mm. As shown in Figure 7, the settlement at the arch crown of the right tunnel far exceeds the specified value, so it can be determined that the collapse accident has occurred.
Considering the major principal stress variation calculated in Equation 3 earlier, the major principal stress from the finite element simulation model is extracted after excavating the left tunnel. All results are shown in Figure 8A. The collapsed section of the tunnel is at the foot of the slope connecting two peaks. There is a topographic feature of left high and right low in the direction of the tunnel section, so the overall curve is left high and right low. The purple area is the right tunnel place. The maximum major principal stress in the right tunnel area is 0.33MPa, and the minimum is 0.25 MPa when the left tunnel is unexcavated. After the excavation of the left tunnel, the major principal stress of the right tunnel increased. The major principal stress of the right tunnel increases by 0.06 MPa which accounts for 20% of field stress, and the stress increase in the right edge was small, which was 0.023 MPa which accounts for 10% of the field stress. The calculation results are 1.1-1.2 σ 0 , and the numerical simulation results are consistent with the calculation results of Eq. 3. Therefore, the impact needs to be considered after the excavation of the left tunnel.
The crown settlement and tunnel horizontal convergence results are extracted from YK5 + 680 monitoring points, as shown in Figure 8B. When the working face is excavated before YK5 + 676, the deformation at the collapse position is slight. The small-scale deformation occurs in the monitoring point area when the working face is excavated to YK5 + 678. When the working face is excavated to YK5 + 680, a large settlement occurs, and the settlement abruptly changes by 8.2 cm, far surpassing the required 3 cm. At the same time, the tunnel diameter convergence has tremendous changes, which also exceeds the prescribed limit value.
Most existing displacement monitoring is based on singledimensional displacement parameters as instability criteria. Still, the apparent prediction parameters and criterion system of onedimensional displacement cannot fully and accurately describe the mechanism and law of evolution of deformation in the process of tunnel excavation. Therefore, the major principal stress variation curve of the left and right tunnels extracted from the landslide location mileage is shown in Figure 8C. With the excavation of the tunnel, the major principal stress increases. When the right tunnel is excavated to YK5 + 678, the construction of workers leads to over-excavation. The curve in the diagram also shows that the stress has a downward mutation, and the stress release occurs (Xu et al., 2021). The stress is changed from −0.02 MPa to -0.06 MPa. When excavated to YK5 + 680, the crown settlement and horizontal convergence caused mutation, and the mutation value is 9 cm and 8 cm. Also, the major principal stress changed from −0.06 to 0.04 MPa, and the tunnel surrounding rock stress changed from compression to tension. When the displacement increases from 2 to 11 mm, the surrounding rock of the tunnel produces large deformation. The left tunnel is still under pressure compared with the YK5 + 680 position. After excavation, the maximum principal stress decrease of ZK5 + 680 is from −0.04 to 0.1 Mpa and reaches a stable state after the excavation of ZK5 + 688.
In the essence, the whole process of collapse is caused by the coupling of three factors. When the tunnel is constructed to the collapse area, the soil property of the collapse location is poor. The previous excavation of the left tunnel impacts the collapse area, making the right tunnel's stress increase. After the workers' construction of over-excavation, the collapse accident is induced. The accident is directly connected to the ground surface, resulting in surface collapse. Nevertheless, the influence of the three factors on the collapse has not been confirmed, and further analysis is required to find out the main reason for the collapse accident.

Analytical Method
The experimental scheme was designed with the distance of two tunnel faces, the surrounding rock strength of the collapse section, and different over-excavation heights. To gain a better understanding of the critical factors influencing the deformation of the surrounding rock, the impact of three parameters on the collapse accident is assessed. This section will quantify the variables that contributed to the collapse accident. In this study, the influence of overburden thickness is excluded. The main factors affecting the stability of the surrounding rock are determined by combining the grey relation analysis and the entropy weight method. Considering the different judgment bases of surrounding rock properties, this study only changes the elastic modulus of soil according to the requirements in the specification. Grey relation analysis has low requirements for working conditions, so nine schemes are set up, and the test scheme is shown in Table 3. The key factors affecting the stability of tunnel surrounding rock are determined by the tunnel crown settlement and horizontal convergence as grey relation analysis evaluation indexes.

Grey Relation Analysis Method
The grey relation analysis method is applied to determine the effect of multiple factors to the collapse. By determining the geometric shape similarity between the reference data column and several comparison data columns, the tightness of their connection is judged. Also, the relation degree between their comparison sequence and the reference sequence is calculated. Then, the results are sorted according to size. The specific steps of the method are as follows.
First, the grey relation analysis method is used to calculate the relation coefficient. Classify the sorted data and divide the corresponding matrix. The reference matrix is shown in Eq. 4: where X 0 is the reference matrix. Then, list the data of the reference matrix. The comparison matrix is shown in Eq. 5: where Y i is the comparison data matrix. Because the comparison data matrix and the reference data matrix units are different, to determine the relation coefficient, the relation coefficient between the data is usually used to unify the data, so that each sequence is dimensionless. The normalization method is as follows: Therefore, considering the influence of the comparison sequence on X 0 for reference sequence X i , the relation coefficient can be determined as follows: where Δ max is the maximum difference of level 2, Δ min is the minimum difference of level 2, and Δ ij is the absolute difference. In this study, ρ is the resolution coefficient. Usually, ρ is chosen as 0.5; and the relation coefficient matrix is obtained after Eq. 7 is processed.
As indicated in Equation 8, the relation coefficients of each parameter are averaged to yield the relation degree of each parameter.

Entropy Weight Method
The entropy weighting method is a form of objective weighting. It calculates the index's entropy weight based on the change in each index. After that, the index is corrected using the entropy weight method. In comparison to the original grey relation analysis, the entropy weight approach can confirm the weighting effect of many indices and improve the relation coefficient's accuracy. As a result, the matrix is further processed using the entropy weight method.
Since the measuring units for various indicators are not standardized, the data must first be normalized. That is, transform the index's absolute value to its relative value. The standardizing equation is as follows: where the p ij is the normalized data. Then, using the standardized data, calculate the index's entropy. The Entropy value is critical for the subsequent weight calculation.
where e j is the entropy value. Next, the weight of each parameter is determined using the entropy calculated previously.
where w j is the weight of each parameter. Finally, the grey relation analysis and entropy methods are combined to yield the final modified grey relation coefficient matrix.
where ψ is the modified grey relation coefficient matrix.
The main factors affecting the tunnel deformation are calculated by comparing the grey relation degrees of overexcavation height, surrounding rock property, and the distance between the left and right tunnel faces.

Result Analysis
Taking the over-excavation height, the distance of the left and right tunnel faces, and the elastic modulus of the surrounding rock as control variables, the tunnel constructions are simulated using the finite element method. Extract the tunnel crown settlement and convergence results, and the roof deformation is positive upward, negative downward, horizontal convergence using the absolute value judgment. All the results of tunnel deformation under different variables are shown in Figure 9. Taking the surrounding rock strength as the control variable, the specific curve is shown in Figures 9A,B. When the tunnel face is progressively advanced, the crown settlement and convergence at YK5 + 680 is increasing, ultimately attaining a stable state in these cases. The maximum crown settlement increases from 0 to 3.0 cm, and the horizontal convergence rises from 0 to 1.25 cm. Before the tunnel is excavated to YK5 + 676, the deformation at YK5 + 680 is almost 0 cm. The crown settlement increases by 2.5 cm after excavating the soils of YK5+680. The horizontal convergence occurs the most significant deformation after YK5 + 680 excavation. The deformation amplitude reaches 1.5 cm in three cases. It can be seen that when the over-excavation height is 1.8 m, the crown settlement and convergence reach the maximum, and the maximum obtained is 3.0 and 1.2 cm. The deformation gap is not evident under the three conditions. The settlement of 1.8 m over-excavation height increases by 0.4 cm, and the convergence increases by 0.2 cm compared with 0.6 m. Therefore, the over-excavation height affects the tunnel  deformation, but the tunnel deformation is not apparent with the rise of over-excavation height. Taking the surrounding rock as the control variable, the specific curve is shown in Figures 9C,D. The settlement of the tunnel crown under different conditions is from 0.6 to 1.5 cm, and the convergence is from 0.3 to 0.75 cm. The deformation of the surrounding rock under 0.2 GPa strength is the greatest when compared to 2 and 1 GPa surrounding rock. The crown settlement reaches 1.5 cm, and the tunnel convergence is 0.75 cm when the tunnel face is excavated to YK5 + 680. When the elastic modulus of the surrounding rock changes from 2 to 1 GPa, the maximum settlement of the roof increases by 0.23 cm, from 0.6 to 0.83 cm. When the elastic modulus of the surrounding rock changed to 0.2 GPa, tunnel deformation increased by 0.67 cm, from 0.83 to 1.5 cm. Then, the tunnel deformation is slowly deformed under the excavation disturbance of the tunnel face. The tunnel deformation is more sensitive to different surrounding rock properties than different over-excavation heights.
While the distance between the left and right tunnel faces are taken as control variables, the tunnel deformation is shown in Figures 9E,F. In the three cases, the maximum deformation of the tunnel crown is 3.0 cm, the convergence is 1. cm, the minimum is 1.48 cm, and the convergence is 0.6 cm. When the distance between the left and right tunnel faces changes from 4-8 m, the tunnel deformation decreases, the crown settlement decreases from 3.0 to 2.5 cm, and the maximum convergence decreases from 1.2 to 1.1 cm. However, when the tunnel face distance changes from 8 to 16 m, the tunnel deformation amplitude is more significant than 4-8 m. The crown settlement decreases from 2.5 to 1.4 cm, and the maximum convergence decreases to 0.75 cm. Compared with the deformation under different surrounding rock strength, the difference of tunnel deformation is not apparent, but it is higher than that of tunnels with different over-excavation heights.
Analyzing all the figures, the distance between the left and right tunnel faces, surrounding rock strength, and overexcavation height all affect the crown settlement and horizontal convergence. The image results show that the variation range of tunnel deformation is small when the overexcavation height is changed. When the surrounding rock properties and the distance between the left and right tunnel faces are altered, the deformation amplitude of the tunnel rises. Nonetheless, the factor which has the most significant impact requires careful identification. Statistics of the results under different conditions are shown in Table 4, using grey relation analysis and entropy weight method to judge.
According to the equations of the chapters earlier, the relation degree between the reference sequence and the comparison sequence is gained. The relation degree value is sorted after calculating which factors are the main influencing factors and relative secondary factors of collapse risk. If r 1 <r 2 the comparison sequence X 1 has a more significant impact on the reference sequence. The final relation degrees of the three factors on crown settlement and tunnel convergence are shown in Tables 5.
The entropy method is used to calculate the weight of crown settlement and horizontal convergence. According to the calculation results, the information entropy of crown settlement is 0.88 and the horizontal convergence is 0.906. Hence, the proportion of crown settlement is 0.56 and the proportion of horizontal convergence is 0.44。According to equation 12, the grey correlation coefficients of the aforementioned factors are combined to obtain the final grey correlation coefficient. All of these calculation results are shown in Figure 10. The relation degree corrected using the entropy method is more objective and accurate. It can be observed that tunnel deformation is influenced by three factors, and their combined influence on tunnel deformation is greater than 0.6. Over-excavation height is the most important element impacting tunnel deformation, with a relation degree of 0.858. This suggests that tunnel collapse is mostly caused by an over-excavation incident. The tunnel's deformation amplitude under various over-excavation heights is, likewise, the greatest. The surrounding rock is the second, with a relation degree of 0.721. The tunnel deformation is also influenced by the distance between the left and right tunnel faces, and the relation degree is 0.7. The other two elements, as can be seen, are also involved in the collapse disaster.

CONCLUSION
A severe collapse occurred at Huangjiakuang tunnels in Weihai City, Shandong Province, which led to a 27.5 m 2 subsidence area on the ground, resulting in delays and economic losses. In this study, the main influencing factors to collapse were studied by combining the advanced geological prediction of working face, excavation characteristics, and construction disturbance. The entire process of collapse to the surface is inversed using numerical simulation software. Based on the data obtained by simulation, the three factors are further analyzed using grey relation analysis and entropy weight method. Some conclusions and numerical insights are summarized as follows: 1. Combined with the advanced geological prediction, geological survey data, and construction records, the whole process of tunnel collapse is inversed by numerical simulation. Before the collapse, the stress in the collapse area is increased to 1.1-1.2 times the initial value. When the collapse happened, the maximum crown settlement in the collapse area is 11 cm, and the surrounding rock bears a tensile stress. Finally, the soil is damaged by tension and evolves into collapse. As a result, it is determined that the collapse is caused by a combination of excavation disturbance and an over-excavation accident in a weak surrounding rock condition. 2. The control variables of surrounding rock property, overexcavation height, and the distance between the left and right tunnel faces are studied, and the corresponding tunnel deformation charts are obtained. It is demonstrated that the three factors quantitatively affect the tunnel deformation. The three factors are quantitatively analyzed using grey relation analysis and entropy weight method and their grey correlation degree exceeds 0.6. The relation degree of over-excavation height to deformation is 0.858. The relation degree of the surrounding rock is 0.721, and the left and right tunnel distance is 0.700. Consequently, the main reason for judging the collapse accident is the tunnel over-excavation.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary Materials; further inquiries can be directed to the corresponding author.