Simulation of Surface Settlement Induced by Parallel Mechanised Tunnelling

Abstract

Mechanised tunnelling is extensively utilised for twin tunnel construction, particularly in urban areas. A common challenge encountered during this construction method is the occurrence of surface settlement (SS) induced by tunnelling activities. The integrity of nearby structures can be compromised by SS, making it imperative to accurately quantify and mitigate this phenomenon. Several methods for determining SS exist, including empirical formulas and laboratory studies. However, these methods are often constrained by specific soil types and are time-consuming. Moreover, crucial parameters such as tunnel operational factors and construction stages are often omitted from empirical formulas. Given these limitations, this paper aims to address these challenges by employing 3D numerical analysis to simulate tunnelling-induced SS in twin tunnels. This approach takes into account tunnel geometry, construction sequencing, soil properties, and tunnelling operational factors. By incorporating data from in-situ and laboratory tests conducted on the ground, engineering soil parameters are established as inputs for the numerical analysis. The simulated SS results obtained from the 3D numerical analysis are compared with field measurements of SS taken from available ground surface settlement markers. The transverse SS pattern derived from the numerical analysis closely mirrors the field measurements. Additionally, SS values above the first and second tunnels are compared with field measurements, resulting in coefficient of determination (R²) values of 0.94 and 0.96, respectively. The utilisation of the 3D numerical modelling approach enables the customizable mitigation strategies for managing the SS with project-specific parameters such as tunnel geometry, geotechnical engineering factors, and tunnelling operational variables. This will help plan and construct more sustainable tunnels with minimal effects on the ground and residential areas.

1. Introduction

Tunnelling is one of the most significant transport solutions where the overlying population does not need to be displaced. With the advancement of current technologies, mechanised excavations to construct underground spaces such as tunnels have become popular. Although mechanised tunnel excavations have been widely used all around the world, considerable safety considerations during construction need to be ensured. Several problems can be encountered during tunnel excavations, such as tunnel face instability, excessive wear of the cutter head, and excessive surface settlement (SS) [1,2]. Among them, SS induced by mechanised tunnelling is still one of the common issues in tunnel construction [3,4]. The excessive SS during and after tunnelling projects has an adverse impact on the existing structures. Therefore, it is important to estimate the settlements caused by the tunnelling to minimise the effects on the existing structures. Several methods can be used to estimate the SS due to tunnelling, which is empirical or semi-empirical, laboratory-based, and numerical analysis. One of the pioneer empirical formulas to determine the SS due to tunnelling was proposed by Peck [5], which shows that the volume loss is from the radial deformations along the tunnel perimeter directly linked to the SS. Based on the field observations on the site and the simplification of the equation proposed by Litwiniszyn [6], Peck [5] suggested Equation (1) to predict SS induced by tunnelling for soft clay. The settlement in this equation is calculated based on the pattern of SS caused by loss of ground and approximated by a Gaussian probability curve.

$$S=S_{max}{e}^{\frac{-x^2}{2i^2}}$$

(1)

where $S$ is the surface settlement in the transverse section at a specific distance, $x$ is the distance from the centreline of the tunnel, and $i$ is the point of inflection (settlement through). In this way, the maximum SS can be defined using Equation (2):

$$S_{max}=\frac{V_s}{ \sqrt{2 x \pi } x i}$$

(2)

where V_s is the volume loss of the soil (m³/m), and can be expressed by Equation (3):

$$V_s=\frac{Volume Loss \left(VL\right)(\%)}{100}\left(\frac{\pi D^2}{4}\right)$$

(3)

The main inputs for the calculation of the maximum SS are affected by volume loss of the soil (VL) and settlement trough. These volume losses are affected by the type of ground, tunnel geometry, and ground condition [7,8]. Hence, many researchers have carried out their investigations into different types of ground conditions with various ranges of VL, as summarised in

Table 1. Tunnelling methods with various ground conditions of VL.

Author (s)	Ground Condition	VL (%)	Method of Tunnelling
Attewell and Farmer [9]	London clay	1.44	Hand excavation shield tunnelling
O’Reilly and New [10]	London clay	1.0–1.4	Open face shield-driven tunnels
Mair and Taylor [11]	Stiff clay	1.0–2.0	Open face method
	Stiff clay	0.5–1.5	NATM
	Sand	0.5	Closed face Tunnelling Boring Machine
	Soft clay	1.0–2.0	Closed face Tunnelling Boring Machine
Wan et al. [12]	London clay	0.8	Earth Pressure Balance Machine (EPBM)
Amir and Mohammad [13]	Graded gravel to silt/clays	0.2–0.7	EPBM
Le et al. [14]	Sandy	<0.2 to 2.4	EPBM

. In addition, several researchers have also proposed various i equations for different ground conditions, as presented in

Table 2. Summary of empirical formulas for estimation of settlement trough width.

Authors	Empirical Formula	Variable Definition	Ground Condition	Tunnelling Excavation Method
Peck [5]	$i=R\prime {(\frac{Z}{2R\prime })}^n$	R’ is the radius of the tunnel. Z is the tunnel depth below ground level. n is a constant parameter dependent on soil type (0.8–1).	Various types of soils	Open cutting excavation
O’Reilly and New [10]	$i=kZ$	k is the constant parameter dependent on the soil type. Z is the tunnel depth below ground level.	Various types of soils	Shield tunnelling
Mair et al. [15]	$i=Z[0.175+0.325(1-\frac{Z\prime }{Z}$)	Z’ is the depth of the calculated settlement trough from the surface settlement. Z is the tunnel at a depth below ground level.	Clay	Centrifuge model test
Loganathan and Poulos [16]	$i=R\prime \left(1.15\right){(\frac{Z}{2R\prime })}^{0.9}$	R’ is the radius of the tunnel. Z is the tunnel at a depth below ground level.	Clay	Tunnelling machine
Wang et al. [17]	$i=m (R+Z\tan \left(45-\frac{\Phi }{2}\right))$	m is the influence coefficient of the width. R’ is the radius of the tunnel. Z is the tunnel depth below ground level. Φ is the internal friction angle of the soil.	Sand	Laboratory
Zhu et al. [18]	$i=0.51\mathrm{Z}+0.48$	Z is the tunnel depth below ground level.	Sand and Clay	Shield machine

.

According to Table 2, VL and i have a wide range of values, which can be due to different methods of tunnel construction and ground conditions. It can also be seen that most studies are carried out for clay material. In addition, no tunnelling operational parameter is included in the previous empirical formulas. Loganathan and Poulos [16] proposed a semi-empirical equation (i.e., Equation (4)), which is based on the gap parameter proposed by Lee et al. [19] and the exponential function to model the nonlinear ground movement. This formula has been checked with the five case studies. From the findings, four tunnels that were constructed in stiff clay showed good agreement; however, soft clay depicted a much lower estimation of the ground loss. Nevertheless, the method proposed by Loganathan and Poulos [16] is limited to clay and based on an elastic solution.

$${\epsilon }_{xy}=\frac{4gR+g^2}{4R^{\prime 2}} e^{\frac{1.38x^2}{{\left(Z+R^{\prime }\right)}^2}+\frac{0.69z^2}{Z^2}}$$

(4)

where R′ is the radius of tunnel, Z is the tunnel depth below ground level, z is the depth of the point, g is the gap parameter, and ε_xy is the equivalent ground loss component at the tunnel soil interface due to the ground movements at point x, z.

Aside from the empirical and semi-empirical equations, 2D finite element analysis has been widely used to determine SS induced by tunnelling. Some of these 2D simplified techniques are the contraction method by Vermeer and Brinkgreve [20], the stress reduction method by Addenbrooke et al. [21], and the modified grout pressure method by Surarak [22]. The contraction method of the analysis works with the introduction of the correction factor (which represents the value of reduction and excavated area) as a predefined uniform radial inward strain. However, this assumption may not reflect the actual field displacement because, based on the centrifuge modelling study carried out by Mair [23], the finding shows that little ground displacement happens at the tunnel invert. In addition, the stress reduction method uses an “unloading factor” to stimulate the 3D tunnelling effects in the 2D numerical model, but this method requires trial-and-error on the unloading factor to calculate/determine the SS induced by tunnelling. Additionally, the modified grout pressure method [22] is the method that considers the tunnel operation parameter for the analysis. Although this method can include the operational parameters (face pressure and grout pressure), it does not consider the 3D mechanised tunnelling sequence in the analysis. Every stage of the mechanised tunnelling construction influences the ground SS. The stages of the tunnelling such as excavation using face pressure, advancement of the tunnelling machine, installation of the lining, injection of the tail void grout, and hardening of the grout, affect the surrounding soil area [24,25]. Apart from that, the type of soil model used in the numerical analysis is important; therefore, proper in-situ and laboratory data interpretation is required to stimulate the closest possible ground conditions at the site. In 3D numerical analysis, the tunnelling construction can cause SS at different stages of the construction. Mathew and Lehane [26] deduced the inconsistency of the calculations using volumetric contraction in comparison with field measured settlements due to the complexities of the actual tunnelling procedures. In comparison with a single tunnel, twin tunnels are more complex and more difficult to analyse as more parameters e.g., the spacing between two tunnel and different tunnelling operational parameters are required to consider in the analyses.

Most studies [3,27] investigated these issues for single tunnels; however, limited studies are available for twin tunnels using earth pressure balance (EPB) shield machines. In the recent numerical tunnelling analysis carried out by Islam and Iskander [28], the effect of geometric parameters and construction sequence on the SS was investigated using MIDAS GTS NX software 2021 (v1.1). However, this case study is not for parallel twin tunnelling and the outputs were not compared with the field measurements. The 3D numerical model is not limited to certain types of ground in comparison with the empirical and semi-empirical methods. In addition, the 3D numerical model requires a shorter period of time to determine the SS than the laboratory works. The 2D numerical methods (contraction and stress reduction methods) do not consider operational parameters in the analysis, whereas modified grout pressure does not consider the construction stages of the tunnelling works such as the advancement of the tunnelling machine. The consideration of important factors affecting the SS due to tunnelling in the analysis is crucial to providing the closest possible results with the actual measurements.

In light of the above discussion, this study considers the holistic approach of 3D modelling shield tunnelling works in the finite element by incorporating the sequence of mechanised tunnel construction. This sequence includes different stages such as excavation, machine advancement, lining installation, grouting, and hardening. Additionally, most previous research, especially empirical formulas, is solely applicable for a single tunnel. Many existing empirical formulas are tailored to single tunnel scenarios and do not account for modern tunnelling machine processes. This work extends the analysis to twin tunnels constructed using mechanised tunnelling, which presents a more complex scenario due to the interactions between adjacent tunnels. The utilisation of the 3D numerical modelling approach provides a sustainable solution for tailoring mitigation strategies to effectively manage SS while accounting for project-specific parameters such as tunnel geometry, geotechnical engineering factors, and tunnelling operational variables. The 3D numerical model approach empowers the customization of SS mitigation strategies. This customization ensures not only the effectiveness of the mitigation strategies but also the adherence to sustainability principles by minimizing adverse impacts on the environment and adjacent structures. By identifying potential SS issues at various construction stages through early estimation, timely adjustments and proactive measures can be implemented, thereby mitigating risks, and reducing the likelihood of detrimental effects on the surroundings.

2. Field Measurement and Input Parameter for Numerical Analysis

The Mass Rapid Transit (MRT) project is located in the Klang Valley, Malaysia, which comprises of total distance track of 52.2 km including the elevated works, underground works and one deport of the MRT system with the purpose of public transport for the citizens. In this project, the total distance of the twin tunnels is 13.5 km. The twin tunnel is two tunnels that can be arranged into four configurations: parallel, stacked, perpendicular crossing, and offset arrangements. The tracked alignment, a combination of parallel and stacked underground tunnels, is from Jalan Ipoh to Bandar Malaysia. Only parallel tunnelling is used for this study, as limited information is only available for stacked underground tunnels. The tunnel is traversed through limestone, granite, alluvium, and the Kenny Hill formation, as shown in Figure 1.

The project used mechanised tunnelling for tunnel excavation, including EPB and variable density tunnelling (VDT) machines with 6.7 m diameter for different ground conditions. VDT is the new mechanised tunnelling machine that was developed to mitigate the risk of blowouts and sinkholes in variable and mixed ground conditions. Whereas EPB was selected for the homogenous and cohesive soils. It uses the excavated material to stabilise the tunnel face with constant pressure. The selected section for the 3D numerical analysis is located between the alluvium and Kenny Hill. Several early geologists [29,30] have used diluvium as a term to describe “alluvium-like” materials comprised of sands, gravels, boulders, and clays on the hard rock in all landscape positions. The alluvium is modern sediments deposited into the ground as a result of recent river activity [31]. Since this soil material is deposited by the river, it can be anticipated that the ground is soft for the alluvium. Furthermore, the Kenny Hill formation in Kuala Lumpur has a thickness of 1200–1500 m of clastic sedimentary rocks and is located in the west and south [32]. It can be stated that this formation is for weathered rock and residual soil [33]; therefore, this formation is anticipated as a firm and hard material. This project encounters two different soil hardness materials: soft (alluvium) and hard (Kenny Hill).

Figure 1. The geological formation and location of the twin tunnels investigated in this research [34].

Figure 1. The geological formation and location of the twin tunnels investigated in this research [34].

2.1. Field Measurement of SS

Throughout the chainage of the tunnelling, the overburden of the soil at the tunnel crown is in the range of 6 m to 42 m. For the back analysis of the SS due to tunnelling, five cross sections of the tunnel that traverse through the soil with a sufficient array of ground settlement markers were selected. In this project, the SS values were measured using a total station from the ground settlement marker as shown in Figure 2.

The ground settlement marker for this project is an instrument with a steel rod length of 0.95 m embedded in the soil that is covered by a steel plate of 150 mm by 150 mm, which acts as a marker label to ease the identification of the location of the markers. The protrusion of a 0.05 m steel rod above the marker label was used to measure the level of the change. The reduced levels of the tops of the rods of settlement markers based on the nearest benchmark (BM) or temporary BM were measured prior to the commencement of the tunnelling works. The baseline readings were set after taking readings for several days before the commencement of the tunnelling works. The settlement was calculated by taking the change in the reduced level of the top of the rod relative to the base readings and previous readings. In general, the readings of the settlement markers were taken once per day during the construction period of the tunnelling. To evaluate the settlement from the field data, the settlement values were collected when both tunnels had passed by the settlement array, as illustrated in Figure 3, and only the final readings were collected as the SS. There are approximately 304 subsurface investigation (SI) boreholes along the area of the tunnelling chainage. In the numerical analysis, the nearest sections to the SI boreholes were utilised to establish the geological profile of the model.

2.2. Interpretation of Geotechnical Parameters

To carry out the numerical analysis of any geotechnical application, the accuracy of the soil-structure interaction problems is crucial as the soil is the weakest material compared to others such as concrete and steel. Hence, the analysis is strongly dependent on the soil constitutive model to reflect the actual condition of the ground.

The soil profile for each model comprises two layers of soil and rock. The tunnels are bored through the soil layers (alluvium and Kenny Hill). Hence, the rock layer does not have an impact on the analysis. The range of soil stresses and strains are important in the analysis; hence, this shall be considered in the numerical analysis. For the finite element analysis, two common soil constitutive models were used by many researchers in relevant studies, which are Mohr-Coulomb (MC) [25,35] and Hardening Soil (HS) [36,37]. Gerhard [37] compared the MC and HS of the tunnelling finite element analysis and found that the MC model can provide higher settlement values in comparison with the measurements on the site, while the HS model predicts SS values that are closer to the actual measurements using the 2D grout pressure method. In addition, Hejazi et al. [38] found the SS induced by tunnelling is strongly influenced by the soil constitutive model. From their findings, MC SS appears unrealistic, which is due to the stiffness modulus in the MC model being constant at certain points of strain. Islam and Iskander [28] mentioned that MC is not able to stimulate the behaviour of the soil, especially during unloading and tunnel excavation. However, they stayed with the MC model due to the simplicity of the parameters. Moller [39] carried out both analyses using 2D and 3D finite element analysis for the second Heinenoord slurry shield using the HS model, and the findings showed the simulated SS is close to the site’s actual settlement values. His analysis results are aligned with another study conducted by Likitlersuang et al. [36], where they reported very close readings to the measured settlement using the HS model. Choon [40] emphasised that the soil constitutive models play a crucial role in the analysis results and concluded that tunnelling involves unloading and reloading the surrounding soils. Hence, the HS model was recommended by Choon [40] for use in the assessment of SS. In addition, Janin et al. [41] stated that the HS model can distinguish the stiffness modulus for the primary loading and unloading, which makes this numerical analysis closer to the measured SS results.

The soil properties can be classified by the strain levels into three categories, which are very small strain level (stiffness modulus is constant in elastic range), small strain level (stiffness modulus varies non-linearly with the strain), and large strain level (soil is close to the failure, and the soil stiffness is relatively small) where tunnelling works can be considered as large strain level [42]. The illustration of the categories of the strain level can be seen in Figure 4a. The HS model was created with the concept of plasticity, where the total strains of the soil are calculated using a stress-dependent stiffness that is different for both loading and unloading. The MC model indicates that the soil behaviour is linearly elastic and perfectly plastic. In addition, the HS model can consider loading and unloading stages, which can reflect the actual situation of tunnel construction for the removal of soil (unloading), followed by tunnel lining construction (loading). Therefore, the HS model is a more suitable soil model for the tunnel numerical analysis of tunnelling compared to MC [43]. An illustration of an indication of the difference between the MC and HS is shown in Figure 4b.

The soil parameter of effective strength is retrieved by carrying out consolidated isotopically undrained (CIU) tests from the closest boreholes at different depths. In addition, several Menard Pressuremetre Tests (PMTs) have been carried out as per BS 5930 [44] throughout the chainage of the tunnelling. The PMT was placed inside the borehole after the Standard Penetration Test (SPT-N) was taken. The depth of PMT was located at approximately the same depth as SPT-N, as this was conducted to ensure the data from SPT-N could be used as a reference for the comparison of similar lithologies. The corrected SPT-N was correlated with the pressuremetre modulus, E_m. The previous researchers [45,46] carried out the correlation of the corrected SPT-N of N₆₀ with the elastic modulus from the PMT. Some of the E_m empirical formulas are presented in

Table 3. Summary of the Em empirical equations to determine the elastic modulus of the soil.

Author (s)	Ground Condition	Method of Tunnelling	R2
Ohya et al. [45]	$E_m=19.3{\left(N_{60}\right)}^{0.63}$	Clayey soil	0.39
Yagiz et al. [46]	$E_m=388.67\left(N_{60}\right)+4554$	Silty Clay	0.83
Cheshomi and Ghodrati [47]	$E_m=N_{60}-2.67$	Silty Clay	0.85
Naseem and Jamil [48]	$E_m=15.214N_{60}+89.276$	Sandy soils	0.88

according to different types of soil. The performance prediction of these formulas is presented based on the coefficient of determination (R²).

A total of 20 and 41 PMTs were carried out at the alluvium and Kenny Hill formations along the chainage of the tunnel, respectively. After carrying out the removal of outliers using the method of the interquartile range rule for a total of 61 PMTs, the corrected STP-N₆₀ is plotted against E_m for 2 different formations, as shown in Figure 5.

From the SPT-N, the selected section for the tunnelling analysis with a borehole is correlated with the elastic modulus according to the type of soil. However, $E_m$ is different from the soil elastic modulus; hence, conversion of the calculation is required. The difference could be due to several reasons, such as the difference in the range of radial, the stress surrounding the borehole wall, the possibility of disturbance during drilling and the installation of probes, and the computation of the modulus with the assumption of cylinder length as infinite [49]. Therefore, correction of the elastic modulus, $E$, of the soil is introduced to the actual value from $E_m$ [50], which is defined as follows:

$$\alpha =\frac{E_m}{E}$$

(5)

where α is a Menard’s factor. In the HS model, there are three crucial elastic modulus inputs for the analysis, which are namely: secant stiffness modulus corresponding to 50% of the ultimate deviatoric stress, $E_{50}$; Oedometer modulus, $E_{Oed}$; and unloading or reloading modulus, $E_{ur}$. Secant stiffness from the drained triaxial test at 50% of the ultimate deviatoric stress $E_{50}$ is the confining stress that depends on stiffness modulus for primary loading where the amount of stress is dependent on the power of $m$.

$$E_{50}=E_{50}^{ref}{\left(\frac{c^{\prime }\cos {\Phi }^{\prime }-{\sigma }_3^{\prime }\sin {\Phi }^{\prime }}{c^{\prime }\cos {\Phi }^{\prime }+P^{ref}\sin {\Phi }^{\prime }}\right)}^m$$

(6)

where p^ref is a reference stress of 100 kPa, c’ is effective cohesion, and ${\Phi }^{\prime }$ is effective friction angle. The stiffness is affected by minor effective principal stress. The power $m$ governs the stress dependency. As reported by Janbu [51], the sand content of $m$ is approximately 0.5. Other than this input, the reference oedometer modulus, $E_{oed}^{Ref}$, acts as the parameter that controls the magnitude of the plastic strains that arise from the yield cap ${\epsilon }_v^{pc}$. The oedometer modulus can be defined as follows:

$$E_{Oed}=E_{Oed}^{Ref}{\left(\frac{c^{\prime }\cos {\Phi }^{\prime }-{\sigma }_1^{\prime }\sin {\Phi }^{\prime }}{c^{\prime }\cos {\Phi }^{\prime }+P^{ref}\sin {\Phi }^{\prime }}\right)}^m$$

(7)

where $E_{oed}^{Ref}$ is approximately similar to $E_{50}^{Ref}$ [52] as illustrated in Figure 6. In addition, the stress-dependent stiffening modulus for unloading and reloading can be calculated as follows:

$$E_{ur}=E_{ur}^{Ref} {\left(\frac{c^{\prime }\cos {\Phi }^{\prime }-{\sigma }_3^{\prime }\sin {\Phi }^{\prime }}{c^{\prime }\cos {\Phi }^{\prime }+P^{ref}\sin {\Phi }^{\prime }}\right)}^m$$

(8)

For the ease of numerical computation, several researchers [40,53] have used $E_{50}^{Ref}$ = 3$E_{ur}^{Ref}$ as per recommended by the previous researchers [52,54]. The input of E₅₀ is correlated with the range from 0.06E to 0.26E with 0.5q_f (where E is the elastic modulus of the soil and q_f is the deviatoric stress from the secant stiffness modulus reduction curves from static torsional and triaxial shear data on clays and sands) [55]. Based on the interpretation,

Table 4. Summary of the stiffness and effective strength of the soil used in the numerical analysis.

Model	Effective Cohesion, c’ (kPa)	Friction Angle (°)	${\mathbf{E}}_{50}^{\mathbf{R}\mathbf{e}\mathbf{f}}$ (MPa)	${\mathbf{E}}_{\mathbf{o}\mathbf{e}\mathbf{d}}^{\mathbf{R}\mathbf{e}\mathbf{f}}$ (MPa)	${\mathbf{E}}_{\mathbf{u}\mathbf{r}}^{\mathbf{R}\mathbf{e}\mathbf{f}}$ (MPa)	Type of Formation
1	2	29	9.3	9.3	27.9	Alluvium
2	3	32	26.0	26.0	78.0	Alluvium
3	2	28	6.1	6.1	18.3	Alluvium
	5	32	11.6	11.6	34.8	Kenny Hill
4	2	30	7.0	7.0	21.0	Alluvium
5	1	25	3.5	3.5	10.5	Alluvium

summaries the geotechnical input parameters for these five models.

A total of 36 and 35 bulk density tests are carried out on the alluvium and Kenny Hill formations across the tunnelling alignment at various depths, respectively. The statistical input of bulk density is summarised in

Table 5. Summary of the inputs of the bulk density of different ground types.

Type of Ground	Bulk Density (kN.m−3)
	Mean	Maximum	Minimum	Standard Deviation
Alluvium	19	21	16	1.05
Kenny Hill	20	22	15	1.51

. In addition, the groundwater level recorded from the closest borehole with the readings recorded for morning and evening throughout the SI works will be used as part of the analysis, and these readings are taken as an average of the groundwater level. The groundwater level is one of the important parameters because when the groundwater level is lowered, it can cause a reduction in the pore pressure, which indirectly increases the surface settlement [56]. Figure 7 shows an illustration of how the lowering of groundwater causes the settlement.

2.3. Tunnel Boring Machine (TBM) Operational Parameters

Various operational factors of TBMs, including face pressures, rates of penetration, and pressures of grouting the tail void, jointly impact the surface settlement caused by mechanised tunnelling operations. TBMs cause surface settlement mainly by moving soil during excavation, resulting in empty spaces where the neighboring ground eventually subsides. While the TBM constructs the tunnel, it can also shift the ground, prompting settlement. The interplay between the tunnel lining and the nearby soil can trigger certain adaptations and shifts in the ground, indirectly resulting in surface settlement. Numerous studies [3,24] have demonstrated that face pressure and grouting pressure are the two important contributors to the SS caused by mechanised tunnelling. Hence, in this study, these two parameters were considered to be part of the finite element analysis. The EPB shield of the face pressure operates with the cutter head rotating, and cutting tools scrap the ground from the tunnel face while additives are injected into the material. The amount of the excavated soil is controlled and transported from the shield, and the tunnel face can be supported by the soil stored in the chamber. Hence, the face pressure stores in the chamber are crucial to maintaining the stability of the excavation. Different conditions of face pressure result in different SS scenarios. Several studies have found that increasing face pressure causes a decrease in settlement [3,57]. A different ideal scenario for the face pressure that affects the SS is illustrated in Figure 8. Case 1 shows that when the face pressures are equivalent to the overburden pressure, there are no impacts on the ground surface. Case 2 depicts the face pressures being less than the overburden pressures and the ground settling. Lastly, case 3 illustrates when the face pressures are more than the overburden pressures, heaving can be seen. Nevertheless, the amount of face pressure exceeded or lesser than the required pressures is subjected to the soil properties.

The shield diameter is larger than the external diameter of the installed lining, which causes a gap between the segments and the excavated soil mass [58]. This gap difference is known as a tail void or annulus, and it easily occupies the inward deforming soil during the shield and tail skin assembly process. When the tunnelling machine advances, the void can be continually filled with grout in one of two ways: by pumping the grout through pre-cast concrete lining segments or injecting the grout via grout ports at the back of the tail skin. The injection was carried out under high pressure in the tail void formed by the shield tail’s back. Thus, this is known as tail void grouting pressure. With high-pressure injection, the void can be filled within a short period of time with grouting material.

Other than these two parameters, another operating parameter that can be modelled into the numerical analysis is jack-force thrust. The distribution of the jack thrust in this study is uniform across the theoretical models. The function of this thrust force is to stimulate the actual condition of the mechanised tunnelling to move forward. The face pressure and grouting pressure, with the average input of the recorded face pressure and grouting pressure, are adopted in the analysis. To avoid the complication of the analysis, the value is rounded off to the nearest two significant values, and these numbers are summarised in

Table 6. Summary of the TBM operational parameters.

Model	Tunnel Bound	Average Face Pressure (kPa)	Average Grouting Pressure (kPa)	Average Force Thrust (kPa)
1	1st	140	230	4700
	2nd	180	370	2900
2	1st	290	350	4400
	2nd	230	390	3800
3	1st	190	170	5400
	2nd	200	160	5800
4	1st	190	220	5800
	2nd	210	250	5800
5	1st	220	410	4400
	2nd	250	270	5200

.

3. Construction of Numerical Model

3.1. General Characteristics

The numerical analysis was carried out using MIDAS GTX NX. In this section, the geometry dimension of the numerical model for the analysis will be further explained, the input for the grouting materials for the analysis will be stated, and the boundary condition of the model will be described.

Before the commencement modelling of the numerical analysis, the extended dimensions of the models affect the analysis required to be identified. In this analysis, the overall geometry of the numerical model was based on the recommendations by Alagha and Chapman [59] and Ruse [60]. Where C is the tunnel cover, D is the tunnel diameter, the recommendation for the minimum numerical model width is 2D, the depth of the model is C + 3D/2, and the length is 3D. Figure 9 shows a view of the models with various dimensions (x, y, and z directions). The overall adopted numerical model dimension is 115 m (width) $\times$ 115 m (depth) $\times$ 20 m (thickness). In this numerical model, the diameter of the tunnel is 6.7 m with a segmental lining width of 1.5 m. The tail void grout thickness is 150 mm. Two types of grouting were used in the analysis which are fresh and hardened grouting. The properties of the grouting materials are tabulated in

Table 7. The inputs of the fresh and hardened grouting.

Type of Material	Unit Weight, γ (kN.m−3)	Young’s Modulus, E (MPa)
Fresh grout	15	7.5
Hardened grout	15	15

. In addition, the degree of freedom for the numerical model was constrained in the x-direction, whereas in the y-direction was on both sides. Lastly, only nodes at the bottom were constrained in the z-direction. The top nodal along the ground surface was not constrained to allow the movement of the ground during the tunnelling process.

3.2. Construction Stages

Several stages have been involved in stimulating tunnel construction. In general, the process of tunnelling includes soil excavation, the installation of segment lining, and the application of external forces. At the beginning of the process, the self-weight and groundwater were activated to simulate the initial stage to establish the geostatic stresses in the soil mass. Subsequently, excavation was carried out by deactivating the soil elements with face pressure on the excavation face and, at the same time, activating the shield shell, which represents the tunnelling machine moving forward (Figure 10a). This process of excavation was repeated until the total length of the shield (equal to six segmental linings) was moved into the soil mass and the force thrust of tunnelling machine was activated to indicate the movement of tunnelling machine, and the segmental lining was installed at this stage (Figure 10b). Subsequently, the tail void grouting pressure was activated with fresh grout material (Figure 10c). Finally, the fresh grout was converted to hardened grout (Figure 10d) to show the later stage of the project. The tunnelling operational procedures of the numerical analysis for the tunnelling machine are shown in Figure 11. The different stages involved in tunnel construction, which encompasses excavation involving face pressure, lining installation, and grouting, exert an influence on SS. During excavation, voids are formed that cause the surrounding soil to settle. The installation of the lining can displace the ground, contributing to settlement as the tunnel’s supporting structure interacts with the soil. The grouting process, aimed at filling the gaps between the tunnel and the excavated area, affects the soil’s radial distribution, consequently impacting SS. Each of these phases brings about distinct settlement effects, and the overall sequence and coordination of these stages play a pivotal role in effectively managing settlement levels.

4. Results and Discussion

In this section, the output SS from the numerical analysis is compared with field measurements of SS, and the interpretation of the results will be further discussed. A total of five numerical models were analysed and compared with the measured settlement at the site. A total of five models’ geometries can be found in

Table 8. Summary of the tunnel geometry models for the numerical analysis.

Model	Tunnel Depth (m)	Pillar Width (m)	Location of Maximum SS
1	10	6	Between 1st and 2nd tunnel
2	13	10	1st bored tunnel
3	15	11	1st bored tunnel
4	15	19	1st bored tunnel
5	15	12	2nd bored tunnel

. One of the tunnel geometries that impacts the SS between twin tunnels is the pillar width. Based on Table 8, it can be observed that Model 1, with the shortest pillar width, exhibits the highest SS between the tunnels. This implies that both the 1st and 2nd bored tunnels contribute to the maximum SS occurring between them.

The field measurements of SS directly above each bored tunnel closely align with the results obtained from the numerical analysis, as presented in

Table 9. Summary of the SS results through numerical analysis.

Model	Sequence of Boring	SS Obtained by MIDAS (mm)	Actual SS (mm)	Percentage of Difference (%)
1	1st	9.8	9.3	5.4
	2nd	7.3	7.5	2.7
2	1st	9.2	9.3	1.1
	2nd	4.5	3.5	28.6
3	1st	13.6	11.6	17.1
	2nd	12.6	10.7	17.8
4	1st	17.2	18.3	6.0
	2nd	16.1	16.5	2.4
5	1st	17.5	18.5	5.4
	2nd	19.0	22.1	14.0

. A comparison between the numerical analysis and the actual measurements is conducted in terms of percentage differences. The disparities are predominantly below 18% across most models, with the exception of Model 2 of the second bored tunnel, where a notably higher percentage difference of 28.6% is observed. In comparison with the actual measurement, only a minor difference of 1 mm is observed between the actual value and the analysis. In general, the difference between measured and actual measurement is in the range of 0.5 mm to 3.1 mm.

Additionally, the SS values obtained from measurements and numerical analysis for both the first and second bored tunnels are graphically represented in Figure 11, with coefficient of determination (R²) values of 0.94 and 0.96 for the first and second bored tunnel SS, respectively. These high R² values suggest strong correlations between the measured and computed SS values for both tunnels.

The SS output of the numerical analysis for each model is visualised alongside the corresponding field measurement, as depicted in Figure 12. Observing the patterns within Figure 12, it can be seen that the 3D numerical analysis of the tunnel construction process effectively approximates the measured SS.

The output of this analysis has shown close similarities between the prediction settlement using MIDAS software and the site measurement. Model 1 has the shortest pillar width (6 m) in comparison with all the models. From this model, it can be seen that the 1st and 2nd bored tunnels do not have much difference (2 mm); however, the maximum SS is found located in between these two tunnels with the shape of a single symmetric settlement trough, which is similar to the findings as per Peck [5] and Likitlersuang et al. [36].

Model 2 depicts the unique trend of SS, where the maximum SS is found at the 1st bored tunnel instead of the 2nd bored tunnel. In the actual scenario of excavation for twin tunnels, the first bored tunnel excavates the undisturbed soil. However, the process of excavation for the first bored tunnel disturbed the soil surrounding it. Hence, this will cause the second tunnel’s soil to be disturbed. Terzaghi [61] was the first researcher who published a paper related to the field data of surface settlements above twin tunnels. His findings showed that the settlements above the second tunnel are larger than the first tunnel. However, based on the result of the analysis and the site-recorded SS, the second bore tunnel does not have the highest SS. In the past, most of the tunnelling operations for twin tunnels were carried out using compressed air [62]. In addition, several researchers [63,64] have found that the settlements above the second tunnels were larger than the first tunnels. The finding that the second bored tunnel has the maximum SS is applicable for open-face tunnelling because the ground reaction is influenced by tunnel geometry and geological conditions.

It is obvious that the SS values above the second tunnel could be lower. This is because the tunnelling operational parameters are considered contributing factors to the settlement [65]. Kannangara et al. [66] also found a similar trend in the outputs, where the observed second excavated tunnel has a lower SS when the face and grouting pressure are controlled at the same or higher margin compared with the first excavation. Although Model 3 has an approximately similar pillar width to Model 2 (different by 1 m) with a deeper depth, the 1st and 2nd bored tunnel SS do not have many differences. This could be due to the similarities in the operation parameters between the two tunnels and the deeper tunnel depth. Hence, these two bored tunnel SS of Model 3 have more consistent SS compared to Model 2.

Model 5 exhibits the highest settlement in comparison to the remaining models due to the lowest soil stiffness of the ground. It can be inferred that when twin tunnels have a larger spacing of 3D, there is no interaction between the tunnels [28]. This is evident from Model 4, where the maximum surface settlement at the 1st excavated tunnel is likely due to slightly low face pressure. Consequently, it can be deduced that within this range of pillar width, there is no impact from the first bored tunnel on the area of the second tunnel.

Several researchers have also carried out the analysis of tunnelling using MIDAS software and reported very close findings to the measured settlement for the different scenarios of the analysis as tabulated in

Table 10. Summary of tunnelling analysis using MIDAS software for various scenarios.

No	Author (s)	Number of Numerical Models	Comparison of the Numerical Analysis with Site Measurement
1	Zhou et al. [67]	1	The site measured maximum settlement = 9.72 mm at x = 12.25 m. The numerical result of maximum settlement = 9.94 mm at x = 11.25 m. The errors between the numerical results and measured data are small.
2	Ahmed et al. [3]	4	The maximum difference between the measured and numerical analysis is 25%.
3	Ayasrah et al. [68]	5	The highest longitudinal settlement is less than a 3 mm difference.

.

In previous studies, the most common method to determine the surface settlement induced by tunnelling is the empirical formula. This empirical formula is dependent on the VL, which varies according to the type of ground and installation tunnel. However, as far as the authors know, no literature on VL for the soil in Malaysia has been published. The use of empirical formulas is limited to certain ground conditions and installation methods. Chakeri and Ünver [69] compared the empirical formula with the measured settlement, and from their studies, the calculated and measured settlement diverge by more than two times. In addition, they also carried out the 3D numerical analysis; however, the 3D numerical analysis was carried out in the lithology of the ground for the Gungoren formation, which mainly consists of very stiff clay, hard clay, and dense sand. Furthermore, Choon [40] used the 2D finite element of the contraction method with various VL inputs to retrieve a similar SS as measured at the site. Furthermore, several studies [70,71,72] have utilised centrifuge laboratory modelling for sandy soil materials. However, these laboratory results have not been compared with field measurements, and they have not taken tunnelling operational parameters into consideration. Therefore, by utilising 3D numerical analysis, it becomes possible to integrate three crucial factors: tunnel geometry, geotechnical soil properties, and tunnelling operational parameters. This integration, combined with the incorporation of tunnelling machine construction processes, allows for the accurate determination of SS. This approach mitigates the constraints associated with alternative methodologies, as discussed previously. Consequently, this demonstrates that this method offers superior advantages compared to other approaches.

5. Limitations and Future Works

In this study, the number of field measurements is limited due to the focused investigation on the alluvium and Kenny Hill formation areas. Therefore, expanding the scope to include more field measurements from various locations, particularly within the soil area, with closer intervals would enhance the development and validation of the model. Additionally, if the budget allows, considering real-time monitoring for tunnelling-induced settlement measurements could be beneficial.

It is important to note that the geological profile is established based on the closest SI boreholes with a diameter of 100 mm, which implies a potential non-homogeneity in the studied geological profile. Furthermore, the input of groundwater levels relies on recorded SI works, which could be influenced by field conditions.

While 3D numerical modelling is time-consuming and unsuitable for real-time analysis during the operation of long tunnels, situations involving geological formations and tunnel operational parameters may arise. Therefore, new computational approaches such as machine learning and intelligent techniques could find applications in this context. Combining these approaches with established theories and empirical equations in relevant domains could yield a more generalised technique for predicting settlement induced by twin tunnels. Moreover, these techniques would offer practicality and relevance to engineers in the field.

6. Conclusions

The findings from the 3D numerical analysis of parallel tunnels using the MIDAS software, incorporating mechanised tunnelling construction stages and a hardening soil model, demonstrate an acceptable level of prediction accuracy for SS when compared with field measurements. Consequently, the industry may consider adopting this approach for SS assessment. Furthermore, the 3D numerical modelling approach enables the customization of mitigation strategies for SS based on specific project parameters, aiming to achieve a more sustainable tunnel construction. The key findings of this investigation can be summarised as follows:

(1)

The 3D numerical analysis produced SS above the tunnel crown of twin tunnels has shown R² = 0.94 and R² = 0.96, respectively for the first and second bored tunnels with the actual field measurements and the largest difference settlement of 3.1 mm.

(2)

The second bored tunnel does not consistently exhibit the largest SS because it can be influenced by other factors, such as tunnel geometry, geotechnical soil properties, and tunnel operational parameters. For example, a narrower pillar width can lead to higher SS between tunnels, as demonstrated in Model 1. Similarly, areas with lower soil stiffness could also result in elevated SS. Additionally, insufficient face pressure can contribute to increased settlement.

(3)

It can be stated that no effect of the first bored tunnel on the second bored tunnel area at a distance of equal or more than 3D between the tunnels has been observed.

(4)

The interpretation of the elastic modulus from the field pressuremetre test and SPT-N can be used as geotechnical soil stiffness parameter inputs in the numerical analysis. The interpreted relationship of the pressuremetre modulus with corrected SPT-N is as follows: $E_m=2.3 N_{60}+11$ (alluvium) and $E_m=1.4 N_{60}+18$ (kenny Hill).

(5)

Three main input parameters, namely tunnel geometry, engineering ground parameters, and tunnel operational parameters, considered in the numerical analysis yield results that closely align with the field measurements.