Numerical Model Validation for Detection of Surface Displacements over Twin Tunnels from Metro Line 1 in the Historical Area of Seville (Spain)

Citation: Bahri, M.; Mascort-Albea, E.J.; Romero-Hernández, R.; Koopialipoor, M.; Soriano-Cuesta, C.; Jaramillo-Morilla, A. Numerical Model Validation for Detection of Surface Displacements over Twin Tunnels from Metro Line 1 in the Historical Area of Seville (Spain). Symmetry 2022,14, 1263. https:// doi.org/10.3390/sym14061263 Academic Editors: Quan Jiang, Gennady Kolesnikov, Vladimir Siunev and Maria Zaitseva Received: 20 May 2022 Accepted: 15 June 2022 Published: 18 June 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). symmetry S S Article Numerical Model Validation for Detection of Surface Displacements over Twin Tunnels from Metro Line 1 in the Historical Area of Seville (Spain) Maziyar Bahri 1, Emilio J. Mascort-Albea 1,* , Rocío Romero-Hernández 1, Mohammadreza Koopialipoor 2, Cristina Soriano-Cuesta 1and Antonio Jaramillo-Morilla 1 1Departamento de Estructuras de la Edificación e Ingeniería del Terreno, Escuela Técnica Superior de Arquitectura, Instituto Universitario de Arquitectura y Ciencias de la Construcción, Universidad de Sevilla, Avenida de Reina Mercedes 2, 41012 Sevilla, Spain; [email protected] (M.B.); r[email protected] (R.R.-H.); [email protected] (C.S.-C.); [email protected] (A.J.-M.) 2 Department of Civil and Environmental Engineering, Amirkabir University of Technology, Tehran 15914, Iran; mr[email protected] *Correspondence: [email protected] Abstract: In order to solve connectivity problems in metropolitan areas, the development of underground metro lines constitutes an unquestionable requirement. However, the construction work thereof encounters unfavourable circumstances when surface excavations must be carried out that cross historical areas of the city, due to the need to control surface movements. The design of the metro in the city of Seville (Spain) from 2004 to 2006 provides a representative example of this situation and triggered major upheavals that exerted repercussions on historical buildings. For these reasons, the excavation stages of Line 1 of this metro have been simulated by numerical methods using FLAC3D software and validated with the results provided by the real conditions. Consequently, various surface settlements have been evaluated by taking not only variates of the main parameters that characterise the soil of Seville, but also of the various load situations and excavation conditions. Notable results have been achieved through calibration of 54 variants of the same model corresponding to Line 1, and their comparison with the real results obtained in nine critical areas of the itinerary. The results obtained have made it possible to determine the effects of excavation on the subsoil of the city of Seville with great accuracy, since the percentage error of calculated vertical surface movements varies from 0.1% to 5.3%. Keywords: soft soils; historical cities; critical areas; statistical analysis; FLAC3D; simulation of symmetrical excavations 1. Introduction Today, 55% of the world’s population lives in cities and highly densified metropolitan environments; an upward trend thereof has been predicted by major international organisations [ 1 ]. The obvious problems resulting from this population concentration materialise as traffic congestion, air pollution, and greenhouse gas emissions, which have led to a global climate emergency. For these reasons, the use of the subsoil holds the key to sustainability strategies for urban environments, through underground mobility [ 2 ] and the use of underground urban space [3,4]. Nevertheless, the construction of underground tunnels inherently results in varying degrees of ground movement towards excavation, thereby generating settlements of the ground surface [ 5 ]. This type of excavation can exert a detrimental effect on adjacent structures and pose environmental hazards if the produced horizontal and vertical movements are not controlled [ 6 ]. Such circumstances can be detected through the identification of critical areas in urban environments [7]. Symmetry 2022,14, 1263. https://doi.org/10.3390/sym14061263 https://www.mdpi.com/journal/symmetry Symmetry 2022,14, 1263 2 of 29 This problem is accentuated in historical cities due to the remarkable number of old buildings with highly vulnerable structures and construction systems [ 8 ]. In this type of context, the limitation of movements caused by underground excavations can become a key factor when operations are carried out in environments in close proximity to monumental complexes of great heritage value [ 9 , 10 ]. Consequently, there are numerous examples of recent research that have investigated historical structures that have been damaged during mechanised excavation [ 11 – 14 ]. These issues were revealed during the excavation work for Line 1 of the Seville metro (Spain) undertaken between 2004 and 2006. In this case, the excavated tunnels crossed the consolidated city through points of its main historical area, some of which were in close proximity to its UNESCO World Heritage Site. The damage caused by the movements produced was significant at certain critical points along the route [ 15 ]. This explains the need for strategies to be established to control the movements caused on the surface by means of predictive methods that enable the critical areas of future excavations in historical urban environments to be determined. Among the main calculation methods developed in recent years, the use of empirical methods, capable of generating simplified formulations, stands out. Additionally noteworthy is the use of numerical models, the different variants of which are currently widely used. In both cases, one of the most illustrative results is the calculation of the surface settlement curve, which can reflect the surface influence of the excavation [ 16 ]. Highly accurate values of the settlement curve can be obtained by using sophisticated numerical models that are capable of considering a notable set of parameters for their calculation. These models are based on numerical methods that have similarities and differences through a common mathematical approach [ 17 ]. These include the use of Finite Difference Method (FDM), Finite Element Method (FEM), and Finite Volume Method (FVM). Recent advances in computational technologies have made it possible to implement such methods for the modelling of complex geotechnical problems [ 18 ]. Consequently, through the development of different software, it is possible to determine the calculation of surface settlement curves through the application of FDM [ 19 ], FEM [ 20 ], and FVM [ 21 ] through different approaches and computer software. The accuracy of the results obtained shows that the simulated cases can be very specific, and that the range of parameters and criteria to be used is strongly conditioned by the conditions of each investigation. Consequently, one of the most important issues in the design of numerical calculation models is the selection of the parameters and calculation conditions utilised in obtaining the surface settlement curve. In this respect, the main factors to be considered are those related to the geometry of the tunnels, the excavation processes and the properties of the soil. The most influential tunnel design conditions in the simulation excavation process in numerical modeling are related to the geometry of the lining. For instance, curve lining has to deal with more stress and strain. Hence, the surface settlement thereof is greater than in other places. Regarding this issue, the analysis of the interaction and settlement in a curved section in mechanised excavation constitutes a highly relevant research topic [ 22 – 24 ]. Another main geometric factor for the calculation of tunnels is the dimension of tunnels. In these situations, conditions such as the depth, diameter, and distance of the tunnels constitute the essential parameters for this question in numerical modeling [ 25 , 26 ]. Initially, the excavation depth of tunnels is a decisive factor in the generation of surface settlement, and a threshold value of 16 m depth is determined for their mitigation [ 27 ]. Similarly, the existence of shallow tunnels is considered from excavation depth values of 10 m [ 28 ]. Concerning the diameter factor, Chakeri et al. [ 29 ] showed that mechanised excavation of a tunnel with a large diameter leads to increased maximum surface settlement. In this respect, Zheng et al. [ 30 ] studied the interaction between large-diameter tunnels and small-diameter tunnels, and illustrated the difference in bending moment and stress interaction of different diameters. Finally, the separation between tunnels should be indicated as a further factor to be considered. In this respect, Suwansawat and Einstein [ 31 ] studied the relationship between the maximum values of the surface settlements provided by the Gaussian curves Symmetry 2022,14, 1263 3 of 29 and the distance between tunnels. Similarly, recent research illustrates that the spacing of twin tunnels is also a relevant factor when the separation thereof is less than 20 m [ 32 – 34 ]. All mechanised excavation of TBM methods follows steps through a procedure that can be simulated in numerical modelling. In the estimation of the surface settlement induced by mechanised excavation using EPB techniques, the stages of the excavation should be simulated in a highly accurate way since the tunnelling-related risk scenarios. Additionally, the weight of EPB, and the properties of the lining are also implemented in the numerical models. By bearing all these factors in mind, high-precision surface settlement and interaction parameters can be determined over twin tunnels in shallow urban areas by simulation in the numerical models, through investigations that have been conducted to simulate EPB excavations and the interaction effect thereof in underground areas [ 35 – 41 ]. Face pressure and grout pressure are the main effective factors of mechanised excavation in EPB-TBM. In this regard, very recent studies illustrate the influence of these factors on the instability of the face of excavation and on tension interaction with soil [ 42 , 43 ]. Furthermore, Forsat et al. [ 32 ] showed that surface settlement has an indirect relationship with face pressure and grout pressure of EPB. There are also other factors in the excavation process that can be considered in the simulation results. Certain properties of the EPB machinery, such as the weight of its components and the lifting forces, have been included in research focused on the simulation of surface settlements [ 40 , 44 ]. Alongside excavation factors, the lining precast must be simulated precisely. In this case, the interaction between segments, the interface stiffness, and the axial joint of the lining have additionally been considered in underground excavation in recent research [45,46]. Soil properties and constitutive models of soil are important factors to consider in the numerical simulation of tunnel excavations. In this respect, Mohr–Coulomb models have been widely used in designing underground excavation in soil [ 8 , 18 , 41 ]. Several researchers have used a modified Cam–Clay model for the development of their numerical simulations [ 47 – 49 ], while others have analysed numerical models using the Plastic-Hardening (PH) model. This model is usually adopted as the constitutive model of soft clay and was established for tunnelling excavation problems [ 50 – 52 ]. Likewise, groundwater conditions must not be overlooked as a clear influencing factor. Shivaei et al. [ 53 ] studied both the drainage condition of the grout layer in mechanised excavation and the difference in surface settlement by changing the soil permeability and water table level. They concluded that with lowering the water table, the surface settlement decreased. Complementarily, Mu et al. [27] investigated the soil parameters and ground settlement during the excavation for the metro with EPB, and they declared that the rate of soil friction angle and soil cohesion have an indirect relationship with ground displacement. Finally, it must be mentioned that the interaction between tunnel lining and surface condition acts on the loading of the tunnel lining and ground displacement. One of the most evident factors involves the burden loads from buildings and traffic, which have been widely analysed [ 8 , 54 – 56 ]. Concerning this parameter, Chakeri et al. [ 29 ] divided the surface load into traffic load and building loads and applied all these loads as surcharges to the surface of the 3D numerical model, by considering 20 kN/m 2 and 30 kN/m 2 , respectively. Similarly, the geometric, structural, and constructional properties of the preexisting building that is located above ground can also be considered through the use of stiffness factors [57,58]. In this research, surface settlement changes during the range of excavation properties and soil properties were studied during the underground excavation by EPB in Seville. One of the innovations of this paper compared to previous research is studying the surface settlement based on the real soil and excavation properties of Seville. Thus, this research aims to evaluate the predictive potential of the results from numerical model calculations of surface settlement in historic urban environments caused by EPB-TBM twin tunnel excavation processes. In this way, and in anticipation of the future execution of Lines 2, 3, and 4 of the Seville metro, which will affect the city even more extensively, this research raises awareness for the need to develop predictive strategies that will enable detection of Symmetry 2022,14, 1263 4 of 29 those more critical areas. For this reason, a numerical model based on the monitoring data gathered during the construction of Seville metro Line 1 has been designed to simulate the executed excavations. This paper is therefore organised as follows: an introduction regarding the selected parameters and the characteristics of the case study (Section 2) are first given, followed by an explanation of the processes developed for the simulation of the numerical model (Section 3), and a discussion of the results attained regarding the validation of the model and the statistical potential of the data obtained for the detection of vulnerability zones in historical cities (Section 4). Conclusions are then drawn, with a view to the future lines of research generated by this work (Section 5). 2. Case Study: Critical Points from Line 1 of the Metro of Seville The city of Seville, with a metropolitan population of approximately 1,500,000 inhabitants and an area of 4905.04 km 2 [ 59 ], is the fourth largest metropolitan area in Spain and the capital of the autonomous region of Andalusia. Historically conditioned by the Guadalquivir River, Seville contains one of the most extensive and valued historical sites in Europe and hosts a UNESCO World Heritage Site in its urban centre [ 60 ], as well as numerous heritage buildings and archaeological remains. In this respect, its metropolitan conditions imply the need to improve its territorial connectivity through the development of an underground connection network project [ 61 ]. Nevertheless, since the second half of the 20th century, the development of this work has experienced numerous interruptions and delays due to a variety of causes, including the historical nature of the city and the conditions of its subsoil [62]. 2.1. Geotechnical Site Characterisation The metropolitan territory of Seville is located on soft alluvial soils found in the valley of the Guadalquivir River. In this respect, the city of Seville has a flat topography and is situated at an average height of 7 m above sea level. Local and regional geotechnical studies agree on the determination of a typical stratigraphic profile to define the composition of the subsoil in the study area [ 63 , 64 ], and this approach was assumed in the design models used in the construction of Line 1 of the Seville metro [65]. These studies also indicate that there is marked variability in the ranges of values that characterise the defined strata, and assume that the sequence of the layers indicated undergoes variations as a consequence of the different modifications made to the course of the Guadalquivir River as it passed through the city of Seville. By means of the above assumptions and based on the geotechnical studies derived from the different project phases of the Seville underground network, it can be affirmed that the stratigraphy employed in the development of this research shows the mechanical and resistant parameters as indicated in Table 1. Table 1. Ranges of soil properties along the route of Line 1 of the Metro in Seville. Soil Unit Density (kg/m3)Young’s Modulus (MPa) Poisson’s Ratio Cohesion (kPa) Friction Angle (Degrees) Permeability (cm/s) Min. Max. Min. Max. Min. Max. Min. Max. Min. Max. Min. Max. Filling 2.10 1.80 7.0 12.6 0.2 0.6 2.0 6.0 24.0 26.0 - Clay 2.20 2.00 29.0 45.0 0.2 0.4 4.0 10.0 30.0 40.0 2×10−92×10−7 Sand 2.20 1.80 20.0 29.0 0.2 0.4 17.0 50.0 26.0 29.0 2×10−52×10−2 Gravel 2.20 2.00 40.0 70.0 0.2 0.4 0.0 0.0 35.0 39.0 5×10−65×10−2 Marl 2.2 1.90 70.0 100.0 0.2 0.4 32.0 44.0 26.0 28.0 2×10−93×10−9 In this respect, it is worth highlighting the great importance of the thickness of the layer of anthropic fills in the main areas of the historical city, which can record building remains more than two millennia old. This layer is composed of heterogeneous materials whose behaviour is difficult to control. In the same way, it is common for the strata relating to clays and sands to appear intermixed, and with usually medium and low values of Symmetry 2022,14, 1263 5 of 29 resistance. The resistance increases with the presence of the gravel strata at a greater depth, and additionally, the presence of loams of high resistance and very low permeability constitute the geological substratum that forms the base of these alluvial soils. Finally, it should be pointed out that the proximity to the river, the low average height and the existence of numerous aquifers mean that the presence of groundwater is usually in very close proximity to the surface. 2.2. Seville Metro Line 1 Main Difficulties Based on the aforementioned conditions, the general design of the Seville metro network was drafted at the beginning of the 21st century [ 66 ]. The complete proposal comprised a total of four metro lines, the layout of which is organised as follows: Lines 1 and 2 cross the metropolitan area in an east-west direction, by going to the Guadalquivir River, with the particularity that Line 2 crosses the central area of the city of Seville completely and Line 1 is displaced to the south. In a complementary manner, Line 3 has a north-south orientation and Line 4 has an annular and closed geometry (Figure 1). Twenty years later, the implementation of the project is still incomplete and only Line 1 has been executed (years 2004 and 2006). This experience has illustrated the main difficulties faced by the project and the challenge it still poses today. Symmetry 2022, 14, x FOR PEER REVIEW 5 of 29 In this respect, it is worth highlighting the great importance of the thickness of the layer of anthropic fills in the main areas of the historical city, which can record building remains more than two millennia old. This layer is composed of heterogeneous materials whose behaviour is difficult to control. In the same way, it is common for the strata relating to clays and sands to appear intermixed, and with usually medium and low values of resistance. The resistance increases with the presence of the gravel strata at a greater depth, and additionally, the presence of loams of high resistance and very low permeability constitute the geological substratum that forms the base of these alluvial soils. Finally, it should be pointed out that the proximity to the river, the low average height and the existence of numerous aquifers mean that the presence of groundwater is usually in very close proximity to the surface. 2.2. Seville Metro Line 1 Main Difficulties Based on the aforementioned conditions, the general design of the Seville metro network was drafted at the beginning of the 21st century [66]. The complete proposal comprised a total of four metro lines, the layout of which is organised as follows: Lines 1 and 2 cross the metropolitan area in an east-west direction, by going to the Guadalquivir River, with the particularity that Line 2 crosses the central area of the city of Seville completely and Line 1 is displaced to the south. In a complementary manner, Line 3 has a north-south orientation and Line 4 has an annular and closed geometry (Figure 1). Twenty years later, the implementation of the project is still incomplete and only Line 1 has been executed (years 2004 and 2006). This experience has illustrated the main difficulties faced by the project and the challenge it still poses today. Figure 1. UNESCO World Heritage Sites and planned network lines for the Metro of Seville. Figure 1. UNESCO World Heritage Sites and planned network lines for the Metro of Seville. The twin tunnels of Seville metro Line 1 have a diameter of 6.35 m and a separation distance of 6.00 m, with a longitudinal development that allows a maximum slope of 4% in the steepest sectors. Based on these geometrical conditions, the tunnels were excavated to a shallow depth using an EPB-TBM LOVAT machine with a maximum thrust of 40,000 kN Symmetry 2022,14, 1263 6 of 29 and a maximum torque of 600 kN × m. During this process, the lining of the tunnel was reinforced with precast concrete, and each circular ring of the tunnel lining comprised six quadrilaterals with a thickness of 0.25 m and a length of 1.6 m [ 67 ]. Additionally, the most peripheral sections of the route were designed to lie above ground and the excavation of the various stations was carried out using reinforced concrete screens [ 65 ]. In this respect, the main adversities encountered during the excavation work were related to the underground crossing of the river and the main civil infrastructures, the interruptions in excavation, the influence of groundwater, the difficulties caused by the existence of archaeological remains in the areas of the historical city, and the many movements caused in numerous buildings in close proximity to the excavation route [15]. 2.3. Selection of Critical Areas The aforementioned circumstances made it possible to identify the main problems to be solved for the future construction of Lines 2, 3, and 4 of the Seville metro. In this way, a set of critical areas belonging to the completed Line 1 has been established that are representative of the problems addressed in this research (Figure 2). Symmetry 2022, 14, x FOR PEER REVIEW 6 of 29 The twin tunnels of Seville metro Line 1 have a diameter of 6.35 metres and a separation distance of 6.00 metres, with a longitudinal development that allows a maximum slope of 4% in the steepest sectors. Based on these geometrical conditions, the tunnels were excavated to a shallow depth using an EPB-TBM LOVAT machine with a maximum thrust of 40,000 kN and a maximum torque of 600 kN × m. During this process, the lining of the tunnel was reinforced with precast concrete, and each circular ring of the tunnel lining comprised six quadrilaterals with a thickness of 0.25 metres and a length of 1.6 metres [67]. Additionally, the most peripheral sections of the route were designed to lie above ground and the excavation of the various stations was carried out using reinforced concrete screens [65]. In this respect, the main adversities encountered during the excavation work were related to the underground crossing of the river and the main civil infrastructures, the interruptions in excavation, the influence of groundwater, the difficulties caused by the existence of archaeological remains in the areas of the historical city, and the many movements caused in numerous buildings in close proximity to the excavation route [15]. 2.3. Selection of Critical Areas The aforementioned circumstances made it possible to identify the main problems to be solved for the future construction of Lines 2, 3, and 4 of the Seville metro. In this way, a set of critical areas belonging to the completed Line 1 has been established that are representative of the problems addressed in this research (Figure 2). (a) (b) Figure 2. (a) Restoration works in sector C–D during 2004 (b) Critical areas and excavation directions of the route of Line 1 (Critical areas are annotated as A–F). An analysis of critical areas has long been recognised as a necessary step in the study of ground displacement. Six zones have been selected as the starting point for this research, and corresponding numerical models were developed and validated in these locations using real data. These areas are in densely populated areas adjacent to historical locations or metro stations where real-time ground displacement was monitored, and these critical points demonstrate the validity of the numerical model. Nervión (A), San Bernardo (B), San Sebastián (C), Puerta de Jerez (D), Plaza de Cuba (E), and Parque de los Príncipes (F) are the selected critical areas. Moreover, it should be borne in mind that the Seville metro tunnel was excavated through different soil layers. The tunnel’s depth and the type of soil in which it is excavated are listed in Figure 3b. The thickness of the soil layer varies in each critical area, Figure 2. ( a ) Restoration works in sector C–D during 2004 ( b ) Critical areas and excavation directions of the route of Line 1 (Critical areas are annotated as A–F). An analysis of critical areas has long been recognised as a necessary step in the study of ground displacement. Six zones have been selected as the starting point for this research, and corresponding numerical models were developed and validated in these locations using real data. These areas are in densely populated areas adjacent to historical locations or metro stations where real-time ground displacement was monitored, and these critical points demonstrate the validity of the numerical model. Nervión (A), San Bernardo (B), San Sebastián (C), Puerta de Jerez (D), Plaza de Cuba (E), and Parque de los Príncipes (F) are the selected critical areas. Moreover, it should be borne in mind that the Seville metro tunnel was excavated through different soil layers. The tunnel’s depth and the type of soil in which it is excavated are listed in Figure 3b. The thickness of the soil layer varies in each critical area, which has been considered in the various models. Three points were chosen on the surface to be monitored; two were placed along the tunnel’s central axis, and the third was placed at the centre point between the tunnels showing a symmetrical geometry (Figure 3a). Symmetry 2022,14, 1263 7 of 29 Symmetry 2022, 14, x FOR PEER REVIEW 7 of 29 which has been considered in the various models. Three points were chosen on the surface to be monitored; two were placed along the tunnel’s central axis, and the third was placed at the centre point between the tunnels showing a symmetrical geometry (Figure 3a). (a) (b) Figure 3. (a) Monitored surface points. (b) Soil layer distributions for the critical areas chosen. Finally, it has been considered that several previous studies have employed the Mohr–Coulomb and Cam–Clay constitutive models to describe the soil layer [24,46,68– 71]. Despite the Hardening Soil model (HS) being recommended for the prediction of the behaviour of soils during geotechnical excavations [8,72], the Mohr–Coulomb constitutive model was employed in this research to simplify the numerical calculation. Nevertheless, the precast lining was assumed to have an elastic constitutive model. The properties of the soil layers, the segmental lining, and the EPB characteristics are summarised in Tables 2 and 3. Table 2. Selected soil properties for the research critical areas. Soil Unit Density (kg/m3) Young’s Modulus (MPa) Poisson’s Ratio Cohesion (kPa) Friction Angle (Degrees) Permeability (cm/s) Filling 2000 9.8 0.30 4.90 25 - Clay 2100 29.4 0.30 9.81 32 1.01 × 10−8 Sand 2000 24.5 0.35 19.60 28 9.5 × 10−3 Gravel 2100 58.8 0.30 0.00 37 9.09 × 10−4 Marl 2005 83.3 0.35 39.2 28 2.83 × 10−9 Table 3. EPB’s shield, tailing of the EPB properties. Material Density (kg/m3) Young’s Modulus (MPa) Poisson’s Ratio Thickness (m) External Diameter (m) Length (m) EPB’s shield 7840 2 × 105 0.20 0.10 6.3 10.0 The precast lining 2500 29 × 103 0.20 0.25 6.0 1.6 Injected grout 1800 - - 28.00 1.7 - Hardened grout 1800 - - 103.00 71.0 - Figure 3. (a) Monitored surface points. (b) Soil layer distributions for the critical areas chosen. Finally, it has been considered that several previous studies have employed the Mohr– Coulomb and Cam–Clay constitutive models to describe the soil layer [ 24 , 46 , 68 – 71 ]. Despite the Hardening Soil model (HS) being recommended for the prediction of the behaviour of soils during geotechnical excavations [ 8 , 72 ], the Mohr–Coulomb constitutive model was employed in this research to simplify the numerical calculation. Nevertheless, the precast lining was assumed to have an elastic constitutive model. The properties of the soil layers, the segmental lining, and the EPB characteristics are summarised in Tables 2and 3. Table 2. Selected soil properties for the research critical areas. Soil Unit Density (kg/m3) Young’s Modulus (MPa) Poisson’s Ratio Cohesion (kPa) Friction Angle (Degrees) Permeability (cm/s) Filling 2000 9.8 0.30 4.90 25 - Clay 2100 29.4 0.30 9.81 32 1.01 ×10−8 Sand 2000 24.5 0.35 19.60 28 9.5 ×10−3 Gravel 2100 58.8 0.30 0.00 37 9.09 ×10−4 Marl 2005 83.3 0.35 39.2 28 2.83 ×10−9 Table 3. EPB’s shield, tailing of the EPB properties. Material Density (kg/m3) Young’s Modulus (MPa) Poisson’s Ratio Thickness (m) External Diameter (m) Length (m) EPB’s shield 7840 2×1050.20 0.10 6.3 10.0 The precast lining 2500 29 ×1030.20 0.25 6.0 1.6 Injected grout 1800 - - 28.00 1.7 - Hardened grout 1800 - - 103.00 71.0 - 2.4. Dimensions of the Model and Boundary Conditions The model’s dimensions were critical in this study, since they needed to be sufficiently large to avoid boundary effects in numerical simulations. Previous researchers have employed the following values to estimate the model’s minimum dimensions: H + 4D for the height of the model, H + 3D for its length, and 6H for the width of the model [ 68 ], where H denotes the depth of the model’s axis, and D denotes its diameter. FLAC3D was utilised to create a complete three-dimensional layout of the twin tunnels of the first line of the Seville metro, which measures 55 m in height, 86 m in width, and 40 m in length. The 3D finite-volume model of the tunnel and excavation area was designed with height Z, Symmetry 2022,14, 1263 8 of 29 length Y, and width X. Additionally, the tunnel was excavated in the direction of Y in the model’s centre (Figure 4). Symmetry 2022, 14, x FOR PEER REVIEW 8 of 29 2.4. Dimensions of the Model and Boundary Conditions The model’s dimensions were critical in this study, since they needed to be sufficiently large to avoid boundary effects in numerical simulations. Previous researchers have employed the following values to estimate the model’s minimum dimensions: H +4D for the height of the model, H +3D for its length, and 6H for the width of the model [68], where H denotes the depth of the model’s axis, and D denotes its diameter. FLAC3D was utilised to create a complete three-dimensional layout of the twin tunnels of the first line of the Seville metro, which measures 55 metres in height, 86 metres in width, and 40 metres in length. The 3D finite-volume model of the tunnel and excavation area was designed with height Z, length Y, and width X. Additionally, the tunnel was excavated in the direction of Y in the model’s centre (Figure 4). Figure 4. Main dimension parameters of the simulated models, representing twin tunnelling. Each of the six geometries of underground simulation has a different number of nodes. Nevertheless, the size ratio of these meshes throughout the geometry remains the same. The model of San Bernardo (B), with 260,640 zone numbers and 271,350 grid points, covers most of the zone and grid points. On the other hand, the model of Plaza de Cuba (E) and Parque de los Príncipes (F) has fewer zone and grid points, with 242,080 zone numbers and 252,234 grid points, respectively. The surface of the model is left free, thereby allowing for the measurement of vertical movement on the surface, while considering the boundary condition. The Y-Z-plane boundary nodes (X = −49 and X = 37) and the X-Z-plane boundary nodes (Y = 0 and Y = 40) are fixed in the normal direction, while the lower boundary is fixed in the vertical direction. It should be noted that the mesh density was increased around the tunnel and decreased as one moved away from the tunnel. This would improve the model’s accuracy and speed [40,46]. The 1 centres of the first and second tunnels are located at X = 0 and X = −12, respectively. Therefore, X = −6 is located in the twin tunnels’ centre point: these points are highlighted in Figure 2. Study points are assigned to the surface of each geometry at X = 0, −6, and −12 to investigate surface settlement in this study. On the surface of the numerical model, the uniform traffic load (live load) and the building load (dead load) are considered. Additionally, in each model, a different load is applied (20,70, and 130 kPa). Thus, 130 kPa is assumed in the densely populated area, while 70 kPa is assumed in the normally populated area of Seville, and 20 kPa is recommended for locations with the lowest rate of buildings and traffic. To this end, it has been assumed that this method simplifies the creation of a model [53]. Two distinct methods are considered when determining the weight of the buildings [29,57]. Moreover, pore pressure is used to simulate the water condition in excavation simulations. Figure 4. Main dimension parameters of the simulated models, representing twin tunnelling. Each of the six geometries of underground simulation has a different number of nodes. Nevertheless, the size ratio of these meshes throughout the geometry remains the same. The model of San Bernardo (B), with 260,640 zone numbers and 271,350 grid points, covers most of the zone and grid points. On the other hand, the model of Plaza de Cuba (E) and Parque de los Príncipes (F) has fewer zone and grid points, with 242,080 zone numbers and 252,234 grid points, respectively. The surface of the model is left free, thereby allowing for the measurement of vertical movement on the surface, while considering the boundary condition. The Y-Z-plane boundary nodes (X = − 49 and X = 37) and the X-Z-plane boundary nodes (Y = 0 and Y = 40 ) are fixed in the normal direction, while the lower boundary is fixed in the vertical direction. It should be noted that the mesh density was increased around the tunnel and decreased as one moved away from the tunnel. This would improve the model’s accuracy and speed [40,46]. The 1 centres of the first and second tunnels are located at X = 0 and X = − 12, respectively. Therefore, X = − 6 is located in the twin tunnels’ centre point: these points are highlighted in Figure 2. Study points are assigned to the surface of each geometry at X = 0, − 6, and − 12 to investigate surface settlement in this study. On the surface of the numerical model, the uniform traffic load (live load) and the building load (dead load) are considered. Additionally, in each model, a different load is applied (20, 70, and 130 kPa ). Thus, 130 kPa is assumed in the densely populated area, while 70 kPa is assumed in the normally populated area of Seville, and 20 kPa is recommended for locations with the lowest rate of buildings and traffic. To this end, it has been assumed that this method simplifies the creation of a model [ 53 ]. Two distinct methods are considered when determining the weight of the buildings [ 29 , 57 ]. Moreover, pore pressure is used to simulate the water condition in excavation simulations. 3. Simulation Processes Mechanised tunnelling in shallow soil requires careful excavation in stages. This process entails the penetration of the EPB Shield, by applying face pressure to the soil and the injection of grout, and by installing the precast lining. These factors must be considered when simulating mechanised excavation in FLAC3D. Multiple researchers have incorporated these variables into their numerical model [ 73 – 77 ]. The progression of EPB advancement in Seville has been modelled in this study, and an attempt has been Symmetry 2022,14, 1263 9 of 29 made to model it using significant factors from the step-by-step excavation. The following summarises the simulation of the mechanised excavation in Seville. The metro system in Seville consists of twin tunnels, each of which is excavated in 50 steps, giving a total of 100 steps. When twin tunnels were excavated with two EPB machines in the same direction but with a delay, the distance between the shields of the two EPB machines was critical in the simulation [ 78 ]. In contrast, the twin tunnels in the Seville metro were excavated by a single EPB machine in a time-delayed manner. Using this method, one of the tunnels (Tunnel A) was excavated by an EPB machine, and the EPB machine returned in the opposite direction of Tunnel A, subsequently excavating Tunnel B. Thus, it is clear that the surface and ground had time to rebalance between the excavations of Tunnel A and Tunnel B, since the time interval between the two excavations exceeded one whole year, as determined by monitor data. In the end, two tunnels were excavated symmetrically. The twin tunnels were excavated in the Y-axial direction, 6 m apart. The EPB shield penetrated the soil in the first step of the excavation, and face pressure was applied (Figure 5a). Since the length of the EPB machine is 15 m following this portion of the excavation (Figure 5b), soft grout was added behind the EPB shield (Figure 5c). The precast lining was subsequently installed with a length of 1.6 m, as shown in Figure 5d. After one further step of excavation, the grouting pressure was removed, and the grout began to dry. This cycle was repeated until the EPB machine could be removed from the geometry and the tunnel lining was completely implanted (Figure 5e,f). Face pressure, grout pressure, the lining, and the EPB shield are illustrated in Figure 6, whereby dried grout and injected grout are highlighted in red and yellow, respectively. Moreover, the EPB shield is drawn with a blue line. Symmetry 2022, 14, x FOR PEER REVIEW 9 of 29 3. Simulation Processes Mechanised tunnelling in shallow soil requires careful excavation in stages. This process entails the penetration of the EPB Shield, by applying face pressure to the soil and the injection of grout, and by installing the precast lining. These factors must be considered when simulating mechanised excavation in FLAC3D. Multiple researchers have incorporated these variables into their numerical model [73–77]. The progression of EPB advancement in Seville has been modelled in this study, and an attempt has been made to model it using significant factors from the step-by-step excavation. The following summarises the simulation of the mechanised excavation in Seville. The metro system in Seville consists of twin tunnels, each of which is excavated in 50 steps, giving a total of 100 steps. When twin tunnels were excavated with two EPB machines in the same direction but with a delay, the distance between the shields of the two EPB machines was critical in the simulation [78]. In contrast, the twin tunnels in the Seville metro were excavated by a single EPB machine in a time-delayed manner. Using this method, one of the tunnels (Tunnel A) was excavated by an EPB machine, and the EPB machine returned in the opposite direction of Tunnel A, subsequently excavating Tunnel B. Thus, it is clear that the surface and ground had time to rebalance between the excavations of Tunnel A and Tunnel B, since the time interval between the two excavations exceeded one whole year, as determined by monitor data. In the end, two tunnels were excavated symmetrically. The twin tunnels were excavated in the Y-axial direction, 6 metres apart. The EPB shield penetrated the soil in the first step of the excavation, and face pressure was applied (Figure 5a). Since the length of the EPB machine is 15 m following this portion of the excavation (Figure 5b), soft grout was added behind the EPB shield (Figure 5c). The precast lining was subsequently installed with a length of 1.6 m, as shown in Figure 5d. After one further step of excavation, the grouting pressure was removed, and the grout began to dry. This cycle was repeated until the EPB machine could be removed from the geometry and the tunnel lining was completely implanted (Figure 5e,f). Face pressure, grout pressure, the lining, and the EPB shield are illustrated in Figure 6, whereby dried grout and injected grout are highlighted in red and yellow, respectively. Moreover, the EPB shield is drawn with a blue line. Figure 5. Excavation steps timeline: (a) EPB penetrate soil step 1, (b) the entire EPB shield is penetrated, (c) the lining and grouting pressure is applied, (d) the EPB pressure is removed after 2, and the grout hardens, (e) EPB and lining are visible, and (f) after 50 steps of excavation. Figure 5. Excavation steps timeline: ( a ) EPB penetrate soil step 1, ( b ) the entire EPB shield is penetrated, ( c ) the lining and grouting pressure is applied, ( d ) the EPB pressure is removed after 2, and the grout hardens, (e) EPB and lining are visible, and (f) after 50 steps of excavation. 3.1. Numerical Simulation of the EPB’s Shield During mechanised excavation, the stability of the tunnel face is critical, especially in the shallow tunnels. Since the EPB shield is subject to the stress and strain of excavation, it must be simulated precisely [ 25 ]. Face pressure and EPB shape are two parameters that can be investigated for the EPB shield. The EPB face pressure has been simulated by applying a trapezoidal distribution of horizontal stress to the tunnel face; this pressure should prevent the surface of tunnels from Symmetry 2022,14, 1263 16 of 29 the burden loads in the plots of Figures 12 and 13 are chosen at different numbers. The burden load in Point A, E, and F has been assumed to be 130 kPa, while in Point B and D it is chosen at 20 kPa, and in Point C it is assumed to be 70 kPa. Furthermore, these plots are juxtaposing in the Y = 20 of the geometry. Figures 12 and 13 show different ground displacements with changes in the face pressure and the grouting pressure of the mechanised excavation. It can be observed that, at first glance, the surface settlement has an indirect relationship with face pressure and grout pressure. These plots also show the change in ground displacement from a V-shape to a W-shape after the second tunnel excavation. Figure 13 shows that the ground settlement on the surface, which is U-shaped, became wider after the second tunnel had been excavated. Symmetry 2022, 14, x FOR PEER REVIEW 16 of 29 excavation was started from Y = 40 to Y = 0. The lagging distance between the two excavation faces when two EPB machines are excavating consecutively from a different perspective was studied by [86]. The ground displacement at the crown of the tunnel is shown in Figure 12, and the surface settlement profile is presented in Figure 13. Notably, the burden loads in the plots of Figures 12 and 13 are chosen at different numbers. The burden load in Point A, E, and F has been assumed to be 130 kPa, while in Point B and D it is chosen at 20 kPa, and in Point C it is assumed to be 70 kPa. Furthermore, these plots are juxtaposing in the Y = 20 of the geometry. Figures 12 and 13 show different ground displacements with changes in the face pressure and the grouting pressure of the mechanised excavation. It can be observed that, at first glance, the surface settlement has an indirect relationship with face pressure and grout pressure. These plots also show the change in ground displacement from a V-shape to a W-shape after the second tunnel excavation. Figure 13 shows that the ground settlement on the surface, which is U-shaped, became wider after the second tunnel had been excavated. Several studies have addressed changes in the settlement value in various excavation stages: for this notion, the simulation should be studied on the longitude profile of EPB excavation. Research was carried out by [37,46], and Ref. [91] regarding the relationship between the grouting properties and settlement, or regarding the fault effect when crossing the tunnel. In their studies, they focused on the longitude profile of EPB excavation in the moving direction in each stage of excavation. Therefore, Figures 14 and 15 display the progress in the EPB and the change in the ground displacement in various stages of excavations. As shown in these figures, the ground movement was in the critical points. Each critical point was excavated with one condition in these figures to show this progress more explicitly. The burden load in Points A, E, and F has been assumed to be 130 kPa, in Points B and D it is chosen at 20 kPa, and in Point C it is assumed to be 70 kPa. Moreover, the face pressure at Points A, B, and C was assumed to be 250 kPa, and at Points D and E it was presumed to be 150 kPa, while at Point F it was taken at 50 kPa. Figure 12. Profile of the ground displacement at the crown of the tunnel at Y = 20 (Critical areas are mentioned as A–F). Figure 12. Profile of the ground displacement at the crown of the tunnel at Y = 20 (Critical areas are mentioned as A–F). Several studies have addressed changes in the settlement value in various excavation stages: for this notion, the simulation should be studied on the longitude profile of EPB excavation. Research was carried out by [ 37 , 46 ], and Ref. [ 91 ] regarding the relationship between the grouting properties and settlement, or regarding the fault effect when crossing the tunnel. In their studies, they focused on the longitude profile of EPB excavation in the moving direction in each stage of excavation. Therefore, Figures 14 and 15 display the progress in the EPB and the change in the ground displacement in various stages of excavations. As shown in these figures, the ground movement was in the critical points. Each critical point was excavated with one condition in these figures to show this progress more explicitly. The burden load in Points A, E, and F has been assumed to be 130 kPa, in Points B and D it is chosen at 20 kPa, and in Point C it is assumed to be 70 kPa. Moreover, the face pressure at Points A, B, and C was assumed to be 250 kPa, and at Points D and E it was presumed to be 150 kPa, while at Point F it was taken at 50 kPa. Symmetry 2022,14, 1263 17 of 29 Symmetry 2022, 14, x FOR PEER REVIEW 17 of 29 Figure 13. Profile of the ground displacement on the surface at Y = 20 (Critical areas are mentioned as A–F). Figure 14. Diagram of the surface settlement in the middle of the first twin tunnel based on step excavation (Critical areas are mentioned as A–F). Figure 14 shows the advancement of the EPB into the geometry from Y = 0 for the excavation of the first tunnel in six critical points. When the entire length of the EPB machine is advanced into the soil, the surface starts to settle dramatically (L_EPB_1_10 in Figure 13. Profile of the ground displacement on the surface at Y = 20 (Critical areas are mentioned as A–F). Symmetry 2022, 14, x FOR PEER REVIEW 17 of 29 Figure 13. Profile of the ground displacement on the surface at Y = 20 (Critical areas are mentioned as A–F). Figure 14. Diagram of the surface settlement in the middle of the first twin tunnel based on step excavation (Critical areas are mentioned as A–F). Figure 14 shows the advancement of the EPB into the geometry from Y = 0 for the excavation of the first tunnel in six critical points. When the entire length of the EPB machine is advanced into the soil, the surface starts to settle dramatically (L_EPB_1_10 in Figure 14. Diagram of the surface settlement in the middle of the first twin tunnel based on step excavation (Critical areas are mentioned as A–F). Symmetry 2022,14, 1263 18 of 29 Symmetry 2022, 14, x FOR PEER REVIEW 18 of 29 Figure 14). In the subsequent steps, the EPB is immediately applied to the grouting pressure and the precast lining. Therefore, the ground displacement speed decreases (L_EPB_1_10 in Figure 14). L_EPB_1_30 in Figure 14 indicates the EPB position after 30 m of excavation. Finally, line L_EPB_1_40 in Figure 14 illustrates the last step of the excavation whereby the EPB is excluded from the model, and the grouting pressure has been added to the last part of the model. It is obvious that the consolidated grouting and the implantation of the lining precast control the ground movement. Even though the EPB cutter head is highlighted in each section with the same colour, the EPB is absent in the last section of the model. Nevertheless, Figure 15 indicates the ground movement during the excavation of the second tunnel. Since these tunnels were excavated in opposite directions, their progress is reversed. Moreover, the first tunnel graphs show the excavation progress from left to right (from 0 to 40). The second tunnel has been excavated from right to left (from 40 to 0). Figures 14 and 15 indicate the surface settlement in the middle distance of the twin tunnels (x = −6), and they show the surface settlement graph when EPB has penetrated 10, 20, 30, and 40 m into the models. Additionally, Figure 16 shows the diagram of the surface settlement after 50 steps of excavation, and the ground displacement is steadied. Figures 14 and 15 indicate that the ground movement decreased subsequent to the first 30 m of the excavation, and most of the significant ground displacement occurred in the first 30 m. Notably, the first 10 m of the excavation can be attributed to the EPB length. For the next 20 metres, considerable ground displacement has occurred. This progress was similar to all critical points with different excavation properties. Figure 15. Diagram of the surface settlement in the middle of the second twin tunnel based on step excavation (Critical areas are mentioned as A–F). The ground displacement increased due to the increase in the burden load, and therefore they have a direct relationship. Figure 16 demonstrates this principle in the Y longitude after having excavated the model completely and having removed the grouting pressure. Those lines with the same colour show the surface settlement with the same burden load. After the excavation of the second tunnel, the amount of ground displacement increased. The ground displacement also increased dramatically when the burden load Figure 15. Diagram of the surface settlement in the middle of the second twin tunnel based on step excavation (Critical areas are mentioned as A–F). Figure 14 shows the advancement of the EPB into the geometry from Y = 0 for the excavation of the first tunnel in six critical points. When the entire length of the EPB machine is advanced into the soil, the surface starts to settle dramatically (L_EPB_1_10 in Figure 14). In the subsequent steps, the EPB is immediately applied to the grouting pressure and the precast lining. Therefore, the ground displacement speed decreases (L_EPB_1_10 in Figure 14). L_EPB_1_30 in Figure 14 indicates the EPB position after 30 m of excavation. Finally, line L_EPB_1_40 in Figure 14 illustrates the last step of the excavation whereby the EPB is excluded from the model, and the grouting pressure has been added to the last part of the model. It is obvious that the consolidated grouting and the implantation of the lining precast control the ground movement. Even though the EPB cutter head is highlighted in each section with the same colour, the EPB is absent in the last section of the model. Nevertheless, Figure 15 indicates the ground movement during the excavation of the second tunnel. Since these tunnels were excavated in opposite directions, their progress is reversed. Moreover, the first tunnel graphs show the excavation progress from left to right (from 0 to 40). The second tunnel has been excavated from right to left (from 40 to 0 ). Figures 14 and 15 indicate the surface settlement in the middle distance of the twin tunnels (x = − 6), and they show the surface settlement graph when EPB has penetrated 10, 20, 30, and 40 m into the models. Additionally, Figure 16 shows the diagram of the surface settlement after 50 steps of excavation, and the ground displacement is steadied. Figures 14 and 15 indicate that the ground movement decreased subsequent to the first 30 m of the excavation, and most of the significant ground displacement occurred in the first 30 m. Notably, the first 10 m of the excavation can be attributed to the EPB length. For the next 20 m , considerable ground displacement has occurred. This progress was similar to all critical points with different excavation properties. Symmetry 2022,14, 1263 19 of 29 Symmetry 2022, 14, x FOR PEER REVIEW 19 of 29 increased. The average ground displacement at Point E, at longitude direction was 5.8 cm when the burden load was 20 kPa. The burden load was then increased to 70 kPa, and the average settlement reached 7.84 cm. Finally, when the burden load was 130 KPa, the average ground displacement was 10.45 cm. These results show the burden load effect when the mechanised properties are the same, and the tunnel is excavated with constant face pressure (250 KPa) and grouting pressure (300 kPa). As already mentioned, Points E and F have the greatest ground displacement. During the excavation of those areas, various geotechnical issues arose, and the process of excavation was halted [92]. Moreover, Figure 16 clearly shows the changing surface during the excavation of the first and second tunnels with the same excavation properties. It should be borne in mind that Figure 16 shows line in X = −6, whereas Figure 2 shows this point location on the surface. Figure 16. Longitudinal profile with the same face pressure and changing the burden load after excavation of both tunnels (Critical areas are mentioned as A–F). 4.4. Statistical Analysis on Results of Extra Models Nine simulations were developed based on the real situation, and 45 models were developed for additional studies with changes made in the excavation properties and boundary conditions, thereby producing 54 models. By relaying these extra models, the surface settlement has been investigated for a wide variety of excavation properties and geometry. A statistical study has been carried out for those surface settlements mentioned in the previous section. This statistical study reveals the surface settlement relation with excavation properties and geometry. As shown in Figures 11 to 16, as the excavation properties and burden load in all critical points changed, ground displacement also changed. Furthermore, the critical points in this paper present several major differences in their soil layers, water table, and geometry. In one case, Point F of the Seville metro was excavated into the sand layer at 10 m, while Point B was excavated into marl at 14 m. Moreover, the water table has different levels in these critical points. The water table in Point D was located 1 m below the surface; in other points, this water table was approximately 4–5 m from the surface. Figure 17 shows the various surface settlements that include different face and burden pressures. Figure 17 sorts the height of the surface settlement from the Figure 16. Longitudinal profile with the same face pressure and changing the burden load after excavation of both tunnels (Critical areas are mentioned as A–F). The ground displacement increased due to the increase in the burden load, and therefore they have a direct relationship. Figure 16 demonstrates this principle in the Y longitude after having excavated the model completely and having removed the grouting pressure. Those lines with the same colour show the surface settlement with the same burden load. After the excavation of the second tunnel, the amount of ground displacement increased. The ground displacement also increased dramatically when the burden load increased. The average ground displacement at Point E, at longitude direction was 5.8 cm when the burden load was 20 kPa. The burden load was then increased to 70 kPa, and the average settlement reached 7.84 cm. Finally, when the burden load was 130 KPa, the average ground displacement was 10.45 cm. These results show the burden load effect when the mechanised properties are the same, and the tunnel is excavated with constant face pressure (250 KPa) and grouting pressure (300 kPa). As already mentioned, Points E and F have the greatest ground displacement. During the excavation of those areas, various geotechnical issues arose, and the process of excavation was halted [ 92 ]. Moreover, Figure 16 clearly shows the changing surface during the excavation of the first and second tunnels with the same excavation properties. It should be borne in mind that Figure 16 shows line in X = −6, whereas Figure 2shows this point location on the surface. 4.4. Statistical Analysis on Results of Extra Models Nine simulations were developed based on the real situation, and 45 models were developed for additional studies with changes made in the excavation properties and boundary conditions, thereby producing 54 models. By relaying these extra models, the surface settlement has been investigated for a wide variety of excavation properties and geometry. A statistical study has been carried out for those surface settlements mentioned in the previous section. This statistical study reveals the surface settlement relation with excavation properties and geometry. As shown in Figures 11–16, as the excavation properties and burden load in all critical points changed, ground displacement also changed. Furthermore, the critical points in this paper present several major differences in their soil layers, water table, and geometry. In one case, Point F of the Seville metro was excavated Symmetry 2022,14, 1263 20 of 29 into the sand layer at 10 m, while Point B was excavated into marl at 14 m. Moreover, the water table has different levels in these critical points. The water table in Point D was located 1 m below the surface; in other points, this water table was approximately 4–5 m from the surface. Figure 17 shows the various surface settlements that include different face and burden pressures. Figure 17 sorts the height of the surface settlement from the highest value down to lowest: the greatest surface settlements are highlighted in red and low surface settlements are pink. Symmetry 2022, 14, x FOR PEER REVIEW 20 of 29 highest value down to lowest: the greatest surface settlements are highlighted in red and low surface settlements are pink. Figure 17. Surface settlement for 54 simulations with different EPB face pressure (first suffix number) and burden load (second suffix number). Table 6 shows the average surface settlement across all the critical points, where they have the same burden load and face pressure. In this part, the difference in geometry has been waived, and the study has focused solely on face pressure and burden load, whereby all simulations that excavated with the same face pressure and burden load are located in the same group and the average of surface settlement is reported. Table 6 illustrates the lowest average of surface settlement when the face pressure and burden loads are 250 kPa and 20 kPa, respectively. On the other hand, the highest average of surface settlement is for 50 kPa face pressure and 130 kPa burden load. The highest and lowest average surface settlement is 0.121866 m and 0.049275 m, respectively. The percentage difference between the highest and lowest average of surface settlement is 84.83%. In the simulations where burden load is assumed at 20 kPa, the average surface settlement increased by 50.9% when the face pressure decreases from 250 kPa to 50 kPa. Moreover, when the face pressure is assumed to be 250 kPa, the burden load increases from 20 to 130 kPa and the average surface settlement reaches 68.4%. Table 6. Average surface settlement. Burden Load (kPa) Face Pressure (kPa) Average Surface Settlement (m) 20 50 0.074351 150 0.056472 250 0.049275 70 50 0.095415 150 0.072697 250 0.064137 130 50 0.121866 150 0.093389 250 0.082971 Figure 17. Surface settlement for 54 simulations with different EPB face pressure (first suffix number) and burden load (second suffix number). Table 6shows the average surface settlement across all the critical points, where they have the same burden load and face pressure. In this part, the difference in geometry has been waived, and the study has focused solely on face pressure and burden load, whereby all simulations that excavated with the same face pressure and burden load are located in the same group and the average of surface settlement is reported. Table 6illustrates the lowest average of surface settlement when the face pressure and burden loads are 250 kPa and 20 kPa, respectively. On the other hand, the highest average of surface settlement is for 50 kPa face pressure and 130 kPa burden load. The highest and lowest average surface settlement is 0.121866 m and 0.049275 m, respectively. The percentage difference between the highest and lowest average of surface settlement is 84.83%. In the simulations where burden load is assumed at 20 kPa, the average surface settlement increased by 50.9% when the face pressure decreases from 250 kPa to 50 kPa . Moreover, when the face pressure is assumed to be 250 kPa, the burden load increases from 20 to 130 kPa and the average surface settlement reaches 68.4%. Figure 18 shows the percentage changes between average surface settlement with changing the face pressure and burden load. The percentage changes are evaluated by V2−V1 V1 × 100, where V denotes different numbers. The negative and positives values in this table indicate an increment or decrement of surface settlement. It has been observed that the average surface settlement increases to 147.32% if face pressure changes from 250 kPa to 50 kPa and the burden load increases from 20 kPa to 130 kPa. In contrast, the lowest Symmetry 2022,14, 1263 21 of 29 percentage changes occur in various models where either the face pressure or the burden load remains constant. Furthermore, if the face pressure and burden load change from 150 kPa and 20 kPa to 50 kPa and 130 kPa, the surface settlement increases to 115.8%. Table 6. Average surface settlement. Burden Load (kPa) Face Pressure (kPa) Average Surface Settlement (m) 20 50 0.074351 150 0.056472 250 0.049275 70 50 0.095415 150 0.072697 250 0.064137 130 50 0.121866 150 0.093389 250 0.082971 Symmetry 2022, 14, x FOR PEER REVIEW 21 of 29 Figure 18 shows the percentage changes between average surface settlement with changing the face pressure and burden load. The percentage changes are evaluated by V2−V1 V1 × 100, where V denotes different numbers. The negative and positives values in this table indicate an increment or decrement of surface settlement. It has been observed that the average surface settlement increases to 147.32% if face pressure changes from 250 kPa to 50 kPa and the burden load increases from 20 kPa to 130 kPa. In contrast, the lowest percentage changes occur in various models where either the face pressure or the burden load remains constant. Furthermore, if the face pressure and burden load change from 150 kPa and 20 kPa to 50 kPa and 130 kPa, the surface settlement increases to 115.8%. Figure 18. Percentage change from average surface settlement. Figures 19 and 20 show the separate effects of the points and the face pressure and burden load effects on surface settlement. Figures 19 and 20 reveal that the surface settlement in Points A, B, and C differs from that in Points D, E, and F. Surface settlement in Points A, B, and C is generally lower than in Points D, E, and F. The surface settlement in Points A, B, and C changed from 0.037 to 0.089 m. However, surface settlement for Points D, E, and F changed from 0.059 to 0.169 m. Furthermore, varying the span is estimated to be 0.052 m for Points A, B, and C and 0.109 m for Points C, D, and E. This difference is related to the geometry of these points. Points A, B, and C are located in gravel and marl, while Point D is located in three different soil layers (clay, sand, and gravel), and Points E and F are located in the sand layer. Figures 19 and 20 support the previous assumption that the surface settlement has a direct relationship with the burden load. These assumptions are shown in all the points. However, the changing surface settlement in Points C, D, and E is considerably more Figure 18. Percentage change from average surface settlement. Figures 19 and 20 show the separate effects of the points and the face pressure and burden load effects on surface settlement. Figures 19 and 20 reveal that the surface settlement in Points A, B, and C differs from that in Points D, E, and F. Surface settlement in Points A, B, and C is generally lower than in Points D, E, and F. The surface settlement in Symmetry 2022,14, 1263 22 of 29 Points A, B, and C changed from 0.037 to 0.089 m. However, surface settlement for Points D, E, and F changed from 0.059 to 0.169 m. Furthermore, varying the span is estimated to be 0.052 m for Points A, B, and C and 0.109 m for Points C, D, and E. This difference is related to the geometry of these points. Points A, B, and C are located in gravel and marl, while Point D is located in three different soil layers (clay, sand, and gravel), and Points E and F are located in the sand layer. Symmetry 2022, 14, x FOR PEER REVIEW 22 of 29 scattered. Moreover, surface settlement in Points D, E, and F is much more scattered when either the face pressure is 50 kPa or the burden load is 130 kPa, which comes from the soil layer effect. Figure 19. The relationship of burden load with surface settlement. Figure 20. The relationship of face pressure with surface settlement. Below, the models have been divided into four categories based on the excavated soil. The X axial in Figure 21 is shown, with an increase in Young’s modulus the surface settlement decreased after excavation of both tunnels in comparison with lower Young’s modulus. Moreover, excavation into the sand and clay with lower Young’s modulus leads simulation to more ground displacement. Figure 21 shows the tunnel depth effect on the simulation. As shown, most of the models were developed within the 14 m to 12 m range and only 18 models were developed at the depth of 10 m. However, Figure 21 indicates that the distribution of the tunnel excavated 10 m deep was greater than in other points. Furthermore, Figure 21 illustrates the effect of Young’s modulus for the tunnel excavated at different depths. Tunnels were placed into the 13 m depth and excavated in two Figure 19. The relationship of burden load with surface settlement. Symmetry 2022, 14, x FOR PEER REVIEW 22 of 29 scattered. Moreover, surface settlement in Points D, E, and F is much more scattered when either the face pressure is 50 kPa or the burden load is 130 kPa, which comes from the soil layer effect. Figure 19. The relationship of burden load with surface settlement. Figure 20. The relationship of face pressure with surface settlement. Below, the models have been divided into four categories based on the excavated soil. The X axial in Figure 21 is shown, with an increase in Young’s modulus the surface settlement decreased after excavation of both tunnels in comparison with lower Young’s modulus. Moreover, excavation into the sand and clay with lower Young’s modulus leads simulation to more ground displacement. Figure 21 shows the tunnel depth effect on the simulation. As shown, most of the models were developed within the 14 m to 12 m range and only 18 models were developed at the depth of 10 m. However, Figure 21 indicates that the distribution of the tunnel excavated 10 m deep was greater than in other points. Furthermore, Figure 21 illustrates the effect of Young’s modulus for the tunnel excavated at different depths. Tunnels were placed into the 13 m depth and excavated in two Figure 20. The relationship of face pressure with surface settlement. Figures 19 and 20 support the previous assumption that the surface settlement has a direct relationship with the burden load. These assumptions are shown in all the points. However, the changing surface settlement in Points C, D, and E is considerably more scattered. Moreover, surface settlement in Points D, E, and F is much more scattered when either the face pressure is 50 kPa or the burden load is 130 kPa, which comes from the soil layer effect. Below, the models have been divided into four categories based on the excavated soil. The X axial in Figure 21 is shown, with an increase in Young’s modulus the surface Symmetry 2022,14, 1263 23 of 29 settlement decreased after excavation of both tunnels in comparison with lower Young’s modulus. Moreover, excavation into the sand and clay with lower Young’s modulus leads simulation to more ground displacement. Figure 21 shows the tunnel depth effect on the simulation. As shown, most of the models were developed within the 14 m to 12 m range and only 18 models were developed at the depth of 10 m. However, Figure 21 indicates that the distribution of the tunnel excavated 10 m deep was greater than in other points. Furthermore, Figure 21 illustrates the effect of Young’s modulus for the tunnel excavated at different depths. Tunnels were placed into the 13 m depth and excavated in two different Young’s moduli; tunnels excavated in the lower Young’s modulus experienced greater surface settlement. Notably, tunnels excavated into sand are located 10 m deep. These simulations had more scatter results to estimate surface settlement than other simulations. In Figure 22, points with the greatest surface settlement were Points E and F. Additionally, the burden load was 130 kPa, and EPB face pressure was 50 kPa. Symmetry 2022, 14, x FOR PEER REVIEW 23 of 29 different Young’s moduli; tunnels excavated in the lower Young’s modulus experienced greater surface settlement. Notably, tunnels excavated into sand are located 10 m deep. These simulations had more scatter results to estimate surface settlement than other simulations. In Figure 22, points with the greatest surface settlement were Points E and F. Additionally, the burden load was 130 kPa, and EPB face pressure was 50 kPa. Figure 21. Boxplots showing surface settlements for different values of Young’s Modulus: Changing soil properties and tunnel depth. Surface settlement for Points D, E, and F has more dispersion than Points A, B, and C. The variance for Points A, B, C, D, E, and F is 2.42 × 10−4, 2.06 × 10−4, 2.23 × 10−4, 6.15 × 10−4, 1.01 × 10−3 and 1.10 × 10−3. Figure 22 shows the dispersion of surface settlement in each point. Figure 22. Boxplots showing surface settlements for different locations. The water table models are divided into three categories containing nine models, whereby (Point D) water table is placed at a depth of 1 m, and the other models of the water table are placed at approximately 4 and 5 m in depth. Figure 23 presents these criteria. On the other hand, models where the water table was placed at a depth of 1 m are Figure 21. Boxplots showing surface settlements for different values of Young’s Modulus: Changing soil properties and tunnel depth. Symmetry 2022, 14, x FOR PEER REVIEW 23 of 29 different Young’s moduli; tunnels excavated in the lower Young’s modulus experienced greater surface settlement. Notably, tunnels excavated into sand are located 10 m deep. These simulations had more scatter results to estimate surface settlement than other simulations. In Figure 22, points with the greatest surface settlement were Points E and F. Additionally, the burden load was 130 kPa, and EPB face pressure was 50 kPa. Figure 21. Boxplots showing surface settlements for different values of Young’s Modulus: Changing soil properties and tunnel depth. Surface settlement for Points D, E, and F has more dispersion than Points A, B, and C. The variance for Points A, B, C, D, E, and F is 2.42 × 10−4, 2.06 × 10−4, 2.23 × 10−4, 6.15 × 10−4, 1.01 × 10−3 and 1.10 × 10−3. Figure 22 shows the dispersion of surface settlement in each point. Figure 22. Boxplots showing surface settlements for different locations. The water table models are divided into three categories containing nine models, whereby (Point D) water table is placed at a depth of 1 m, and the other models of the water table are placed at approximately 4 and 5 m in depth. Figure 23 presents these criteria. On the other hand, models where the water table was placed at a depth of 1 m are Figure 22. Boxplots showing surface settlements for different locations. Surface settlement for Points D, E, and F has more dispersion than Points A, B, and C. The variance for Points A, B, C, D, E, and F is 2.42 × 10 −4 , 2.06 × 10 −4 , 2.23 × 10 −4 , Symmetry 2022,14, 1263 24 of 29 6.15 ×10−4 , 1.01 × 10 −3 and 1.10 × 10 −3 . Figure 22 shows the dispersion of surface settlement in each point. The water table models are divided into three categories containing nine models, whereby (Point D) water table is placed at a depth of 1 m, and the other models of the water table are placed at approximately 4 and 5 m in depth. Figure 23 presents these criteria. On the other hand, models where the water table was placed at a depth of 1 m are found to suffer from more ground displacement. However, outliers reported in Figure 23 show more surface settlement in Points E and F. Therefore, this research presents certain limitations regarding generalisations on the water table and surface settlement. Symmetry 2022, 14, x FOR PEER REVIEW 24 of 29 found to suffer from more ground displacement. However, outliers reported in Figure 23 show more surface settlement in Points E and F. Therefore, this research presents certain limitations regarding generalisations on the water table and surface settlement. Figure 23. Boxplot showing surface settlement for three different depths of water table. 5. Conclusions During the excavation of the tunnels for the metro of Seville, engineers faced many issues, with surface displacement setting a tremendous challenge. This study strives to highlight its location and the main reason the surface displacement increased during excavation. • The ground displacement in crowded locations should be controlled and measured step by step, particularly where twin tunnels in the crowded location are excavated into sand. For instance, in this study, Points E and F were placed in sand. The average surface settlement for Points A, B, and C was −0.058, −0.057, and −0.06 m, respectively. However, the surface settlement average in Points D, E, and F were −0.093, −0.101, and −0.105 m, respectively. Moreover, the average surface settlement in Point F is 1.8 times more than in Point A, and the average surface settlement in Point D is 1.5 times greater than that in Point C. • Face pressure and grout pressure have an indirect relationship with the surface settlement, and this study shows with increasing these pressures the surface settlement is decreased; it is notable that this pressure must be measured in a decent range because increasing the EPB face pressure and grouting pressure can cause swelling on the surface or hydraulic fracture in the soil. • Excavating tunnels at greater depths or in stronger soil types (high Young’s modulus) would decrease the surface settlement in the tunnels of the metro of Seville. It can be observed that, in Points F and E, the shallower and weaker the soil is at the excavation site, the higher the variance. The range of surface settlement was changing approximately between 0.04 and 0.08 m, where Young’s modulus was 83.3 and 58.8. On the other hand, surface settlement was changing between 0.06 and 0.16 m when the Young’s modulus was set at 24.5. However, changing the depth of the tunnel could impact the results of this study. It is recommended that for tunnel excavation in the soft soil of Seville in the future, with increasing the depth of tunnel and grout pressure or assuming the ground treatment, the surface settlement would be reduced. • The water table and pore pressure of the soil also commands high priority. The existence of a shallow water table in gravel leads to a high amount of ground displacement into the crown of the tunnel and increases the surface settlement. For instance, Figure 23. Boxplot showing surface settlement for three different depths of water table. 5. Conclusions During the excavation of the tunnels for the metro of Seville, engineers faced many issues, with surface displacement setting a tremendous challenge. This study strives to highlight its location and the main reason the surface displacement increased during excavation. • The ground displacement in crowded locations should be controlled and measured step by step, particularly where twin tunnels in the crowded location are excavated into sand. For instance, in this study, Points E and F were placed in sand. The average surface settlement for Points A, B, and C was − 0.058, − 0.057, and −0.06 m, respectively. However, the surface settlement average in Points D, E, and F were − 0.093, − 0.101, and − 0.105 m, respectively. Moreover, the average surface settlement in Point F is 1.8 times more than in Point A, and the average surface settlement in Point D is 1.5 times greater than that in Point C. • Face pressure and grout pressure have an indirect relationship with the surface settlement, and this study shows with increasing these pressures the surface settlement is decreased; it is notable that this pressure must be measured in a decent range because increasing the EPB face pressure and grouting pressure can cause swelling on the surface or hydraulic fracture in the soil. • Excavating tunnels at greater depths or in stronger soil types (high Young’s modulus) would decrease the surface settlement in the tunnels of the metro of Seville. It can be observed that, in Points F and E, the shallower and weaker the soil is at the excavation site, the higher the variance. The range of surface settlement was changing approximately between 0.04 and 0.08 m, where Young’s modulus was 83.3 and 58.8. On the other hand, surface settlement was changing between 0.06 and 0.16 m when the Young’s modulus was set at 24.5. However, changing the depth of the tunnel could impact the results of this study. It is recommended that for tunnel excavation in the soft soil of Seville in the future, with increasing the depth of tunnel and grout pressure or assuming the ground treatment, the surface settlement would be reduced. Symmetry 2022,14, 1263 25 of 29 • The water table and pore pressure of the soil also commands high priority. The existence of a shallow water table in gravel leads to a high amount of ground displacement into the crown of the tunnel and increases the surface settlement. For instance, the excavation process in Seville was halted for months in Point D for this reason. Moreover, ground treatment is suggested for dealing with the location where the water table is high. It is also recommended to decelerate EPB until the grouting has enough time to harden. • In tunnels excavated at a lower depth, surface settlement would be greater than in tunnels excavated at a deeper depth. The result was clear tunnels located on 10 m of surface that have more surface settlement than other tunnels. However, the soil properties have to be considered too. • In these simulations, when tunnels are excavated in the E and F, the burden load, face pressure, and grout pressure were 130, 50, and 100 kPa, respectively. The surface displacement at both points was 0.16 m. This amount of surface displacement is too high. Therefore, it is recommended when tunnels want to be excavated in these areas and the burden load is 130 kPa. The face pressure and grout pressure would be increased to 250 kPa and 300 kPa, respectively. Then the numerical study in FLAC3D shows that the surface displacement is reduced to approximately 0.10 m. Therefore, this sample shows the importance of face pressure and grout pressure in the underground excavation with EPB and controlling the surface displacement. Finally, an increase in EPB pressure is also recommended for future excavation for other metro lines in crowded areas with a high burden load. It is clear that all parameters of excavation and soil properties should be considered next to each other because these parameters have an impact on each other. Author Contributions: Conceptualisation, M.B., E.J.M.-A., R.R.-H. and A.J.-M.; Data curation, M.B.; Funding acquisition, M.B. and E.J.M.-A.; Investigation, M.B. and E.J.M.-A.; Methodology, E.J.M.-A. and R.R.-H.; Resources, R.R.-H. and C.S.-C.; Software, M.B.; Supervision, E.J.M.-A., R.R.-H., C.S.-C. and A.J.-M.; Validation, E.J.M.-A., R.R.-H., M.K., C.S.-C. and A.J.-M.; Visualization, M.K.; Writing—original draft, M.B.; Writing—review & editing, E.J.M.-A., M.K. and C.S.-C. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Not applicable. Acknowledgments: This research has been supported by Itasca Consulting Group through its Itasca Educational Partnership (IEP) Programme, under whose framework the FLAC3D software licence used for the research has been provided. Therefore, the authors would like to thank the ITASCA consulting group for their support and guidance. Additionally, the authors would like to thank Pedro Arozamena-Cagigal for his assistance and support in this research. Conflicts of Interest: The authors declare no conflict of interest. References 1. Department of Economic & Social Affairs of United Nations. World Urbanization Prospects: The 2018 Revision; United Nations: New York, NY, USA, 2019; Volume 12. 2. Williams, M.; Edgeworth, M.; Zalasiewicz, J.; Waters, C.N.; Steffen, W.; Wolfe, A.P.; Minter, N.J.; Cearreta, A.; Gałuszka, A. ; Haff, P.; et al. Underground metro systems. A durable geological proxy of rapid urban population growth and energy consumption during the Anthropocene. In The Routledge Companion to Big History; Routledge, Taylor & Francis Group: London, UK, 2019; pp. 434–455. 3. Qiao, Y.-K.; Peng, F.-L.; Sabri, S.; Rajabifard, A. Socio-environmental costs of underground space use for urban sustainability. Sustain. Cities Soc. 2019,51, 101757. [CrossRef] 4. Li, X.; Li, C.; Parriaux, A.; Wu, W.; Li, H.; Sun, L.; Liu, C. Multiple resources and their sustainable development in Urban Underground Space. Tunn. Undergr. Space Technol. 2016,55, 59–66. [CrossRef]