Nonlinear Response of Tunnel Portal under Earthquake Waves with Different Vibration Directions

Tunnel portal sections often suffer serious damage in strong earthquake events. Earthquake waves may propagate in different directions, producing various dynamic responses in the tunnel portal. Based on the Galongla tunnel, which is located in a seismic region of China, three-dimensional seismic analysis is conducted to investigate the dynamic response of a tunnel portal subjected to earthquake waves with different vibration directions. In order to simulate the mechanic behavior of slope rock effectively, an elastoplastic damage model is adopted and applied to ABAQUS software by a self-compiled usermaterial (UMAT) subroutine. Moreover, the seismic wave input method for tunnel portal is established to realize the seismic input under vertically incident earthquake waves with different vibration directions, e.g., S waves with a vibration direction perpendicular or parallel to the tunnel axis and Pwaves with a vibration direction perpendicular to the tunnel axis. The numerical results indicate that the seismic response and damagemechanisms of the tunnel portal section are related to the vibration direction of the earthquake waves. For vertically incident S waves running perpendicular to the tunnel axis, the hoop tensile strain at the spandrel and arch foot and the hoop shear strain at the vault and arch bottom are the main contributors to the plastic damage of the tunnel. The strain is initially concentrated around the tunnel foot and spandrel, before shifting to the tunnel vault and bottom farther away from the tunnel entrance. For vertically incident S waves running parallel to the tunnel axis, very large hoop shear strain and plastic damage appear at the tunnel haunches. This strain first increases and then decreases with distance from the tunnel entrance. For vertically incident P waves running perpendicular to the tunnel axis, the maximum damage factor of the slope rock and the maximum plastic strain of the tunnel are significantly lower than for S waves. Moreover, with increasing distance from the tunnel entrance, the plastic damage to the tunnel lining rapidly decreases.


Introduction
The thousands of mountain tunnels that have recently been built in China play an important role in local development. Such tunnels may pass through highly active seismic areas, especially in Southwest China. Several strong earthquakes in recent decades, some mountain tunnel portals suffered particularly serious damage to some mountain tunnel portals, which is just less critical than that of the fault [1][2][3][4]. Slope instability, rock collapse and lining cracks are triggered by strong earthquakes, introducing numerous challenges to the safe operation of these tunnels. Thus, the tunnel portal is one of the weaker components of the tunnel structure in the event of an earthquake [5,6].
There has been considerable research on the seismic damage suffered by tunnel portals, mainly through model tests, numerical simulations, and in-situ investigations. Sun et al. [7] performed model tests to investigate the seismic response of a tunnel portal, and analyzed the interface interaction characteristics of the surrounding rock and the tunnel lining. To study the failure mode of a tunnel portal under strong earthquakes, Zhao et al. [8] adopted a rigid box model performed shaking table tests of the tunnel portal section. Xu et al. [9] conducted further shaking table tests to study the seismic response of the tunnel portal under different seismic scenarios. Compared with model tests, numerical simulations require less time and labor resources, and are capable of simulating a larger range of environmental and boundary conditions. Wu et al. [10] studied the damage mechanisms of a tunnel portal numerically, while Liu et al. [11] conducted a numerical study of the seismic response of the tunnel portal at the Dianzhong water diversion project. Importantly, Wang et al. [12] proved that the hydrodynamic pressure has a significant adverse effect on the seismic capacity of the tunnel portal. Other excellent research has been reported by Geni [13] and Yu et al. [14], Yu et al. [15] and so on. These preliminary works make us have more knowledge of the seismic performance and failure mechanism of the tunnel portals.
Because of the smaller overburden depth, the tunnel portal has a seismic response that is strongly affected by the slope around the structure. The slope rock is often severely weathered, reducing its ability to restrain the portal structure compared with other tunnel parts that are buried deeper underground [12]. Therefore, a precise and appropriate constitutive model is essential in simulating the failure strength, plastic deformation, damage softening, and other seismic mechanical behavior of slope rock. However, previous works have only considered heavily simplified constitutive models, such as linear elastic models, plastic models based on the Mohr-Coulomb criterion, and so on. Moreover, earthquake waves propagate through the ground in different vibration directions. Although many researchers have studied the seismic response of deep-buried tunnels under earthquake waves propagating in different vibration direction [16][17][18][19], the seismic response and damage characteristics of tunnel portals under such conditions have not been fully considered.
In this study, an elastoplastic damage model is integrated into the ABAQUS software through a self-compiled user material (UMAT) subroutine, allowing the mechanical behavior of the slope rock to be simulated. First, the accuracy of the model is verified. Based on the Galongla tunnel project, a dynamic time-history model is then constructed to simulate the seismic response. A viscous-spring boundary is used to absorb the wave radiation. In addition, a wave input method for the slope rock is established, and the applicability and accuracy of this method are verified. Vertical incident seismic waves in different vibration directions are simulated by the proposed method, and the response of the portal structure of Galongla tunnel is numerically studied. Finally, the seismic response characteristics and damage mechanism of the tunnel portal are discussed.

Galongla Tunnel Project
Galongla tunnel is a critical part of the Zhamo highway in Tibet, the only highway leading to Motuo County in China. the local geography is shown in Fig. 1a. The geological condition of the tunnel is very complex. As shown in the geological profile in Fig. 1b, the rock mass along the tunnel can be divided into Grades II, III, and IV according to the Chinese rock classification code. The tunnel portal is loose Holocene deposits consisting of welded tuff, gravel, silt, and clay, as shown in Fig. 1c. The cross-section of the Galongla tunnel is horseshoe shaped, as shown in Fig. 2. The tunnel structure consists of rock bolts, a primary lining, a waterproof layer, and a secondary lining, and was excavated and constructed by the New Austrian tunneling method.  Galongla tunnel is located in the Himalayan fault zone, which is subjected to frequent geological activity, i.e., an earthquake-prone area. According to the seismic data issued by the China Earthquake Administration, the basic intensity degree is VIII. Therefore, the tunnel, especially the tunnel portal, faces a serious earthquake risk. To investigate the seismic damage mechanics of the tunnel portal part, the results of dynamic numerical analysis are described in the following sections.

Elastoplastic Damage Model for Rocks and Realization with the ABAQUS Software
Based on the nonlinear unified strength criterion, Du et al. [20] established a triaxial elastoplastic damage model for rock mass. The plastic part of the model employs the nonlinear unified strength criterion as the yield function in the effective stress space. The hardening law takes the equivalent plastic shear strain as the hardening parameter, and considers the effect of mean stress on the hardening rate. For the damage part, the evolution law is assumed to have an exponential form so as to consider the degradation and softening behavior of the rock mass. In this study, the constitutive model is implemented in the ABAQUS software through a self-compiled UMAT subroutine to simulate the mechanical behavior of the rock mass [21].

Description of the Elastoplastic Damage Model
The following nonlinear unified strength criterion developed by Lu et al. [22] is adopted as the yield function of the plastic part: where κ is the hardening variable;σ i and M β are the characteristic stress and the failure strength ratio, respectively, which can be expressed as where the reference stress p r takes a value of 1 MPa to give a dimensionless stress; β denotes the parameter for converting the normal stressσ i to the characteristic stressσ i , which is constant for a certain material; ϕ c denotes the internal friction angle.
To describe the rock dilatancy effectively, a non-associated flow rule is adopted. The expression for the plastic potential function is similar to the yield function, with the parameter M β replaced by M g , which is a function of the dilatancy angle ψ: The hardening variable κ is a function of the equivalent plastic shear strain γ p , and has the form where a o is the initial hardening rate; σ c and σ o are the uniaxial compression strength and the three-dimensional tensile strength, respectively; and A γ is the plastic shear strain corresponding to the hydrostatic pressure p, which has a value of σ c /3. The equivalent plastic shear strain γ p can be calculated as where e p ij is the plastic deviatoric strain; ε p ij and ε p v denote the plastic strain tensor and the plastic volume strain, respectively; δ ij is Kronecker's delta.
The damage factor d is introduced to characterize the micro-defects in the surrounding rock material and the development of its internal damage state. There are two main forms of the damage factor, i.e., scalar damage variables and tensor damage variables. The damage variable in scalar form describes the isotropic damage of a material, and is mainly used for the distribution of spherical defects. It is considered that the damage degree of a material is the same in any direction. However, the direction of crack propagation for rock materials often has a certain directionality under the action of external load, and the influence of these cracks on the crosssectional area will be different in each direction. Given the anisotropy of rock damage, many researchers use tensor damage variables to characterize the directionality of crack distribution and propagation. In establishing a rock damage evolution equation, many scholars have established models based on different damage internal variables [23][24][25][26][27][28][29][30][31][32][33]. Considering that the rock damage is softened by the volume expansion caused by the development of internal micro-cracks, this paper uses the volume strain as the internal variable of the evolution of damage. The specified damage variable is positive in compression and negative in tension. In the damage part of the model, the damage variable d is defined as where b 1 is a parameter controlling the slope of the softening curve; d is the equivalent strain, which can be written as where d pv is the volumetric plastic strain rate; χ 1 controls the rate of increase of d ; b 2 is a constant that affects the value of χ 1 .

Relationship between the Stress and Strain
For the elastoplastic damage model, the increment in the Cauchy stress σ ij can be expressed as the effective stressσ ij : The relationship between increments in the effective stressσ ij and the stain ε ij is where ε e ij and ε p ij denote the elastic strain and the plastic strain, respectively; D e ijkl is the elastic stiffness matrix, which has the form The plastic strain increment dε p ij can be obtained by the following flow rule: Substituting Eq. (11) into Eq. (5) yields The consistency condition then gives Substituting Eq. (12) into Eq. (13), we obtain where Combining Eqs. (9), (11), and (14), we find that Combining Eqs. (9), (11), and (16) gives, after some rearrangement, the stress-strain increment relationship: The expressions for ∂f ∂γ p , ∂f ∂σ ij , and ∂g ∂σ ij in Eqs. (15) and (17) are given in the Appendix.

Implementation in ABAQUS Software and Verification
Based on the derived relationship between the stress and strain, the elastoplastic damage model is integrated into the ABAQUS software through a self-developed UMAT subroutine. To verify the accuracy of the implemented model, a three-dimensional (3-D) solid unit cube model measuring 1 m × 1 m × 1 m is built and discretized into a single grid using reduction integral elements (C3D8R). The numerical results are compared with the results of conventional triaxial experiments by Kawamoto et al. [34] and Lin et al. [35]. The constitutive model parameters of the triaxial experiments performed by Kawamoto et al. [34] and Lin et al. [35] are listed in Table 1.
Lin et al. [ [34] and Lin et al. [35]. From this figure, it is clear that the constitutive model accurately simulates the whole process of rock deformation and damage under a complex stress state, especially the strain softening and the changing characteristics of strength with respect to the confining pressure.

Geometry of the Computational Model
To investigate the seismic response of the Galongla tunnel portal, dynamic history simulations were performed using the ABAQUS software. The 3-D finite element model and its dimensions are shown in Fig. 4, in which the slope inclination is set to be 4:5 in height and length along the tunnel axial direction. The Mohr-Coulomb constitutive model was adopted for the nonlinear behavior of the tunnel lining, with material parameters of density ρ = 2500 kg/m 3 , Young's modulus E = 30 GPa, Poisson's ratio υ = 0.22, friction angle ϕ = 58.7 • , and cohesive stress c = 2.3 MPa. The rock medium was simulated by the elastoplastic damage model, implemented in the ABAQUS software as described in Section 3, with material parameters of ρ = 2300 kg/m 3 , In this study, the tunnel portal responses under three different earthquake incident waves are investigated. The first case considers the incident waves to be plane P waves with a vibration direction perpendicular to the tunnel axis. The second and third cases consider the earthquake waves to be vertically incident S waves with vibration directions perpendicular and parallel to the tunnel axis, respectively. The Kobe University earthquake records for rock strata were obtained in similar geological conditions to those of the Galongla tunnel project site. The Kobe waves are real seismic records with rich spectra, as shown in Fig. 5. They have been widely used to study the seismic response of mountain tunnels through numerical simulations and experimental tests. Therefore, the Kobe waves are selected as the input waves in this study. During the numerical analysis, a viscous-spring artificial boundary is employed to simulate the radiation damping effect at the truncated boundary, as introduced in Section 4.2. The seismic input method for the slope site considering different vibration directions is then described in Section 4.3.

Viscous-Spring Artificial Boundary
In recent decades, many artificial boundaries have been developed which can be found in the works of Lysmer et al. [36], Deeks et al. [37], Du et al. [38], and so on. The viscous-spring boundary developed by Du et al. [38] is adopted in this study due to its simple form and high accuracy. The boundary is realized by establishing springs and dampers in three directions, as shown in Fig. 6a. At boundary node l of the 3-D finite element model, the parameters of the springs and dashpots are where the subscripts n and s denote the normal and tangential direction of the boundary plane, respectively. The distance between the wave source and the artificial boundary is denoted by R; λ, G, and ρ denote the Lamé constant, shear modulus, and mass density, respectively; c s is the shear wave velocity in the medium, and c p is the compression wave velocity; a r and b r denote the modified coefficients, which have values of 0.8 and 1.1, respectively. Fig. 6a describes the incident conditions for the truncated region of the slope site. The earthquake wave motions can be converted into equivalent nodal forces, which are applied at the artificial boundary [39][40][41]. Together with the viscous-spring boundary developed by Du et al. [38], the equivalent nodal forces at a given boundary node l can be obtained from Eqs. (19)- (21). The specific forms are as follows.

Seismic Input Method
For node l on the bottom boundary (y = 0), the equivalent nodal forces are where F x , F y , and F z are the equivalent nodal forces in the x, y, and z directions, respectively; u x , u y , and u z are the free field displacements of boundary node l in the x, y, and z directions, respectively;u x ,u y , andu z are the free field velocities; σ yy , σ yx , and σ yz are the free field stresses; A 3l is a quarter of the total area of all boundary elements containing boundary node l, as shown in Fig. 6a For node l on the left boundary (x = 0) and the right boundary (x = L x ), the equivalent nodal forces are where δ is a boundary-dependent parameter, i.e., δ = −1 for the left boundary and δ = 1 for the right boundary; σ xx , σ xy , and σ xz are the free field stresses on the left and right boundaries, respectively.
For node l on the back boundary (z = 0) and the front boundary (z = L z ), the equivalent nodal forces are where δ takes a value of −1 on the back boundary and a value of 1 on the front boundary; σ zz , σ zx , and σ zy are the free field stresses on the front and back boundaries, respectively.
As shown in Eqs. (19)-(21), the free field displacements u x , u y , and u z , velocities ...,u y , anḋ u z , and stresses σ xx , σ yy , σ zz , σ xy , σ zy , and σ zx need to be determined to calculate the boundary equivalent nodal forces. The free field displacements, velocities, and stresses of each boundary node can be determined through free field site analysis. In consideration of the homogeneity and isotropy of the site and the symmetry of the loading conditions, the free field response of the slope site can be analyzed using a two-dimensional (2-D) site model, as shown in Fig. 6b. The displacement time histories, velocity time histories and stress time histories of each boundary node are provided. In the 2-D model for slope seismic analysis, the viscous-spring boundary developed by Du et al. [38] is applied at the bottom and lateral boundaries. At boundary node l of the 2-D model, the spring and dashpot parameters can be written as , C 2ls =b r ρc s In addition, the incident earthquake waves are converted into equivalent nodal forces that apply at each boundary node. For node l on the bottom boundary (y = 0) of the 2-D model, the equivalent nodal forces of incident SV, SH, and P waves are given by where u 0 (t) andu 0 (t) denote the displacement and velocity of the incident waves. A 2l is half of the total length of all boundary elements containing boundary node l in the 2-D finite element model, as shown in Fig. 6b, where A 2l = (L 1 + L 2 )/2.
For node l on the left boundary (x = 0) and the right boundary (x = L x ), the equivalent nodal forces of the incident waves are given by where the boundary-dependent parameter δ takes a value of −1 on the left boundary and a value of 1 on the right boundary; h is the height of the lateral boundary.
In general, the seismic input method for the 3-D model is as illustrated in Fig. 6. First, the free field response of the slope site is numerically simulated by a 2-D site model. Using the displacement, velocity, and stress time histories from the 2-D site model, the equivalent nodal forces of the 3-D model are calculated with Eqs. (19)-(21).

Validation of Seismic Input Method
The proposed seismic input method is now verified. The 3-D slope site model illustrated in Fig. 7 was constructed with material parameters of ρ = 2300 kg/m 3 , E = 0.87 GPa, and υ = 0.42. The input seismic waves take the form of an impulse, with the displacement time history shown in Fig. 8. With the aid of the present seismic input method, the impulse is transmitted to the 3-D slope site model as vertically incident P waves and vertically incident S waves with vibration directions perpendicular and parallel to the x-direction, respectively. In addition, a 2-D slope site model was constructed to provide reference solutions.  Fig. 9. For reasons of brevity, only the displacement contours under vertically incident P waves are shown in this figure. It can be observed from Fig. 9 that the displacement contours of the 3-D model are very similar to those of the 2-D model. Moreover, the displacement in the y-direction at monitoring points A and B is plotted in Fig. 10. From Fig. 10, it is apparent that the 3-D model provides good agreement with the 2-D model. Thus, the incident method of the 3-D slope site proposed in this paper has a high degree of accuracy, and can thus be employed to investigate the seismic response of the tunnel portal.  Fig. 11 shows the damage factor contours of the rock ground near the tunnel portal under S waves with a vibration direction perpendicular to the tunnel axis. The surrounding rock suffers heavy damage, especially near the tunnel-ground interface. The damage to the surrounding rock near the tunnel foot and spandrel is significant near the tunnel entrance. With increasing distance from the entrance, the heavy damage gradually shifts from the tunnel foot and spandrel to the tunnel vault and bottom. The maximum principal plastic strain contours of the tunnel lining under S waves with a vibration direction perpendicular to the tunnel axis is presented in Fig. 12. The dynamic response of the tunnel lining is mainly affected by the tunnel-foundation interactions, the mechanics of which are significantly different from those of upper structures built on open land. Thus, the plastic deformation of the tunnel matches the progress of damage in the surrounding rock (see Fig. 11). In detail, at the tunnel entrance, the plastic strain is concentrated around the tunnel foot and spandrel. With increasing distance from the entrance, the plastic strain shifts from the tunnel foot and spandrel to the tunnel vault and bottom. The strain can be decomposed into the hoop axial strain h , hoop shear strain γ h , and axial tensile strain a , as shown in Fig. 13. Fig. 14 shows the distributions of these strain components at tunnel cross-sections of D = 7.5 m and 61.5 m, respectively, using L-foot to represent the lower left arch foot position of the lining, using L-haunch to represent the left arch waist position of the lining, using L-spandrel to represent the upper left arch shoulder position of the lining, using R-foot to represent the lower right arch foot position of the lining, using R-haunch to represent the right arch waist position of the lining, using R-spandrel to represent the upper right arch shoulder position of the lining, using Bottom to represent the bottom center position of the lining, and using Vault to represent the top center position of the lining. The same goes in the following. Significant hoop axial strain h appears at the tunnel foot and spandrel, and the hoop shear strain γ h around the tunnel bottom and vault is considerable. Though the rock-tunnel model is bilaterally symmetric, the shaking produced by an earthquake is not symmetric, and neither are the vertically incident S waves with a vibration direction perpendicular to the tunnel axis. Thus, h and γ h are not symmetric at the cross-sections. Compared with h and h , the axial tensile strain a is very small, which indicates that the tunnel suffers only slight tensile deformation in the tunnel axial direction. The distributions of hoop axial strain h at the left spandrel and right arch foot along the tunnel axis are given in Fig. 15a, and the hoop shear strain γ h at the tunnel vault and bottom are shown in Fig. 15b. As the axial tensile strain a is very small, the distribution of this component along the tunnel axis is not presented. From Fig. 15a, h at the spandrel and arch foot are dominant at the tunnel entrance, and then quickly decrease to a stable value with distance from the tunnel entrance. Therefore, the plastic strain at the spandrel and arch foot are higher around the tunnel entrance than at other positions, as shown in Fig. 12. Fig. 15b shows that the hoop shear strain first increases and then decreases with increasing distance from the tunnel entrance. Thus, the plastic damage changes to become concentrated around the tunnel vault and bottom, as shown in Fig. 12.  Fig. 16 shows the damage factor contours of the rock ground near the tunnel portal under S waves with a vibration direction parallel to the tunnel axis. The surrounding rock suffers severe damage near the tunnel entrance. Moreover, the damage contours of the rock are symmetrical. Fig. 17 shows the maximum principal plastic strain contour at the tunnel lining under S waves with a vibration direction parallel to the tunnel axis. The plastic strain of the tunnel lining in Fig. 17 is concentrated around the tunnel haunches, and the value increases and then decreases with distance from the tunnel entrance. Moreover, the rock damage and the plastic strain on the tunnel under S waves with a vibration direction parallel to the tunnel axis are different from those in the case of S waves vibrating perpendicular to the tunnel axis. The distributions of the strain components at two tunnel cross-sections for S waves with a vibration direction parallel to the tunnel axis are shown in Fig. 18. This figure shows that the tunnel suffers significant hoop shear strain, but only slight hoop axial strain and axial tensile strain. Moreover, the hoop shear strain is much higher around the tunnel haunches than at other positions. This indicates that the tunnel haunches are vulnerable to S waves with a vibration direction parallel to the tunnel axis. The distribution of the hoop shear strain changes at the same pace as the plastic damage at the haunches (see Fig. 17). Fig. 19 plots the distribution of the hoop shear strain at the haunches along the tunnel axis. From Fig. 19, we see that the hoop shear strain increases at first and then decreases with increasing distance from the tunnel entrance, and changes at the same pace as the plastic strain of the tunnel (see Fig. 17).

Vertically Incident P Waves with Vibration Direction Perpendicular to Tunnel Axis
The seismic response of the tunnel portal under vertically incident P waves with a vibration direction perpendicular to the tunnel axis is now analyzed. The damage factor contours of the rock near the tunnel portal are shown in Fig. 20, and the maximum principal plastic strain contour of the tunnel lining is shown in Fig. 21. The damage contours of the rock, as well as the plastic strain on the tunnel lining, are symmetric and significantly concentrated around the tunnel entrance. With increasing distance from the tunnel entrance, the plastic damage to the tunnel lining quickly decreases. The maximum damage factor of the rock and the maximum plastic strain of the tunnel are significantly lower than in the case of vertically incident S waves.  22 shows the distributions of strain components along two tunnel cross-sections for vertically incident P waves. The hoop axial strain h and hoop shear strain γ h are larger than the axial tensile strain a . Moreover, Fig. 22 indicates that h and γ h are dominant around the tunnel haunches, and change at the same pace as the plastic strain concentrated around the tunnel haunches. Fig. 23 shows the distributions of h and γ h at the tunnel haunches along the tunnel axis. With increasing distance from the tunnel entrance, the hoop axial strain gradually decreases, whereas the hoop shear strain increases at first and then decreases.

Summary and Conclusions
During past earthquake events, tunnel portals have suffered more serious damage than other tunnel sections. Earthquake waves may propagate in different vibration directions, producing a complex dynamic response in the tunnel portal. Based on the Galongla tunnel project, which is located in a seismic region of China, a 3-D numerical model of the tunnel portal was established to study such effects. The seismic mechanical behavior of the slope rock was simulated using an elastoplastic damage model implemented in ABAQUS through a self-compiled UMAT subroutine. The viscous-spring artificial boundary was used to simulate the radiation damping effect at the truncated boundary. Moreover, a seismic input method for the slope site was established for vertically incident earthquake waves with different vibration directions. The damage characteristics of the Galongla tunnel portal under different vibration directions of earthquake waves were then numerically studied. The conclusions are as follows: (1) For vertically incident S waves with a vibration direction perpendicular to the tunnel axis, the hoop tensile strain at the spandrel and arch foot and the hoop shear strain at the vault and arch bottom are the main contributors to the plastic damage of the tunnel. For vertically incident S waves with a vibration direction parallel to the tunnel axis, a very large hoop shear strain occurs at the tunnel haunches. For vertically incident P waves with a vibration direction perpendicular to the tunnel axis, the tunnel is mainly subjected to hoop tensile strain and hoop shear strain. (2) The vibration direction of the earthquake waves significantly influences the dynamic response of the tunnel portal. For vertically incident S waves with a vibration direction perpendicular to the tunnel axis, the hoop axial strain at the spandrel and arch foot are dominant at the tunnel entrance, and then rapidly decrease to a stable value away from the tunnel entrance. The hoop shear strain increases at first and then decreases with increasing distance from the tunnel entrance. For vertically incident S waves with a vibration direction parallel to the tunnel axis, the hoop shear strain increases and then decreases with distance from the tunnel entrance. For vertically incident P waves with a vibration direction perpendicular to the tunnel axis, the hoop axial strain decreases gradually, and the hoop shear strain first increases and then decreases with increasing distance from the tunnel entrance.