Following the work of Randolph et al. (1994), the cone tip resistance q_c is related to the limit pressure p_lim during the spherical cavity expansion. When using a more practical complex constitutive model, numerical analysis is often required for simplified calculations. A spherical cavity expansion was modelled in the Plaxis 2D to simulate the penetration process, following similar procedures described by Xu and Lehane (2008) and Suryasentana and Lehane (2014). Considering the symmetry of spherical cavity expansion, the model dimension is x = 12 m, y = 24 m. According to the studies of Suzuki (2015) and Xu (2007), there are no boundary effects. The initial radius of the spherical cavity a₀ is 0.1 m. A positive volume strain of 10% is applied to the cavity cluster in multiple stages to simulate the gradual expansion. The construction phase of the model is divided into 21 steps. Before calculation, 7 nodes and 8 stress points around the cavity wall are selected as output. Finally, the selected nodes and stress points data are averaged separately to plot the relationship between the limit pressure p_lim and excess pore pressure Δu with the normalized displacement a/a₀. The cone tip resistance q_c is obtained from the following relationship: q c = p lim + 3 [ c ′ + ( p lim − u ) ] tan φ ′ .

In order to correct the relationship proposed by Suzuki (2015), Table 1 listed 15 analysis cases to study the influence of model parameters E 50 ref and φ^’ on the limit pressure p_lim using the HSS model. The influence of the over consolidation ratio is not considered (OCR = 1). The permeability coefficient of soil is 0.02 m/d (2.3×10^-7 m/s), according to the geological prospecting report.

The cone tip resistance q_c can be obtained by substituting the limit pressure p_lim of each group into the formula. Therefore, the normalized cone tip resistance Q, E 50 ref /p’₀ and φ^’ are drawn in Figure 2. The relationship between logQ and log( E 50 ref /p’₀) is almost linear in logarithmic coordinates. According to the results proposed by Suzuki (2015) based on the HS model, the normalized cone tip resistance Q relation based on the HSS model can be obtained as follows:

The influence factors such as particle composition and transport process lead to different soil properties. Therefore, there is inherent variability in soil, independent of the knowledge state of geotechnical properties, and will not decrease with the increase of knowledge (Phoon and Kulhawy, 1999; Wang and Cao, 2013; Wang et al., 2023). E 50 ref is a continuous variable and must be non-negative because of its physical meaning. Therefore, geotechnical parameters are often modeled as logarithmic normal random variables (Wang and Cao, 2013). Assume that E 50 ref is a lognormal random variable with a mean μ and a standard deviation σ, expressed as

in which μ_N and σ_N are the mean and standard deviation of ln E 50 ref , z is a normal random variable with a mean of 0 and SD of 1.

Eq. (1) can be rewritten in a log-log scale as:

in which lnQ is normalized cone tip resistance in a log scale; a, b is the coefficient, ϵ is a Gaussian variable representing the transformation uncertainty.

Combining Eqs. (2) and (3) lead to:

For the given prior knowledge and CPTU data, the probability density function of E 50 ref can be expressed as:

in which P( E 50 ref |μ,σ) is the conditional probability of a given set of μ and σ. Because E 50 ref is lognormally distributed, it can be expressed as:

Using Bayesian formula, P(μ,σ|Data, Prior) can be expressed as:

in which K is a normalized constant independent of μ and σ; Data = {lnQ_i , i=1, 2, …, n_s } is a set of CPTU test data; P(Data|μ,σ) is the likelihood function which is expressed as:

P(μ,σ) is assumed to be a uniform distribution, expressed as:

Substitute Eqs. (6) - (9) into Eqs. (5) to obtain the posterior probability density function of E 50 ref calculated according to the approach proposed by Wang and Cao (2013).

The equivalent samples of E 50 ref reflect the combined information of a prior knowledge and corresponding test data. Three groups of mean shown in Figure 7A are selected to study the influence of mean μ. All of them are uniform distribution, of which the second group has been used in Section 4.1 and will serve as the basis for the research. The Bayesian inversion results are shown in Figure 7B. As seen in the histogram, the final results obtained by μ∈[5.0, 7.0] and μ∈[3.9, 8.0] are closer, peaking at 6.5 MPa and 5.5 MPa in the center of the interval, respectively. The results obtained in the maximum interval range μ∈[1.6, 15.0] differ greatly from the two ranges mentioned above, reaching a peak at 13.5 MPa in the center of the interval. At the same time, the sample distribution ranges of μ∈[5.0, 7.0] and μ∈[3.9, 8.0] are consistent. Samples with a smaller prior range μ∈[5.0, 7.0] are more concentrated, and the sample size at the peak exceeds 7000. The large prior range of the mean value makes the final result deviate greatly from the true value of E 50 ref .

The samples of sensitivity analysis were also used to estimate the probability density function and cumulative distribution function of E 50 ref , as shown in Figure 7C, D. The gray dashed line in the figure is the peak value, from which the results are consistent with the above analysis. Through the equivalent samples estimation, the final results of E 50 ref obtained from μ∈[5.0, 7.0] and μ∈[3.9, 8.0] are 6.2 MPa and 6.9 MPa, with a deviation of 3.1% and 7.8% from 6.4 MPa in laboratory test, respectively. However, the results obtained from μ∈[1.6, 15.0], with a deviation almost doubled, significantly deviating from the acceptable error range of the parameter inversion results. The standard deviation of the three groups is 1.9MPa, 2.1MPa and 2.7MPa, respectively. The range of the mean value of E 50 ref in this inversion is a relatively subjective judgment obtained by statistical analysis of a small amount of data. Its value range is derived from the statistics of E 50 ref test values of soils in different regions. Therefore, better inversion results were obtained when the selected test values μ∈[3.9, 8.0] were chosen for the silty clay in each region. Narrowing this range gives more accurate results, such as μ∈[5.0, 7.0]. And when μ∈[1.6, 15.0], the increase in scope leads to ambiguity of information, resulting in unreasonable results. Therefore, a reasonable range of prior knowledge of soil parameters must be selected in practical application. Statistics of silty clay values in specific areas will greatly increase the reliability of the inversion results.

Similar to the research of μ, the three groups of standard deviation are selected to explore the influence of σ, as shown in Figure 8A. Compared to the μ, the results of σ∈[0.5, 2.0], σ∈[0.5, 3.0] and σ∈[0.5, 8.0] have little difference in Figure 8B. The peaks are located close to each other, at 7.5 MPa, 6.5 MPa and 6.5 MPa in the center of the interval, respectively. Similarly, as the selected range of σ decreases, the obtained samples become more concentrated and the probability density function becomes steeper, reflecting that the accurate selection of the range of prior knowledge has a greater influence on the posterior distribution.

The equivalent samples are plotted as probability density function and cumulative distribution function, as shown in Figure 8C, D. It can be determined from the peak position that the change of σ range has less effect on the peak of E 50 ref probability density function than μ. However, the change of σ range has a certain influence on the shape of probability density function, which gradually transitions from normal distribution to lognormal distribution. The estimated E 50 ref for the three ranges is 6.9 MPa, 6.9 MPa, and 6.8 MPa, with deviations from laboratory test 6.4 MPa of 7.8%, 7.8%, and 6.3%, respectively. The standard deviation of samples estimation is 1.6MPa, 2.1MPa and 4.6MPa, respectively. From another perspective, the larger σ range makes the sample more discrete and changes more significantly than μ.

PLAXIS was used to establish the numerical model of the Yellow River tunnel passing through the south embankment in this paper. To simplify the calculation, the following assumptions were made in the modelling process: (1) The silty clay stratum is homogeneous and isotropic; (2) Regardless of the longitudinal slope of the tunnel, the tunnel is laid horizontally; (3) The water level of the Yellow River is constant, and the seepage is neglected.

The model dimension is 360 m×150 m×100 m, the overburden depth of the tunnel is 40 m, and the shallowest section is 30 m below the Yellow River. A support pressure of 0.5 MPa is set in front of the tunnel face and increases with depth by 0.01 MPa/m. The tunnel advances in the Y-direction with every 10 m (5 rings) for a construction stage. The eastern line is constructed before the western line. The numerical model is shown in Figure 9.

The HSS model and elastic model were assumed for soil, shield shell and segment, respectively. Bayesian inversion result was used for the value of E 50 ref . The value of E oed ref and E ur ref was obtained from the empirical relationships, expressed as E oed ref = E 50 ref , E ur ref = 3 E 50 ref . And the other parameters were taken from the laboratory test, such as consolidation drainage test, standard consolidation test and resonant column test, or the default values. The Young’s Modulus E of segment is decreased to 80% of its initial value of 35.5 GPa, which is based on the tunnelling experience in China, to take into account the impact of joint on the stiffness of tunnel lining (Feng et al., 2011; Xie et al., 2016). The soil parameters are listed in Table 2.

To further study the deformation caused by shield tunnelling, three monitoring points, D₁, D₂ and D₃, were set up at the top of the embankment section. Meanwhile, two deep settlement monitoring points, S₁ and S₂, were set at 30 m and 34 m below point D₃. The site monitoring point layout is also shown in Figure 9.

The monitoring and simulation results of the D₁, D₂ are shown in Figure 10A. It can be seen that the embankment settlement can be divided into three stages, which are the initial stage, rapid settlement stage and stable stage. According to the monitoring results, the support pressure applied in front of the tunnel face caused the monitoring sites D₁ and D₂ to heave. Compared with the field measurement, the simulation results of shield tunnelling did not show the phenomenon of surface heave. And the settlement started at the beginning of the tunnelling, but the settlement was almost negligible. When the shield tunnelling reached the vicinity of the monitoring point, both the simulated and monitored embankment settlement rates accelerated significantly and a large embankment settlement appeared. In this stage, the maximum settlement predicted by the monitoring and simulation is 5.5 mm and 6.5 mm, respectively, accounting for 65% and 90% of the total settlement. The difference between the two settlement is small, but the settlement proportion is large. This is because in the field measurement, there is the re-consolidation settlement after soil disturbance, and the settlement of this part is second only to the settlement generated by construction disturbance in the project. The rapid settlement stage (Stage 1) at point D₂ is longer, with a maximum settlement of 12 mm, accounting for more than 90% of the total embankment settlement. This indicates that the advanced construction of the east line greatly influences the stacking of the west line settlement, and the final settlement of the west line is greater than the east line. The settlement rate of monitoring points D₁ and D₂ in the rapid settlement stage is about 0.8 - 1.3 mm/d. Figure 10A also shows the results obtained from the non-inverted model, which can be seen to be larger than the measured and inverted results. It indicates that the non-inverted parameters underestimate the stiffness of the soil, resulting in a certain deviation in the calculation results. One issue should be addressed in Figure 10A. When the east line was excavated, the boundary settlement was about -6 mm. rather than 0 mm, as we find out in most field monitoring. The mesh size and the chosen soil model are both to cause this phenomenon. A smaller mesh size is typically required for the advanced soil model to obtain convergence results. However, a small mesh size will reduce computational efficiency and is unsuitable for large models. It is noted that the total settlement of monitoring point D₂ above the west line was higher than D₁ above the east line. This is because when the axis distance between the two tunnels is relatively close, the soil in the middle area of the two tunnels will be affected by two large-diameter shield construction disturbances. In this model, the results of measured maximum settlement at D₁ and D₂ are 8.5 mm and 13 mm, respectively. The numerical results are 7.3 mm and 10.5 mm, respectively.

The comparison of field measurement between monitoring points D₃, S₁ and S₂ in Figure 10B shows that the settlement trend of D₃ was roughly the same as that of S₁ and S₂. Among them, the settlement of the monitoring point S₂ closer to the tunnel (about 3 m) is more extensive, and the maximum settlement reaches -13.5 mm. The value and settlement rate of this monitoring point is greater than those of the other two. From the simulation results, the settlement changes of S₁ and S₂ are the same at the beginning of the tunnel excavation, and the amount of change is small. The settlement at D₃ is about 2 mm larger than the deep settlement. The settlement of the three monitoring points grows quickly when the tunnel is close to an embankment, and the deep settlement is more noticeable than the surface settlement. The maximum settlement of S₂ is slightly larger than S₁, and the ground settlement is the smallest. In terms of the total settlement of each stage, the measured monitoring point S₁, which is located around 7 m above the tunnel, is close to the monitoring point D₃ at the top of the embankment. The maximum settlement is -10.6 mm and -11.4 mm, respectively.

At the same time, the simulated settlement curve is about 10 rings ahead of the field measurement. This is because in numerical simulation, considering the horizontal and homogeneous soil layer, the set support pressure is a constant. In actual construction, due to the dynamic adjustment of construction parameters, the displacement near the tunnel surface caused by construction has been well controlled, and the affected range is smaller than the simulation results, resulting in lag in the measured results. The corresponding conclusion can also be drawn from the fact that the final simulation prediction result is larger than the actual measurement. In general, comparing numerical simulation and field measurement, it can be seen that the simulation results agree with the measurements, proving the viability of the method used to obtain the HSS model parameters from the CPTU data.

The soil disturbance caused by tunnel construction can be significantly reduced by using reasonable support pressure, avoiding engineering risks such as splitting and roof falling in shallow underwater soil sections and river water backflow. The actual support pressure set in the Yellow River tunnel project is mostly between 0.3 MPa - 0.6 MPa. Therefore, six analysis cases of 0.3 MPa, 0.35 MPa, 0.4 MPa, 0.45 MPa, 0.5 MPa and 0.55 MPa were selected.

Figure 11 shows the maximum riverbed settlement under different support pressures. When the support pressure is 0.3 MPa, the influence of soil displacement in front of the tunnel face has caused a large riverbed settlement, the maximum settlement is -8.9 mm. When pressure increases to 0.45MPa, the riverbed settlement is well controlled, and the maximum settlement is -2.0mm. Before 0.45MPa, the decrease rate was almost linear, and after 0.45MPa, the decrease rate slowed down. It can be seen that when the support pressure increases to a certain extent, the influence of support pressure on settlement control is limited. For this model, the optimal supporting pressure is between 0.45MPa and 0.5MPa.

The longitudinal slope of the shield tunnel is mainly controlled within 5^° in the design. However, the longitudinal slope may exceed the design value due to the complexity of geological conditions and the uncertainty of shield attitude control. Cheng et al. (2021) designed the longitudinal slope angle as 15^° to capture the tunnel failure mode more obviously in the small-size model test and verified it with numerical simulation. In most projects, the maximum longitudinal slope angle in practical application is usually less than 10^°. Therefore, +10^°, +5^°, 0^°, -5^°, -10^° were selected for comparative study.

Figure 12 shows the failure mechanism of the tunnel face at different slopes. When the slope angle of the tunnel is +10^°, the soil deformation extends to the ground, resulting in an extensive range of settlement trough. When the slope angle is +5^°, the range of settlement trough is smaller. The failure mechanism of a horizontal tunnel is similar to +5^°. In all cases of shield downslope, the deformation effects did not reach the surface. This condition is mainly determined by the horizontal projection area of tunnel on the ground. It is generally understood that the larger the projection area is, the larger the impact area of the tunnel excavation on the riverbed, which will lead to increased settlement.

Figure 13A shows the transverse riverbed settlement of different slopes at which the buried depth is 2D. The transverse riverbed settlement of the tunnel reaches its maximum value at the middle line of the tunnel. With the tunnel slope from +10^° to -10^°, the maximum transverse riverbed settlement is -25.8 mm, -20.8 mm, -17.0 mm, 13.4 mm and -11.1 mm, respectively. It is obvious that the uphill section will cause greater riverbed settlement than the downhill section.

Buried depth affects the tunnel construction deformation. Under the same condition, different buried depths will cause different stratum displacements. Generally, the shallowest buried depth of the tunnel should be greater than 0.7D (D refers to the outside diameter of the tunnel). And the buried depth of the Jinan Yellow River Tunnel is between 11.2 - 42.3m, about 0.7D - 2.8D. Based on this, six working conditions were selected: 0.7D, 1.0D, 1.5D, 2.0D, 2.5D and 3.0D.

Figure 13B shows the transverse riverbed settlement of different buried depths. Under the same conditions, when the buried depth is less than or equal to 2.0D, the riverbed settlement increases with the increase of the buried depth. Under 0.7D, the supporting pressure of 0.3 MPa has exceeded the sum of soil and water pressure. Most areas are uplifted, and the maximum uplifted is about 5.0mm. The transverse riverbed settlement of the tunnel at 1.0D is small, about -1.6mm at most. With the increase of tunnel buried depth, the maximum riverbed settlement reached -8.5mm at 2.0D. It is worth noting that the above phenomenon was studied only by changing the buried depth. In other words, the greater the buried depth, the smaller the support pressure ratio.

It can also be seen that when the tunnel buried depth is 2.5D and 3.0D, the maximum riverbed settlement no longer increases with the increase of buried deep, and the law is precisely the opposite. This is because when the buried depth is greater than 2.0D, the soil arch effect is generated in the upper soil, preventing the formation deformation caused by tunnel excavation transmitting to the riverbed, and the riverbed deformation is reduced.

Construction of underwater tunnels is challenging because they must withstand greater water pressure than conventional tunnels. The level of water not only affects the structural design of segments and waterproofing and plays a crucial role in the selection of construction parameters (Du et al., 2021). The water condition also restricts the implementation of various pretreatment measures. Before the Wuhan Yangtze River Tunnel, the maximum water pressure of the large-diameter shield built on the Huangpu River in Shanghai did not exceed 0.45 MPa. Subsequently, the maximum water pressure is constantly refreshed with the vigorous development of large-diameter shield tunnels.

Combined with the previous research, when the support pressure is 0.3 MPa, the water level of 3 m, 5 m, 7 m, 10 m, 20 m was selected for analysis. The calculation results are shown in Figure 13C. The riverbed settlement increases with the water level rise. This is because the change in water level in the model is equivalent to the uniform load of different sizes applied on the riverbed. In contrast, the high-water level represents that the model receives a greater vertical load, naturally increasing the riverbed settlement. The maximum transverse settlement caused by the 20 m water level is -54.1mm and -54.7mm, respectively. It indicates that the support pressure adopted cannot support the tunnelling, so the results are not shown.