A geometric model is established, assuming the tunnel is excavated along the positive y-axis, which y = 0 represents the location of the tunnel’s face. The x-axis lies on the ground surface and is perpendicular to the tunnel axis, while the z-axis represents the vertical direction, as illustrated in Figure 1.

Neglecting the elliptical deformation of segment lining resulting from long-term service, and based on the cavity contraction theory within soil proposed by Sen (1951), Mindlin and Cheng (1950), as well as the theory of uniform radial contraction proposed by Sagaseta (1987), it is assumed that the volume loss caused by lining contraction during excavation is V = 2 π R u ε . Pinto and Whittle (2014) derived theoretical solutions for the three-dimensional surface deformation field, which are expressed as Equations 1–3: u x 0 = − V π · 1 − v x x 2 + H 2 · x 2 + y 2 + H 2 − y x 2 + y 2 + H 2 (1) u y 0 = − V π · 1 − v x 2 + z 2 + H 2 (2) u z 0 = − V π · 1 − v H x 2 + H 2 · x 2 + y 2 + H 2 − y x 2 + y 2 + H 2 (3)

In these equations, u x 0 , u y 0 and u z 0 represent the displacement components in the x, y, and z directions, respectively; ν denotes the Poisson’s ratio of the soil or rock mass; u ε is the radial contraction of the tunnel lining; R and H are the tunnel radius and the depth of the tunnel axis, respectively.

For Equation 3, by neglecting soil compressibility and taking Poisson’s ratio of the ground as 0.5, and setting y = − ∞ to represent a stabilized settlement condition, the vertical displacement along the x-direction can be simplified as Equation 4: u z 0 y → − ∞ v = 0.5 = − 2 u ε R H x 2 + H 2 (4)

In fact, Equation 4 is identical to the two-dimensional analytical solution proposed by Verruijt and Booker (1996) for surface deformation induced by tunnel excavation. However, field monitoring data from actual engineering projects have shown that this model does not provide accurate predictions of surface settlement. Therefore, Yang et al. (2020) introduced an improved formulation for the displacement component u z 0 by incorporating correction coefficients in three spatial dimensions. As a result, Equation 3 was modified to Equation 5, and Equation 4 was refined into Equation 6: u z 0 ′ = − λ 1 u ε R · H λ 2 x 2 + H 2 · λ 2 x 2 + λ 3 y 2 + H 2 − λ 3 y λ 2 x 2 + λ 3 y 2 + H 2 (5) u z 0 ′ y → − ∞ v = 0.5 = − 2 λ 1 u ε R H λ 2 x 2 + H 2 (6)

In these equations, λ 1 is the correction parameter for the maximum settlement value, λ 2 is the correction parameter for the settlement trough width, and λ 3 is the correction parameter for the settlement trough length.

The correction principle in the x-axis direction is illustrated in Figure 2. Clearly, the value of V remains unchanged before and after the correction. Therefore: ∫ − ∞ + ∞ ∫ − ∞ + ∞ u z 0 ′ d x d y = ∫ − ∞ + ∞ ∫ − ∞ + ∞ u z 0 d x d y (7)

By solving Equation 7 with respect to the x-axis and the x-y plane, respectively, the results are presented in Equations 8, 9: λ 1 = λ 2 = λ (8) λ 3 = 1 (9)

Thus, the three-dimensional surface deformation induced by tunnel excavation can be corrected and expressed as Equation 10: u z 0 ′ = − λ u ε R · H λ x 2 + H 2 · λ x 2 + y 2 + H 2 − y λ x 2 + y 2 + H 2 (10)

Therefore, the vertical deformation in the x-direction can be simplified as Equation 11: u z 0 ′ y → − ∞ v = 0.5 = − 2 λ ε R 2 H λ x 2 + H 2 (11)

The settlement expression in the y-direction is given by (see Equation 12): u z 0 ′ x = 0 v = 0.5 = − λ u ε R H · y 2 + H 2 − y y 2 + H 2 (12)

In their study, Yang et al. (2020) determined the empirical value of λ through model tests, as shown in Equation 13: λ = 0.514 + 3.356 e − 2.466 R H (13)

The correction principle in the y-direction is illustrated in Figure 3. The correction only amplifies the maximum surface deformation by a factor of λ , without altering the distribution pattern or variation trend of the surface deformation. This implies that the second derivative of the surface settlement distribution curve in the y-direction remains the same before and after the correction. Moreover, in both solutions, the surface settlement directly above the tunnel face corresponds to half of the maximum post-excavation surface settlement.

However, both field monitoring and subsequent numerical simulation studies indicate that in mudstone strata, the assumption of y = 0 representing the position of the tunnel face may not be appropriate due to the significant support provided by the shield shell. Therefore, in addition to the conventional approach for calculating surface deformation induced by shallow-buried shield tunneling in mudstone, this study also adopts an alternative method in which y = 0 denotes the location of tail grouting.

This is not a simple translation of the settlement surface along the negative y-axis by a distance equal to the length of the shield shell l ; rather, it requires the introduction of a parameter y ′ as the coordinate reference along the y-axis. The difference between the two approaches lies in the value of y ′ : the former method uses y ′ = y , while the latter uses y ′ = y + l . For the case study project—Chongqing Metro Line 27—the shield shell length measured on-site is l = 10.8 m , as shown in Figure 3 and Equation 14. u z 0 ′ = − λ u ε R · H λ x 2 + H 2 · λ x 2 + y ′ 2 + H 2 − y ′ λ x 2 + y ′ 2 + H 2 = − λ u ε R · H λ x 2 + H 2 · λ x 2 + y 2 + H 2 − y λ x 2 + y 2 + H 2 or − λ u ε R · H λ x 2 + H 2 · λ x 2 + y + l 2 + H 2 − y + l λ x 2 + y + l 2 + H 2 (14)

The shield machine concurrently performs both tunnel face excavation and support of the surrounding ground. As the segment lining is assembled, the gap at the shield tail is immediately filled with grout. The grout gradually hardens, and the grouting pressure slowly dissipates, thereby providing support to the surrounding strata. The factors involved in this process are illustrated in Figure 4.

In constructing the model, the geological conditions at the site of Chongqing Metro Line 27 were considered. The detailed geological profile is shown in Figure 6. The ground strata were defined as mudstone, overlain by a 5-m layer of plain fill. The established model consists of a total of 393135 finite elements. Material properties used in the simulation are listed in Table 1. Based on a typical section of Chongqing Metro Line 27, a finite element model of shield tunneling was developed using COMSOL Multiphysics, as illustrated in Figure 7. In the model, the positive direction of the y-axis corresponds to the tunnel excavation direction.

Based on the assumptions and analysis described above and considering that the grout hardening process spans ten segment rings as well as the need to minimize boundary effects, the model length along the tunnel axis was set to the equivalent of 50 segment rings. As a result, the numerical simulation domain was defined with dimensions of 45 m × 90 m × 45 m (x × y × z), respectively. The tunnel excavation diameter is D = 8.83 m, making the model length along the tunnel axis approximately 10.23D, and the width and height about 5.11D. These dimensions satisfy the commonly accepted range of 3–5 times the excavation diameter for minimizing boundary influence in tunnel simulations, allowing boundary effects to be reasonably neglected. To facilitate gravity loading, the vertical coordinate of the ground surface (Z) was set to 0 m. Similarly, the tunnel axis was aligned such that the X-coordinate at the centerline is 0 m, and the excavation starting point was assigned a Y-coordinate of 0 m.

The shield advancement rate was assumed to be six rings per day under ideal conditions. At each construction step, all associated processes—such as excavation, segment installation, grouting, and grout hardening—advance by exactly one ring. Therefore, an increment of one in the step parameter T corresponds to 4 h of actual construction time. Considering the time requirements of various construction processes, a total of 68 steps were required to complete the full construction sequence in the numerical model, which covers a tunnel length of 50 segment rings along the tunnel axis. The key construction stages and their corresponding step parameters are summarized in Table 2.

The development of vertical surface displacement at the middle cross-section perpendicular to the tunnel axis (Y = 45 m) is shown in Figure 8. It can be observed that the surface settlement in all models approximately follows a normal distribution, which is consistent with the surface loss theory proposed by Peck. Figure 8a illustrates the vertical surface displacement when the hardening process of the grout is not considered. When the cutterhead reaches the location, the vertical surface displacement is −0.26 mm, and it reaches −2.02 mm as the shield tail passes. Ultimately, the surface experiences a maximum settlement (most negative displacement) of −4.36 mm. The distribution of surface displacement for grout initial setting times of 1 h, 2 h, and 4 h is shown in Figures 8b–d. The maximum settlement (most negative displacement) reaches −5.60 mm, −6.51 mm, and −7.49 mm, which represent increases of 19%, 39%, and 59%, respectively, compared to the case with an initial setting time of 0 h.

As shown in Figure 9, the final vertical displacements along the cross-sections at positions x = 0m, 5m, 10m, 15m, and 20 m on the middle cross-section (Y = 45 m) of the model are plotted. This is done to analyze the impact of different grout initial setting times on the vertical displacement of the surrounding strata. It can be observed that the further the distance from the tunnel center, the more stable the strata become.

When the initial setting time is 1 h, as shown in Figure 9a for the position x = 0 m: The strata above the tunnel exhibit settlement, and the settlement value fluctuates and increases with depth. Below the tunnel, a gradually decreasing heave is observed, with the maximum settlement and heave values at the crown and invert of the lining being −9.62 mm and 25.13 mm, respectively.

At position x = 5m, as shown in Figure 9b, with increasing depth, the vertical displacement of the strata initially experiences settlement followed by heave. In the shallow layer (z > −12 m), the settlement value remains relatively stable at approximately −4.9 mm. Subsequently, the settlement decreases rapidly and transforms into heave, reaching a maximum heave value of 9.29 mm near z = −25 m.

At position x = 10m, as shown in Figure 9c, the vertical displacement of the strata with increasing depth shows the largest settlement at the surface, with a maximum value of −2.96 mm. At z = −21m, the displacement reaches zero and then turns into heave, with the maximum heave value of 3.75 mm near z = −31 m.

At position x = 15m, as shown in Figure 9d, the vertical displacement of the strata initially shows a maximum settlement value of −1.27 mm at the surface, which decreases until it reaches an inflection point at z = −15 m. The settlement then increases again, reaching −0.79 mm at z = −22 m. As the depth increases, the settlement value gradually decreases, and at z < −30m, the vertical displacement of the strata turns into heave, with the maximum heave of 1.1 mm at z = −36.5 m.

At position x = 20m, as shown in Figure 9e, with a grout initial setting time of 1 h, the vertical displacement of the strata first shows heave followed by settlement with increasing depth. The maximum heave of 0.22 mm is observed at z = -9m, after which it decreases. At z = −15m, the heave becomes zero, and at z = −27m, the settlement reaches its inflection point with a maximum settlement value of −0.88 mm.

The development of horizontal strata displacement at a depth of Z = −19.4m, which corresponds to the tunnel axis depth, is shown in Figure 10. The horizontal displacement follows an antisymmetric distribution about the tunnel axis.

Figure 10a illustrates the horizontal displacement development when the grout initial setting time is 0 h. Before the shield tail passes, the horizontal displacement is within 2 mm. As the shield tail passes and the grouting is performed, the horizontal displacement increases towards the tunnel arch waist along both sides of the model, reaching a maximum value of 5.5 mm. Subsequently, the horizontal displacement rebounds and decreases. After 4 h of grouting, the horizontal displacement at the tunnel arch waist decreases to 2.4 mm, and when the grouting pressure dissipates (20 h after grouting), it stabilizes at 2.1 mm. After 40 h of grouting, the horizontal displacement at the tunnel arch waist is 2.3 mm, which is a 58% reduction compared to the horizontal displacement at the time when the shield tail passed and grouting began.

The development of horizontal displacement for grout initial setting times of 1 h, 2 h, and 4 h is shown in Figures 10b–d. During the grout hardening process, especially at the beginning of grouting, the slower the initial setting of the grout, the larger the horizontal displacement at the same moment. After 4 h of grouting, the horizontal displacements at the tunnel arch waist corresponding to the grout initial setting times of 1 h, 2 h, and 4 h are 2.4 mm, 2.9 mm, 3.1 mm, and 3.7 mm, respectively.

The initial setting time used in the field is approximately 1 h. Therefore, the calculation results for a grout initial setting time of 1 h are compared with the data obtained from on-site monitoring. The vertical surface displacement variation at the center of the surface in the model is extracted and compared with the measured vertical surface displacement at the monitoring point above the tunnel axis in the project, as shown in Figure 11.

Clearly, the calculation results are in close agreement with the on-site monitoring data. This indicates that the established numerical model can reasonably reflect the actual strata response observed in the field. Furthermore, the model can converge the errors within a relative tolerance through iterative calculations, supporting the validity of the model.

Based on the variation patterns of vertical surface and strata displacement presented earlier, it can be concluded that vertical displacement magnitude increases as the grout initial setting time increases. Additionally, grouts with slower setting times, due to their delayed effect, reduce the supporting capacity of the strata, weakening the stress coordination ability of the strata. This leads to an increase in the vertical displacement difference between the strata on both sides of the tunnel. The large settlement zone caused by tunnel construction stabilizes within a range of 15 m on either side of the tunnel axis, with more significant increases in settlement values.

Furthermore, for the strata farther from the tunnel axis, the influence of a slower-setting grout layer results in more pronounced surface heave. This suggests a more apparent trend of overall shear deformation of the geotechnical body.

As the grout initial setting time increases, the maximum value of vertical settlement also increases. As seen in Figure 9c, the depth at which the displacement value is zero gradually rises. This indicates that the increased grout initial setting time causes the strata at this depth to transition from settlement to heave. Figure 9e also shows that at certain depths; vertical displacement shifts from settlement to heave due to the delayed hardening process.

At the tunnel axis depth, horizontal displacement shows that the grout initial setting time has a minimal impact on the final horizontal displacement, with differences of less than 1 mm. This is likely due to the supporting effect of the grouting pressure. However, rapidly setting grout significantly reduces the maximum horizontal displacement of the strata immediately after the shield tail passes and grouting occurs, while also accelerating the rebound of horizontal displacement during the grouting process.

Figure 12 shows the final surface vertical displacement for different grout initial setting times. Clearly, as the grout initial setting time increases, the surface vertical displacement significantly increases, and the rate of surface vertical displacement development also accelerates.

Figure 13 shows the distribution of surface vertical displacement along the tunnel axis when the shield cutterhead reaches Y = 72 m (at this point, the grouting layer behind the 24th ring of lining is fully hardened, and the shield tail is at Y = 61.2 m). Clearly, for the four models with different grout initial setting times, the surface settlement reaches its maximum approximately two ring positions (Y = 46.5 m) before the grouting layer is fully hardened. The slower the grout hardens, the worse its control effect on surface settlement. When the grout initial setting times are 0 h, 1 h, 2 h, and 4 h, the maximum surface displacement above the tunnel axis is −4.13 mm, −4.97 mm, −5.83 mm, and −6.77 mm, respectively.

Moreover, the development of surface vertical displacement follows the same pattern in all four models. Before the cutterhead arrives, the surface exhibits slight heave due to the pressure from the cutterhead. As the cutterhead arrives, the surface vertical displacement rapidly develops into settlement. The rate of settlement growth reaches its maximum at the shield tail, with the surface settlement at this point being 0.5 times the maximum settlement value. After reaching the maximum settlement value, the surface settlement begins to rebound.

The geotechnical tests conducted at the metro construction site indicate that the mudstone strata encountered by the tunnel exhibit a high Poisson’s ratio and a low elastic modulus. The low compressibility of the mudstone strata suggests that the error introduced by neglecting the surface vertical displacement caused by strata compression is acceptable. Additionally, as mudstone is a relatively weak rock formation, it is reasonable to approximate it as soil for calculation purposes.

For the Chongqing Metro Line 27, the shield excavation radius is 4.43 m. In this study, the tunnel axis depth is uniformly taken as 15 + 4.43 m. Therefore, the calculated correction factor is 2.43. After analyzing the radial shrinkage data of the segments on-site, it was found that the shrinkage values of the segments are concentrated between 3.5 and 4.5 mm. Four frequently occurring values were selected for the calculation, namely,: 3.6, 3.9, 4.1, and 4.4 mm. Additionally, to investigate which location, between the cutterhead and the shield tail, is more suitable as the point of y = 0 for the calculation, two calculation methods—where y ′ = y and y ′ = y + l —were employed to compute the surface vertical displacement.

Based on this, Figure 14 presents the surface vertical displacement along the y-direction calculated using the two methods, along with the surface deformation results obtained from numerical simulations at the same excavation position. Due to the size limitations of the numerical model, the simulation results do not cover the entire Y range. However, it can be observed that the surface vertical displacement range shown by the numerical simulation is smaller. Compared to the theoretical calculation results, surface vertical displacement in the numerical simulation begins to occur only when the cutterhead is closer, and it stabilizes earlier. This also indicates that the missing numerical simulation data has a negligible impact on the analysis in this study.

Moreover, it is evident from the numerical simulation results in the figure that surface vertical displacement significantly occurs at the shield tail. When using the method y ′ = y + l ., when calculating with the original model, it is more reasonable to add the shield shell length to the value of y. Furthermore, it is observed that the numerical simulation results not only change more rapidly but also show a noticeably larger maximum settlement value. The on-site monitoring data shown in Figure 11 indicates that there is indeed slight heave on the surface before the cutterhead reaches the location. At the same time, the on-site monitoring results are slightly larger than the numerical simulation results, which suggests that the theoretical calculation of the maximum deformation is closer to the actual situation. However, the numerical simulation can predict the surface heave in front of the cutterhead, which highlights the advantage of the simulation.

Figure 15 shows the surface vertical displacement along the x-direction at different locations using two calculation methods, compared with the numerical simulation results. (a), (b), (c), and (d) show the surface vertical displacement in four different states: when the cutterhead reaches the location, when the shield tail passes and grouting occurs, 20 h after grouting, and after the surface vertical displacement stabilizes. It is evident that the numerical simulation results change more rapidly over time. Clearly, the numerical simulation results exhibit a faster rate of change. Additionally, surface vertical displacement calculated using methods y ′ = y and y ′ = y + l gradually converge over time. However, before the deformation stabilizes, the results obtained using y ′ = y + l generally fall between those of the other two methods.

Although this work combines analytical solutions with numerical modeling to quantify how grout setting and hardening affect construction-stage surface settlement in shallow-buried shield tunneling in mudstone strata, and identifies the shield tail as the critical longitudinal reference location for deformation development, several simplifications were adopted to keep the framework computationally efficient and suitable for systematic parametric evaluation. The main boundaries of applicability are summarized as follows.

Time-scale representation and long-term effects. The modeling focuses on construction-stage responses associated with synchronous grouting and early-age grout behavior, and it adopts a construction-step, quasi-static representation to efficiently evaluate the influence of grout setting time (0 h–4 h) (Wang et al., 2025). This approach captures the stabilized response relevant to the short time window of interest, but it does not reproduce detailed short-term transients immediately after grout injection and pressure dissipation, nor does it address long-term rheological effects of mudstone (e.g., creep) that may become relevant over months to years (Wang et al., 2020; Yan et al., 2025). Fully time-dependent analyses and rheology-aware constitutive descriptions would be valuable when transient evolution or longer time scales are the focus.