Dynamic Parameter Calibration of an Analytical Model to Predict Transient Groundwater Inﬂow into a Tunnel

: The use of either analytical or numerical models of groundwater inﬂow into a tunnel, ignoring the excavation process, will result in inaccurate prediction. More researchers have begun to study prediction methods with temporal and spatial variables, in accompaniment with the tunnel excavation progress. The hydrogeological parameters signiﬁcantly affect tunnel inﬂow prediction; thus, it is always worthy to carry out an inversion analysis when the inﬂow rate observation is available for the excavated part of a tunnel. Based on a transient analytical model and the Trust Region Reﬂection (TRR) algorithm, this paper proposes a Dynamic Parameter Calibration (DPC) method to sequentially optimize parameters for a tunnel that is divided into several sectors. The results of two case studies indicate that the ﬁtting effects by DPC are signiﬁcantly improved compared with the empirical trial-and-error method, indicating good optimizations. The proposed scheme can conduct calibration simultaneously for multiple tunnel sectors and also several optimizations for the same sector. The parameter optimization results reﬂect the lithologic heterogeneity of different strata. Parameter sensitivity analysis proves that the hydraulic conductivity K has a greater inﬂuence on water inﬂow calculation than the speciﬁc storage coefﬁcient S .


Introduction
The study of tunnel water inflow has always been a hot and difficult point of concern.In the process of tunnel construction, accurate prediction of water inflow can effectively reduce the risk of water inrush disaster and ensure the safety of tunnel construction and excavation progress.Scholars have systematically summarized the calculation methods of water inflow and introduced in detail the theoretical principles, assumptions, advantages and disadvantages, and applicable conditions of the existing calculation methods of tunnel water inflow [1,2].When using analytical methods and numerical models to calculate the water inflow of tunnels, most of the model-makers ignore the excavation process, which leads to inaccurate prediction of water inflow [3].Based on convolution and superposition principles, Perrochet et al. proposed an analytical model of transient water inflow during tunnel excavation in a subvertical multilayered heterogeneous stratum [4,5].Golian et al. used MODFLOW to simulate the inflow of groundwater as the tunnel boring machine (TBM) gradually advances and compared the predicted value of water inflow with the observed value.The results showed that the presented method could satisfactorily predict the inflow rates as the TBM advances [6].
Water 2023, 15 Hydrogeological parameters are important in the prediction of tunnel water inflow.According to the long-term observation data of groundwater or the observation data of a pumping test, the hydrogeological parameters are determined and the hydrogeological conditions are understood; this is the so-called inverse problem [7].Various inversion methods have been proposed in the last three decades, from the trial-and-error manual calibration to the current complex automatic data assimilation algorithms [8].Perrochet et al. proposed an analytical equation to calculate the water inflow into a tunnel; they also calibrated the hydrogeological parameters by trial-and-error method to fit against inflow measurements [5].Some researchers have calibrated and optimized the model parameters by using the groundwater inverse programs UCODE and PEST based on nonlinear least squares theory [9][10][11].Xu et al. established a large-scale finite element model for a diversion tunnel and optimized the hydrogeological parameters of the aquifer and faults, and the maximum water inflow of the diversion tunnel was predicted in a small-scale model considering tunnel construction [12].The influence of various parameters on the calculation of water inflow were analyzed through different methods; these parameters include the curtain grouting reinforcement cycle, the thickness and hydraulic conductivity of tunnel lining and grouting rings, and so on [13,14].Ma et al. used the measurements of the tunnel water inflow and the pore water pressure to calibrate the hydraulic conductivity, then the optimized parameters were assigned to the program to calculate the water inflow for another sector, and the comparison of the calculation and the observation proved that the Differential Evolution Algorithm performs well in the seepage inverse analysis [15].Scholars have calibrated the parameters of the water inflow model by combining a variety of optimization algorithms and techniques, and they significantly improved the accuracy of water inflow prediction [16,17].
The permeability of a rock mass is a parameter with a very wide range of variability due to a large number of influencing factors, especially the spatial heterogeneity of the discontinuities for fractured rock [18].Sometimes, even despite similar filtration parameters of the aquifer, groundwater flows may be characterized by significant differences in terms of flow paths and flow characteristics of individual streams of water [19].Different factors may have different effects on the calculation of groundwater inflow into a tunnel.Sensitivity analysis helps to identify important parameters that are relatively sensitive to a model, reduce the parameter dimension, and provide support for model parameter optimization and uncertainty analysis [20].Based on the comprehensive study of several hard rock tunnels, the influence of geological parameters on the prediction of tunnel water inflow were investigated, such as Q value, faulting, rock type, and width of weakness zones [21,22].The literature demonstrated that sensitivity analysis is a feasible and effective method to estimate the influence of different hydrogeological parameters on the calculation of water inflow, and some results show that hydraulic conductivity K is the most sensitive parameter which has the greatest influence on the calculation result of water inflow [23,24].Xia et al. assigned different K (hydraulic conductivity)-S s (specific storage) combinations for a high-K zone and simulated the daily flow rate at the intersection of a tunnel and that high-K zone.The sensitivities of the time-varying drainage process to the K-Ss combinations were compared, and the results showed that the discharge process was very sensitive to the hydraulic parameters; the maximum flow rate increased with the value of K as well as S s .By analyzing the linear relationship between the parameters and the decay exponent of flow rate, it was found that the variability of the K value had a more significant effect on the flow rate [3].However, research is still lacking on the automatic calibration of the hydrogeological parameters by some optimal algorithm to achieve a good fitting of calculation versus observation of groundwater inflow rate into a tunnel.
In this paper, an optimization scheme was proposed based on an analytical model of groundwater inflow by Perrochet, and the trust region reflection algorithm was applied for calibration.Taking the Modane exploratory tunnel in France and the Xiema tunnel in China as examples, the key hydrogeological parameters were dynamically calibrated and optimized considering the excavation process of the tunnels.The influence of different parameters on the prediction by the analytical model is investigated according to sensitivity analysis, and the uncertainty of parameters is studied as well.

Perrochet Model and Dynamic Parameter Calibration Process
By using convolution and superposition principles, Perrochet et al. proposed an equation to calculate the total water inflow Q (the dimension is L 3 /T) of the tunnel excavated sectors with time t in a subvertical multilayered heterogeneous stratum (the exposed rock mass of the tunnel is divided into N sectors, as shown in Figure 1) [5]: where t is the time (T), t i is the moment when the ith sector of the tunnel is excavated (T); the subscript i indicates the ith sector of the tunnel, a positive integer of 1 ≤ i ≤ N; L is the length of each sector (L); x is the length of the tunnel excavation sector along the axis (L), 0 ≤ x ≤ L; v is the tunnel drilling speed (L/T); K is the hydraulic conductivity (L/T); S is the specific storage coefficient (1/L); s is the water level drop depth (L); r is the radius of the tunnel (L); and H(L i − x) is the Heaviside step-function.
Water 2021, 13, x FOR PEER REVIEW 3 of 14 for calibration.Taking the Modane exploratory tunnel in France and the Xiema tunnel in China as examples, the key hydrogeological parameters were dynamically calibrated and optimized considering the excavation process of the tunnels.The influence of different parameters on the prediction by the analytical model is investigated according to sensitivity analysis, and the uncertainty of parameters is studied as well.

Perrochet Model and Dynamic Parameter Calibration Process
By using convolution and superposition principles, Perrochet et al. proposed an equation to calculate the total water inflow Q (the dimension is L 3 /T) of the tunnel excavated sectors with time t in a subvertical multilayered heterogeneous stratum (the exposed rock mass of the tunnel is divided into N sectors, as shown in Figure 1) [5]: where t is the time (T), ti is the moment when the ith sector of the tunnel is excavated (T); the subscript i indicates the ith sector of the tunnel, a positive integer of 1 ≤ i ≤ N; L is the length of each sector (L); x is the length of the tunnel excavation sector along the axis (L), 0 ≤ x ≤ L; v is the tunnel drilling speed (L/T); K is the hydraulic conductivity (L/T); S is the specific storage coefficient (1/L); s is the water level drop depth (L); r is the radius of the tunnel (L); and H(Li − x) is the Heaviside step-function.In this paper, the key hydrogeological parameters, namely, the hydraulic conductivity K and specific storage coefficient S in Equation (1) are calibrated for optimization.Other parameters ti, Li, si, vi, ri are set as known variables, their values have been specified right after the tunnel division.The measurement of groundwater inflow rate in the tunnel is denoted as Q obs , which are M dimension time series data.In general, M is much larger than N; in other words, the quantity of water inflow measurements is far more than the number of tunnel sectors.Qi obs is referred to as the time series data of the inflow observation in the ith sector, and mi is the number of measured inflows in the ith sector.In this paper, the key hydrogeological parameters, namely, the hydraulic conductivity K and specific storage coefficient S in Equation (1) are calibrated for optimization.Other parameters t i , L i , s i , v i , r i are set as known variables, their values have been specified right after the tunnel division.The measurement of groundwater inflow rate in the tunnel is denoted as Q obs , which are M dimension time series data.In general, M is much larger than N; in other words, the quantity of water inflow measurements is far more than the number of tunnel sectors.Q i obs is referred to as the time series data of the inflow observation in the ith sector, and m i is the number of measured inflows in the ith sector.
The flow chart of dynamic parameter calibration is shown in Figure 2. When the process goes to the ith tunnel sector, the following calculation will be carried out step by step as described below: 1.
Initial values assignment: The initial values of K i in and S i in are assigned to the parameters of hydraulic conductivity and specific storage for sector i.

2.
Water inflow calculation: This step is denoted as "Perrochet's Equation" in the flow chart in Figure 2. A MATLAB code script for inflow calculation has been written according to Equation (1).After the initial values of K i in and S i in are assigned, the code runs, then the inflow calculation is obtained for the ith sector, which is referred to as Q i cal .It is also a time series including m i values, exactly the same dimension as its paired observation Q i obs .

3.
Parameter calibration: Once the Q i obs and Q i cal are available, the trust region reflection algorithm will be applied to calibrate and optimize K i and S i .The number of iterations in the optimization depends on the stopping condition of the solver, namely, Tolerance, referred to as Tol.When the tolerance is lower than the threshold ε (ε = 1 × 10 −6 ), the stop criterion is met, then the iteration of the solver will be terminated.In this step, the proposed optimal algorithm will run repeatedly until the stop criterion is met; at that time, the optimized hydraulic conductivity K i op and specific storage coefficient S i op are determined as well.4.
Results output: The optimized K i op and S i op will be stored in a database to form parameter sets, while the Coefficient of Determination is calculated as: where SSE is the sum of squared errors between observations Q i obs and calculations Q i cal , and SST stands for the sum of squares.This coefficient is often used to evaluate the fitting effect of optimization.The time-dependent water inflow curves Q i obs -t and Q i cal -t will be plotted as well, to visually display the fitting effect between calculation and measurement.

5.
Next calculating sector: As the calibration goes to the (i + 1)th sector, the entire process will be conducted again, as shown in the flowchart in Figure 2. The optimized parameters of the ith sector will be used for water inflow calculation in step (2), and the parameters of the (i + 1)th sector will be calibrated in step (3).
Therefore, the calibration process is carried out dynamically with the tunnel excavation.Thus, the parameters of hydraulic conductivity K and specific storage coefficient S for each sector can be optimized accordingly.
Considering the uncertainties of parameters, we made a slight change to the above dynamic process to realize simultaneous optimization for multiple tunnel sectors.As a result, the same parameters of a sector will be optimized several times.Let k be the number of tunnel sectors in each optimization step, which should be a positive integer greater than or equal to 1.When the calibration process goes to the ith sector, the total k tunnels will be calibrated including the i, i − 1,... i − k + 1 sectors, while the parameters of the i − k to 1 sectors are unchanged, becoming the optimized values by the last calibration.For example, if k = 2, when the process is going to the ith sector, parameters for both the i and i − 1 sectors will be calibrated at this stage, while no parameters will be changed for the i − 2 and all the preceding sectors.In particular, if k = 1, that is, only one sector is being optimized at a time, then the parameters from i − 1 to 1 remain unchanged.Specifically, if i < k, all excavated sectors should be calibrated.For instance, if k = 3, when i = 1 and 2, the first and second sector must be calibrated; when i ≥ 3, there are three sectors to be calibrated in each stage until the final Nth stage.Figure 3 illustrates this example.It is worth noting that a reference range needs to be set for both hydraulic conductivity K and specific storage coefficient S, which is also the bounded constraint specified in optimization algorithm.Specifically, the optimal value interval of K is set to be 10 −10 -10 1 m/d, while it is 10 −10 -10 −1 m −1 for S.
The trust region reflection algorithm was used for parameter calibration in this paper, referred to as TRR, which is a nonlinear least squares optimal algorithm.The basic idea of TRR is to approximate the objective function f(x) with the quadratic function q(s), which is the mapping of the function f(x) on the neighborhood N around the current point x.x is   It is worth noting that a reference range needs to be set for both hydraulic conductivity K and specific storage coefficient S, which is also the bounded constraint specified in optimization algorithm.Specifically, the optimal value interval of K is set to be 10 −10 -10 1 m/d, while it is 10 −10 -10 −1 m −1 for S.
The trust region reflection algorithm was used for parameter calibration in this paper, referred to as TRR, which is a nonlinear least squares optimal algorithm.The basic idea of TRR is to approximate the objective function f(x) with the quadratic function q(s), which is the mapping of the function f(x) on the neighborhood N around the current point x.x is It is worth noting that a reference range needs to be set for both hydraulic conductivity K and specific storage coefficient S, which is also the bounded constraint specified in optimization algorithm.Specifically, the optimal value interval of K is set to be 10 −10 -10 1 m/d, while it is 10 −10 -10 −1 m −1 for S.
The trust region reflection algorithm was used for parameter calibration in this paper, referred to as TRR, which is a nonlinear least squares optimal algorithm.The basic idea of TRR is to approximate the objective function f (x) with the quadratic function q(s), which is the mapping of the function f (x) on the neighborhood N around the current point x.x is a vector with bounded constraints, and the neighborhood N is called a trust region [25].In order to solve the trusted region problem, the test step size s is calculated by minimizing  3) and (4).If f (x + s) < f (x), the current point x is updated to x + s.This step is considered a success, and the trust area will continue to the next step.Otherwise, if the step is not successful, then x remains unchanged and the region N will be reduced to the next step [25].
where g is the gradient of f (x) for the current x, H is a symmetric matrix of second derivatives, D is a diagonal scaling matrix, ∆ is the trust-region radius > 0, and is the second norm [25].
The general process of unconstrained minimization based on the trust region is as follows: Step 1, to construct a two-dimensional trust region subproblem; Step 2, to solve Equation (4) to determine the trial step s; Step 4, to adjust ∆.Repeat these four steps until convergence to obtain the optimized parameter values.
The trust region reflection algorithm can be realized by using the lsqcurvefit function in the MATLAB (R2021b) software.

Case Studies
In this section, two case studies of two respective tunnels will be investigated.The transient groundwater inflow for both tunnels were calculated by an analytical model proposed by Perrochet et al. [5].In addition, the hydrogeological parameters, namely, the hydraulic conductivity K and specific storage coefficient S, were calibrated by trial-and-error method based on that model.Herein, we will calibrate those parameters by the aforementioned dynamic parameter calibration (DPC) method and compare the results of groundwater inflows and optimized parameter values between the two calibration methods.

Modane Exploratory Tunnel, France
The Modane exploratory tunnel is located in the French Western Alps and is part of the geologic investigation program of the Lyon-Turin high-speed railway project.This tunnel is excavated in metamorphic rocks belonging to a Permian-Triassic sequence of the so-called "Brianzonese Zone" [5].Perrochet et al. [5] applied the calculation model of transient water inflow to divide the 760 m tunnel into 10 sectors according to the permeability of the rock stratum based on the collected detailed data.The hydraulic conductivity K and specific storage coefficient S were inverted by the trial-and-error method.The analytical calculation results were in good agreement with the measured values.The application of this approach for tunnel discharge forecasting should therefore always be associated with a reliability analysis of the hydrogeologic conceptual model [5].
Here, we applied the dynamic parameter calibration (DPC) to optimize parameters K and S; the tunnel sector division is exactly the same as that in the Reference [5], and the parameterization keep same except K and S. The comparison of the time-varying inflow rate between the trial-and-error method and DPC is shown in Figure 4a.The calculated inflow by two methods both presented the dynamic characteristics of the real tunnel inflow.However, the coefficient of determination R 2 P is 0.9555 by the trial-and-error method in the Reference [5], while the R 2 10 is 0.9780 by DPC in this paper (k = 10, all sectors).Because R 2 P < R 2 10 , the fitting effect of DPC is a little better.Figure 4b displays the scatter plots of the measured and calculated inflow.DPC can quickly optimize the parameters by computers, and this is a great advantage over the trial-and-error method.Moreover, DPC is stable; no matter what initial values are assigned, the calibration will converge to the same optimized values, indicating that the initial parameter assignment has little effect on the optimization.
In reality, a tunnel is excavated gradually, or the excavation is a dynamic process, so as the excavation progresses forward, the influence of the parameters of the preceding excavated sectors on the inflow calculation for subsequent unexcavated sectors will weaken.Thus, an appropriate k needs to be carefully considered.Here, we set k = 1, k = 2, and k = 3 to calibrate K and S, respectively, instead of the aforementioned k = 10. Figure 4c presents the time-dependent inflow rate for the three cases, and Figure 4d shows the related scatter plots of the measured and calculated values.
The coefficient of determination R 2 are 0.9569, 0.9604, and 0.9591 for k = 1, 2, and 3, respectively.Therefore, the result R 2 P < R 2 1 < R 2 3 < R 2 2 indicates that the fitting effect of DPCs are still slightly better than the trial-and-error method.It also demonstrates that k affects the DPC results; different k will lead to different parameters optimization and different fitting effects as well.Generally, a larger k achieves better fitting; for example, the biggest R 2 is 0.9780 when k = 10.However, a larger k also results in more computation and is more time-consuming.
The process of dynamic parameter calibration including nine stages and ten sectors is presented in Figure 5 for the case k = 3. Accordingly, the optimization results of each sector in each stage are presented in Figure 6, wherein the optimized parameters are written as Kij and Sij, where the subscript i represents the ith sector as aforementioned, while the superscript j denotes the jth times optimization for this sector.j should be a positive integer and meets 1 ≤ j ≤ k.DPC can quickly optimize the parameters by computers, and this is a great advantage over the trial-and-error method.Moreover, DPC is stable; no matter what initial values are assigned, the calibration will converge to the same optimized values, indicating that the initial parameter assignment has little effect on the optimization.
In reality, a tunnel is excavated gradually, or the excavation is a dynamic process, so as the excavation progresses forward, the influence of the parameters of the preceding excavated sectors on the inflow calculation for subsequent unexcavated sectors will weaken.Thus, an appropriate k needs to be carefully considered.Here, we set k = 1, k = 2, and k = 3 to calibrate K and S, respectively, instead of the aforementioned k = 10. Figure 4c presents the time-dependent inflow rate for the three cases, and Figure 4d shows the related scatter plots of the measured and calculated values.
The coefficient of determination R 2 are 0.9569, 0.9604, and 0.9591 for k = 1, 2, and 3, respectively.Therefore, the result indicates that the fitting effect of DPCs are still slightly better than the trial-and-error method.It also demonstrates that k affects the DPC results; different k will lead to different parameters optimization and different fitting effects as well.Generally, a larger k achieves better fitting; for example, the biggest R 2 is 0.9780 when k = 10.However, a larger k also results in more computation and is more time-consuming.
The process of dynamic parameter calibration including nine stages and ten sectors is presented in Figure 5 for the case k = 3. Accordingly, the optimization results of each sector in each stage are presented in Figure 6, wherein the optimized parameters are written as K ij and S ij , where the subscript i represents the ith sector as aforementioned, while the superscript j denotes the jth times optimization for this sector.j should be a positive integer and meets 1 ≤ j ≤ k.The optimization of the hydraulic conductivity K is relatively stable, namely, the change of K is not significant after several calibrations for the same sector.The range of the K value is 10 −3.26 -10 0.99 m/d in the sectors with relatively developed fractures, such as in sectors 2, 3, 4, 6, 8, and 9, as shown in Figure 6.On the other hand, the range is 10 −10.00 -10 −9.76 m/d in the sectors with weak fracture development, as in sections 1, 5, 7, and 10.
The value range of specific storage coefficient S is 10 −10.00 -10 −1.00 m −1 .In most sectors, the optimization of S is as stable as K.However, when sectors 2-4 and 3-5 were calibrated, S 3 of the 3rd sector changed dramatically, as shown in Figure 6c,d; this may be due to the sudden decrease of the inflow rate in the 5th sector after a steep rise in the 4th sector, seen in Figure 5c,d.Similarly, in Figure 6g,h, the S of the 7th and 8th sectors have increased obviously, in accordance with the significant inflow increase of the 9th sector, as shown in Figure 5g,h.The optimization of the hydraulic conductivity K is relatively stable, namely change of K is not significant after several calibrations for the same sector.The rang the K value is 10 −3.26 -10 0.99 m/d in the sectors with relatively developed fractures, suc in sectors 2, 3, 4, 6, 8, and 9, as shown in Figure 6.On the other hand, the range is 10 10 −9.76 m/d in the sectors with weak fracture development, as in sections 1, 5, 7, and 1 The value range of specific storage coefficient S is 10 −10.00 -10 −1.00 m −1 .In most sec the optimization of S is as stable as K.However, when sectors 2-4 and 3-5 were calibra S3 of the 3rd sector changed dramatically, as shown in Figure 6c,d; this may be due to sudden decrease of the inflow rate in the 5th sector after a steep rise in the 4th sector, in Figure 5c,d.Similarly, in Figure 6g,h, the S of the 7th and 8th sectors have incre obviously, in accordance with the significant inflow increase of the 9th sector, as sh in Figure 5g,h.
In Figure 5f-h, when the calibration process goes to sectors 5-7, 6-8, and 7-9, the coefficient of determination R 2 is a little smaller, and the fitting effect of sectors 6 and not good.Even if the fitting of the 6th sector was greatly improved compared to the sector, it was still not a satisfactory fitting.The water inflow in the 6th sector rose sha in order to make a better fit for the subsequent decline sector; the K value in the 6th se underwent a small decrease after optimization.However, even if the K and S were small in the 7th sector, its dropping trend still could not be fitted well, as seen in Fi 6e-h.In summary, the abrupt change of inflow brings great difficulty to parameter mization.In Figure 5f-h, when the calibration process goes to sectors 5-7, 6-8, and 7-9, the local coefficient of determination R 2 is a little smaller, and the fitting effect of sectors 6 and 7 is not good.Even if the fitting of the 6th sector was greatly improved compared to the 5th sector, it was still not a satisfactory fitting.The water inflow in the 6th sector rose sharply, in order to make a better fit for the subsequent decline sector; the K value in the 6th sector underwent a small decrease after optimization.However, even if the K and S were very small in the 7th sector, its dropping trend still could not be fitted well, as seen in Figure 6e-h.In summary, the abrupt change of inflow brings great difficulty to parameter optimization.

Xiema Tunnel in Chongqing, China
Xiema Tunnel is located in Beibei District, Chongqing City, Southwest China.The stratum occurring on the entrance side is steep, also the west limb of Guanyinxia Anticline, with an inclination angle of 81~87 • .The strata, which the tunnel passes through, successively from new to old include Jurassic Middle Shaximiao Formation (J 2 s) sandstone with mudstone, Xintiangou Formation (J 2 x) sandstone with mudstone, Lower Ziliujing Formation (J 1 zl) sandstone with mudstone, Lower Zhenzhuchong Formation (J 1 z) sandstone with mudstone, Triassic Upper Xujiahe Formation (T 3 xj) sandstone shale, Middle Leikoupo Formation (T 2 l) dolomite limestone, Lower Jialingjiang Formation (T 1 j) limestone dolomite, and Feixianguan Formation (T 1 f ) limestone with mudstone.The geological cross-section along the tunnel is shown in Figure 7a, and Figure 7b presents the curve of measured inflow versus location.According to the above analysis of the geological conditions of Xiema Tunnel and the geological cross-section, the rock formation in the entrance section of Xiema Tunnel is nearly vertical, which satisfies the prerequisite of using the Perrochet Equation to calculate the water inflow.
koupo Formation (T2l) dolomite limestone, Lower Jialingjiang Formation (T1j) limest dolomite, and Feixianguan Formation (T1f) limestone with mudstone.The geolog cross-section along the tunnel is shown in Figure 7a, and Figure 7b presents the curv measured inflow versus location.According to the above analysis of the geological c ditions of Xiema Tunnel and the geological cross-section, the rock formation in the trance section of Xiema Tunnel is nearly vertical, which satisfies the prerequisite of us the Perrochet Equation to calculate the water inflow.Perrochet's analytical equation of transient inflow was used to calibrate the hydra conductivity K and specific storage coefficient S of the tunnel rock mass by the trial-a error method.The results reflect the variation of the time-dependent cumulative inflow the tunnel [26]; the coefficient of determination R 2 X is as high as 0.9284.Here, we calibra K and S by using the DPC method; with the same tunnel division scheme, the whol sectors were optimized simultaneously, and R 2 20 is 0.9910 larger than R 2 X.The optim tion results are shown in Figure 7c.
The results show that K and S changed significantly when the tunnel passed thro the fault impact zone, or near the boundary of two lithologic strata, indicating the cha of permeability and storability of the strata accordingly, as shown in sectors 6-7, 10 and 17-18 in Figure 7c.Perrochet's analytical equation of transient inflow was to calibrate the hydraulic conductivity K and specific storage coefficient S of the tunnel rock mass by the trial-anderror method.The results reflect the variation of the time-dependent cumulative inflow of the tunnel [26]; the coefficient of determination R 2 X is as high as 0.9284.Here, we calibrated K and S by using the DPC method; with the same tunnel division scheme, the whole 20 sectors were optimized simultaneously, and R 2 20 is 0.9910 larger than R 2 X .The optimization results are shown in Figure 7c.
The results show that K and S changed significantly when the tunnel passed through the fault impact zone, or near the boundary of two lithologic strata, indicating the change of permeability and storability of the strata accordingly, as shown in sectors 6-7, 10-11 and 17-18 in Figure 7c.
In particular, because the rock stratum in the fault zone and its influence zone is more fractured, the permeability of the rock stratum is generally greater than that of the adjacent rock stratum.In addition, when the tunnel is excavated to the fault zone and its influence zone, the fault zone and its influence zone may have a great influence on the water ingress in the tunnel.This can be seen from the measured water inflow of the tunnel.When the tunnel is excavated to the Baimiaozi fault (the 18th sector), the water inflow increases slightly, as shown in Figure 7b.According to the inversion results, the permeability coefficient of the 18th sector is significantly higher than that of its adjacent sectors 16-17 and 19-20, as shown in Figure 7c.In conclusion, the inversion results of the proposed optimization scheme in the fault zone and its influence zone (occurrence near vertical) are reasonable.
Excluding sectors where the values of K and S largely varied, there are three types of rock mass that can be classified.Type 1: the range of K is 10 −4.08 -10 −2.50 m/d, and S is 10 −4.52 -10 −1.00 m −1 in sectors where sandstone is interbedded with mudstone and limestone with mudstone; Type 2: in sectors where sandstone and shale mainly occur, K is 10 −5.16 -10 −4.11 m/d, and S is 10 −1.31 -10 −1.00 m −1 ; Type 3: in the carbonate rock strata, K is 10 −8.99 -10 −1.20 m/d, and S is 10 −9.31 -10 −1.00 m −1 .Thus, the overall permeability of the rock mass through the tunnel is sequenced as: Type 2 < Type 1 < Type 3; the storability is: Type 1 < Type 2 < Type 3.
In order to evaluate the sensitivities of parameters K and S, we calibrated these two parameters separately.First, parameter K is calibrated, while S is fixed, keeping unchanged.Second, inversely, we calibrate S while K is fixed.
Based on the results of dynamic calibration when k = 3, the probability density function (PDF) and the cumulative distribution function (CDF) are calculated and plotted, respectively.The probability density of both the lgK value and the lgS are approximately bimodal, as seen in Figure 8.When parameter K is fixed, lgS is assigned as −1.0 m −1 in sectors with relative high storability, while lgS is −10.0 m −1 in lower sectors, and lgS is −2.1 m −1 for others.Likewise, when S is fixed, lgK is −1.0 m/d in sectors with greater permeability, while lgK is −10.0 m/d in sectors with lower permeability, and lgK is −3.4 m/d for the rest.

sonable.
Excluding sectors where the values of K and S largely varied, there are three typ rock mass that can be classified.Type 1: the range of K is 10 −4.08 -10 −2.50 m/d, and S is 10 10 −1.00 m −1 in sectors where sandstone is interbedded with mudstone and limestone mudstone; Type 2: in sectors where sandstone and shale mainly occur, K is 10 −5. 16m/d, and S is 10 −1.31 -10 −1.00 m −1 ; Type 3: in the carbonate rock strata, K is 10 −8.99 -10 −1.20 and S is 10 −9.31 -10 −1.00 m −1 .Thus, the overall permeability of the rock mass through the nel is sequenced as: Type 2 < Type 1 < Type 3; the storability is: Type 1 < Type 2 < Ty In order to evaluate the sensitivities of parameters K and S, we calibrated these parameters separately.First, parameter K is calibrated, while S is fixed, keeping changed.Second, inversely, we calibrate S while K is fixed.
Based on the results of dynamic calibration when k = 3, the probability density tion (PDF) and the cumulative distribution function (CDF) are calculated and plotte spectively.The probability density of both the lgK value and the lgS are approxim bimodal, as seen in Figure 8.When parameter K is fixed, lgS is assigned as −1.0 m sectors with relative high storability, while lgS is −10.0 m −1 in lower sectors, and lgS is m −1 for others.Likewise, when S is fixed, lgK is −1.0 m/d in sectors with greater per bility, while lgK is −10.0 m/d in sectors with lower permeability, and lgK is −3.4 m/ the rest.According to the comparison between these two calibration processes, the fittin fect changed greatly after the 13th sector.Sectors 13-17 were chosen to demonstrate details; the calculated time-varying water inflow curves are shown in Figure 9 and pared with the measured values.
The coefficient of determination R 2 K by the first way is greater than R 2 S by the se way, indicating its better fitting effect over the second way, especially in the sectors w water inflow drops sharply as seen in Figure 9.Moreover, the fitting effect of the se way is significantly reduced from sectors 14 to 16, as shown in Figure 9b,c.In Figur when calibrating sectors 15-17, R 2 S is rather small by the second way, indicating a fitting.According to the comparison between these two calibration processes, the fitting effect changed greatly after the 13th sector.Sectors 13-17 were chosen to demonstrate some details; the calculated time-varying water inflow curves are shown in Figure 9 and compared with the measured values.
The coefficient of determination R 2 K by the first way is greater than R 2 S by the second way, indicating its better fitting effect over the second way, especially in the sectors where water inflow drops sharply as seen in Figure 9.Moreover, the fitting effect of the second way is significantly reduced from sectors 14 to 16, as shown in Figure 9b,c.In Figure 9c, when calibrating sectors 15-17, R 2 S is rather small by the second way, indicating a poor fitting.In summary, the fitting effect of the one-parameter calibration by the first way is similar to the two-parameter methods.However, the one-parameter calibration by the second way is greatly different, especially in sectors where water inflow reduces significantly.Therefore, we concluded that the hydraulic conductivity K is more sensitive compared to the specific storage coefficient S, and K has a greater impact on the dynamic calibration.
Although the proposed dynamic calibration methods can optimize parameters effectively, it is worth noting some limitations.Limitation (1): We regarded the water level drawdown s in the Equation (1) as a known condition; however, the groundwater level monitoring data are actually quite lacking in vicinity of a tunnel engineering site.Thus, it may be a key parameter requiring calibration as important as K and S. Limitation (2): When using the optimization scheme, we can set the range of the optimization result first to ensure the rationality of the optimization result, we call this a parameter constraint.The parameter constraint may affect the optimization and fitting, whose effect were not investigated in this paper.These are the problems that need to be further studied and solved in the future.In addition, aside from parameter calibration and estimation, we will further develop some predictive capabilities based on this approach.First, we can evaluate the parameters uncertainty for excavated sectors, and water inflow prediction with parameter uncertainty could be obtained for unexcavated sectors.In summary, the fitting effect of the one-parameter calibration by the first way is similar to the two-parameter methods.However, the one-parameter calibration by the second way is greatly different, especially in sectors where water inflow reduces significantly.Therefore, we concluded that the hydraulic conductivity K is more sensitive compared to the specific storage coefficient S, and K has a greater impact on the dynamic calibration.
Although the proposed dynamic calibration methods can optimize parameters effectively, it is worth noting some limitations.Limitation (1): We regarded the water level drawdown s in the Equation (1) as a known condition; however, the groundwater level monitoring data are actually quite lacking in vicinity of a tunnel engineering site.Thus, it may be a key parameter requiring calibration as important as K and S. Limitation (2): When using the optimization scheme, we can set the range of the optimization result first to ensure the rationality of the optimization result, we call this a parameter constraint.The parameter constraint may affect the optimization and fitting, whose effect were not investigated in this paper.These are the problems that need to be further studied and solved in the future.In addition, aside from parameter calibration and estimation, we will further develop some predictive capabilities based on this approach.First, we can evaluate the parameters uncertainty for excavated sectors, and water inflow prediction with parameter uncertainty could be obtained for unexcavated sectors.

Conclusions
In this paper, Dynamic Parameter Calibration (DPC) is formed by combining the Perrochet Equation and the Trust Region Reflection (TRR) algorithm, which can quickly obtain the hydrogeological parameters of the constructed sector of the tunnel through the

Conclusions
In this paper, Dynamic Parameter Calibration (DPC) is formed by combining the Perrochet Equation and the Trust Region Reflection (TRR) algorithm, which can quickly obtain the hydrogeological parameters of the constructed sector of the tunnel through the measured water inflow during excavation.The study of two cases proves that this scheme can be used to obtain hydrogeological parameters of tunnels in a subvertical multilayered heterogeneous stratum.The main conclusions are as follows: 1.
The proposed Dynamic Parameter Calibration (DPC) scheme, which can automatically and quickly optimize hydrogeological parameters by computer, has great advantages over the empirical trial-and-error method.The coefficient of determination R 2 by DPC is always greater than that by trial-and-error, which means a better fitting between observation and calculation.

2.
The number of sectors being calibrated in each step can be set arbitrarily in this DPC scheme.This number has a great influence on the fitting effect, specifically, the more sectors being calibrated, the better the fitting.However, more sectors to be calibrated simultaneously means more parameters to be optimized, and this will lead to massive computation and time consumption.Moreover, the scenario with too many sectors to be optimized simultaneously deviates from the reality of dynamic construction and continuous advancement of tunnel engineering.

3.
Considering the uncertainty of parameters, the proposed calibration scheme can also conduct multiple optimizations in the same sector.The first case study of the Modane exploration tunnel shows that the hydraulic conductivity K of the same sector is relatively stable after multiple optimizations, and the change is insignificant.The optimization results of the specific storage coefficient S are relatively stable in most sectors, except in the sectors where the inflow changes dramatically.4.
For the case where the tunnel passes through heterogeneous strata with different lithologies, the proposed DPC scheme is also applicable.The example of Xiema tunnel demonstrates that K and S will change abruptly near the rock boundary and fault influence zone, reflecting the changes of permeability and storability.The parameter sensitivity analysis shows that K is more sensitive to water inflow than S.

Figure 1 .
Figure 1.Three-dimensional (3D) schematic of a tunnel progressively drilled through a subvertical multilayered system: different shades of gray represent different permeability and storability strata (modified after Perrochet et al. [5]).

Figure 1 .
Figure 1.Three-dimensional (3D) schematic of a tunnel progressively drilled through a subvertical multilayered system: different shades of gray represent different permeability and storability strata (modified after Perrochet et al. [5]).

Figure 2 .
Figure 2. Flowchart of the dynamic parameter calibration process.

Figure 3 .
Figure 3. Three-dimensional schematic of the dynamic calibration process when k = 3. Red sectors indicate the parameters being calibrated; Blue sectors indicate that the parameters are fixed or have been calibrated.

Figure 2 .
Figure 2. Flowchart of the dynamic parameter calibration process.

Figure 2 .
Figure 2. Flowchart of the dynamic parameter calibration process.

Figure 3 .
Figure 3. Three-dimensional schematic of the dynamic calibration process when k = 3. Red sectors indicate the parameters being calibrated; Blue sectors indicate that the parameters are fixed or have been calibrated.

Figure 3 .
Figure 3. Three-dimensional schematic of the dynamic calibration process when k = 3. Red sectors indicate the parameters being calibrated; Blue sectors indicate that the parameters are fixed or have been calibrated.

Water 2023 ,
15, 2702 6 of 14 the area N, as shown in Equations (

Figure 4 .
Figure 4. Comparison between the calculated and observed transient groundwater inflow of the Modane tunnel: (a,c) time series curves; (b,d) scatter diagrams.(Water inflow measurements are from Perrochet et al. [5]).

Figure 4 .
Figure 4. Comparison between the calculated and observed transient groundwater inflow of the Modane tunnel: (a,c) time series curves; (b,d) scatter diagrams.(Water inflow measurements are from Perrochet et al. [5]).

Water 2023 , 8 Figure 5 .
Figure 5.Comparison between calculated versus measured inflow rate in the process of dyn parameter calibration (k = 3, R 2 is the coefficient of determination of the optimized sectors).

Figure 5 .
Figure 5.Comparison between calculated versus measured inflow rate in the process of dynamic parameter calibration (k = 3, R 2 is the coefficient of determination of the optimized sectors).

Figure 6 .
Figure 6.Optimized parameters of each sector in all 9 steps by DCP (k = 3).

Figure 7 .
Figure 7. (a) Transverse geological cross-section of Xiema tunnel; (b) measured groundwater in into tunnel during excavation: the dark gray and white color of the tunnel entrance represents length of the different sectors, totaling 20 sectors; and (c) optimized values of K and S for each se (Plots (a,b) were modified after Long [26]).

Figure 7 .
Figure 7. (a) Transverse geological cross-section of Xiema tunnel; (b) measured groundwater inflow into tunnel during excavation: the dark gray and white color of the tunnel entrance represents the length of the different sectors, totaling 20 sectors; and (c) optimized values of K and S for each sector (Plots (a,b) were modified after Long [26]).

Figure 9 .
Figure 9.The comparison between the analytical calculation of the water inflow by DPC against the measured values.(R 2 K and R 2 S are the determination coefficients of two calibration processes, respectively.The shadow on the left side of the dashed line indicates the sectors with a fixed parameter, whereas the right side are the sectors with a parameter to be calibrated).

Figure 9 .
Figure 9.The comparison between the analytical calculation of the water inflow by DPC against measured values.(R 2 K and R 2 S are the determination coefficients of two calibration processes, respectively.The shadow on the left side of the dashed line indicates the sectors with a fixed parameter, whereas the right side are the sectors with a parameter to be calibrated).