Numerical Simulation of Nanosilica Sol Grouting for Deep Tunnels Based on the Multifield Coupling Mechanism

Grouting is an e ﬀ ective technical way for the construction of deep tunnels in unfavorable geological conditions. The ﬂ uid-solid-chemical coupling mechanism of grouting process is analyzed from the following three aspects: in ﬂ uence of physical properties of silica sol on permeability coe ﬃ cient, dynamic changes of porosity and permeability of geotechnical media with seepage pressure, and governing equations for ﬂ ow and mass transfer characteristics. A dynamically changing model for nanosilica sol grouting in deep tunnels is established, considering the changing physical properties of grout and surrounding rock. Based on the Xianglushan Tunnel of Yunnan Water Diversion Project, the temporal and spatial evolution of silica sol grout is studied. The e ﬀ ect characteristics of grouting pressure and initial permeability are clari ﬁ ed. The rationality of this model is veri ﬁ ed by classical Newtonian ﬂ uid grouting theory. The main conclusions: with the molar concentration as the index, the grout range can be divided into the raw grout region and the transition region; with the decrease of the grouting pressure, the growth rate of the normal grouting radius and the axial grouting radius will gradually decrease; due to the mechanical dispersion and molecular di ﬀ usion, the range of the transition region will gradually increase with time. The ratio of the transition region to grouting radius ﬂ uctuates slightly with time under the initial permeability of 5 D. The ﬂ uctuation increases with the decrease of initial permeability, and the average ratio increases with the decrease of grouting pressure. This study can provide theoretical guidance for grouting design of deep tunnel engineering.


Introduction
Due to the influence of crustal movement, human excavation engineering, and weathering erosion, the rock mass contains not only macroscopic defects such as joints and cracks but also microscopic defects such as microcracks and micropores [1][2][3][4]. Grouting is an effective technical way for the construction of deep tunnels in unfavorable geological conditions [5][6][7]. Due to nanoparticle size, low viscosity, high injectability, and friendly to ecological environment, silica sol grouting has been widely used in petroleum engineering, groundwater protection engineering, tunnel construction, nuclear waste storage, liquefiable soil improvement, and other geotechnical fields [8][9][10][11].
Wang et al. [10,12] investigated the rheological properties of nanosilica sol grout and derived an equation of realtime viscosity and permeability length. Fujita and Kobayashi [13] studied the transport mechanism of colloidal silica in unsaturated sand and the effect of charging properties of sand and silica particles by sand column experiment. In order to study the penetration law of silica sol in microcracks with a width of 40-50 microns, fluorescent whitening agent-ultraviolet observation method, microscopic observation method, and hydraulic experiment method were used to jointly predict and verify the penetration length and sealing effect of silica sol [14]. The influence of gel effect on the migration ability of silica sol in liquefied soil was studied by using soil column experiments and a pilot-scale facility [15][16][17]. The migration process of silica sol along the whole length of soil column, the gel capacity of silica sol when it reached the target area, the injection amount of silica sol, and the migration amount of silica sol with the concentration of chloride ion as the indexes were monitored during the soil column experiment. The flow laws of silica sol, epoxy resin, polyurethane, and cement-based grout in onedimensional and two-dimensional fractured rock mass with different grouting modes (constant flow grouting, constant pressure grouting, and constant energy grouting) were studied by theoretical analysis [18]. An analytical model of colloidal silica gelling in groundwaters with varying pH and background electrolyte concentrations was developed and validated [19]; the model was then applied to a hypothetical case study to demonstrate its use in grout design, based on published in situ groundwater data from the Olkiluoto area of Finland. The pulsating grouting technology of silica sol was tested in a 500-meter-deep tunnel at Mizunami Underground Research Laboratory in Japan [20]. Pan et al. [21,22] studied the influence characteristics and mechanism of fragmental size of broken coal mass on the injection regularity of silica sol grouting. Injection regularity of silica sol grouting in broken coal mass generally presented a three-stage regularity. The fragmental size had a great impact on the 2nd and 3rd stages.
In summary, previous experts have achieved fruitful research results about the injection law of silica sol grouting through numerical simulations, column tests, and field experiments. However, there are still some deficiencies in the study of the synergy effect of grouting materialsurrounding rock physical properties on silica sol grouting in deep tunnels, and the study of multifield mechanism for silica sol grouting in deep tunnels is still insufficient.
In this paper, in order to explore the theoretical mechanism of nanosilica sol grouting technology for deep tunnels in unfavorable geological conditions, a dynamically changing model for silica sol grouting in deep tunnels is established, taking into account the changing physical properties of grout and surrounding rock. Based on the project background of Xianglushan Tunnel of Central Yunnan Water Diversion Project, the temporal and spatial evolution of Darcy velocity, viscosity, molar concentration, and fluid pressure of silica sol grouting are analyzed. The effects of grouting pressure and initial permeability of surrounding rock on permeability grouting characteristics of silica sol are explored. The rationality of this model is verified by combining the column theory and sphere theory.

Engineering Background and Grouting Material
Xianglushan Tunnel was located on the head of the Dali I Section, which was the first section of the main channel of the Central Yunnan Water Diversion Project (Figure 1). At the feasibility study stage, the total length of the section was 63.426 km [23]. For Xianglushan Tunnel, the burial depth in the area to the north of the trough valley of the Runan River was 600~1000 m. According to the stress assessment results of tunnel surrounding rock in the geological survey stage, the average weight of rock mass can be calculated as 26.5 kN/m 3 . The vertical stress of the surrounding rock at a buried depth of 1000 m was 26.50 MPa, which was a high level of ground stress. The statistics of conventional water pressure tests results showed that the overall permeability of the rock formations in the study area was weak water permeability (1 ≤ q ðpermeability rateÞ < 10 Lu, 10 −5 ≤ K ðpermeability coefficientÞ < 10 −4 cm/s), and some were medium permeable. The permeability of karst strata was stronger than that of nonkarst strata. Part of the medium permeability is mainly concentrated limestone, dolomitic limestone, argillaceous limestone, and sandstone of Songgui formation. The drilling and blasting method and TBM method were used to construct the Xianglushan Tunnel, in which the TBM section was 35.37 km long. According to the information of exploration and design, the tunnel diameter of downstream TBM section was 9.8 m, in which Class IV and V surrounding rocks accounted for 53.14% (Table 1), and large deformation surrounding rocks accounted for 16.62%. Under the condition of naked tunnel, the length of the main tunnel section of Xianglushan Tunnel was 33.085 km with a specific yield greater than 3 m 3 /m·d in dry season.
Considering environmental protection, ecology, and flow control in the tunnel, in order to avoid more leakages in the tunnel, reduce the impact of pumping and drainage on groundwater environment and tunnel construction, groundwater should be plugged and drained in a way that followed the principle of "first plugging, and limited drainage." When the local hydraulic head of groundwater was less than 60 m, it was better to adopt the full plugging method; when the hydraulic head of groundwater was more than 60 m, drainage and guide mode should be adopted.
Silica sol was selected as sealing grouting material. It was a one-component nanogrouting material, which belonged to colloidal solution, odorless and nontoxic. Silica sol had a typical double-electron layer structure. The dominant size of the silica particles in silica sol was 8-12 nm. When strong electrolyte solution (i.e., catalyst) was added, the doubleelectron layer of colloidal particles was compressed and thinned, and the gel reaction began to take place. The viscosity of colloidal particles increased gradually [25,26]. The catalyst chosen in the experiment was NaCl solution ( Table 2). The gel time of silica sol can be adjusted according to the amount of catalyst. The gel reaction products had a good transmittance and compact microstructure ( Figure 2).

A Dynamically Changing Model for Silica Sol
Grouting in Deep Tunnels Based on the Multifield Coupling Mechanism The injection process in a tunnel is a complex and changeable system problem: (a) geotechnical medium is multiscale, and different stress states of surrounding rocks at different positions in a tunnel lead to uneven development of fracture and pore; (b) geotechnical environment is multi-interface, such as grout-rock, grout-gas, and grout-water interfaces; 2 Geofluids the permeability grouting presents nonlinear characteristics because of the different states, structures, and physical and chemical properties of the interfaces; (c) in the grouting pro-cess, there are many reactions, such as convection, molecular diffusion, mechanical dispersion, adsorption and retention, and gel reaction between silica sol and catalyst. These reactions will change the migration characteristics of fluids [27]. In this paper, the injection process of silica sol in deep tunnels is regarded as a multireaction coupled flow in porous media. The fluid-solid-chemical coupling mechanism of grouting process is analyzed from the following three aspects.      4 Geofluids coefficient of grout is mainly restricted by two aspects: (a) the redistribution of surrounding rock stress field caused by engineering disturbances such as excavation and grouting; (b) the change of viscosity and density of grout caused by self-gel and dilution in groundwater. Self-gel is embodied in the dependence of viscosity on time; adsorption retention and dilution in groundwater are embodied in the dependence of viscosity on concentration and density on concentration.

Three-Dimensional Surface Equation for
Viscosity-Time-Concentration of Grout. For silica sol with mass fraction of 15%, the molar concentration of silica is 2750 mol/m 3 . The grout is prepared with silica sol : catalyst volume ratio of 9 : 1, and the viscosity is measured with an NDJ-5S rotary viscometer [28]. The function expression of fitting viscous time curve is as follows: μ = 7:32 × 10 −13 e t/1309:55 where t is in seconds and μ is in Pa·s.
By using the linear interpolation algorithm [29], the three-dimensional surface equation of grout with respect to time and concentration can be obtained: where the preferable range of t is 0-36000 s and the preferable range of c is 0-2750 mol/m 3 . As shown in Figure 3, a three-dimensional surface diagram is drawn.

Relationship between Density and Molar
Concentration of Grout. The molar concentration of grout (9 : 1) is 2750 mol/m 3 , and the density is 1097 kg/m 3 ; however, the molar concentration of water is 0 mol/m 3 and the density is 1000 kg/m 3 . The molar concentration of silica sol grouting is basically proportional to its density, and the density of grout can be expressed as follows: where ρ is the density of grout, kg/m 3 .

Dynamic Changes of Porosity and Permeability of
Geotechnical Media with Seepage Pressure. In the process of geotechnical grouting, the seepage pressure of grouting fluid will change the spacing or pattern of geotechnical particles, so that the actual porosity, permeability, and other physical parameters will dynamically change with the change of the seepage pressure, thereby affecting the seepage process. When neglecting the temperature effect, the relationship between the dynamic changes of porosity and permeability in geotechnical media under seepage pressure can be given as follows [30]: where φ is the dynamic porosity; k is the dynamic permeability; φ 0 is the initial porosity; k 0 is the initial permeability; ε v is the volumetric strain; Δp = p − p 0 (p 0 is the initial pressure; p is current pressure); E s is the bulk modulus of geotechnical media. Combined with Equations (2), (3), and (5), a variable permeability coefficient model of silica sol grouting in the tunnel can be constructed: where K is the permeability coefficient.

Governing Equations for Flow and Mass
Transfer Characteristics 3.3.1. Governing Equations for Grout Flow. The governing equations of grout flow can be expressed by Darcy's law [2].
where u is the Darcy velocity, m/s; ∇D is the unit vector in gravity direction; D is the vertical coordinates.

Governing Equation for
Mass Transfer. Silica sol grouting approximately satisfies transport of diluted species in porous media [11]. Assuming that only convection, molecular diffusion, and mechanical dispersion are considered and that adsorption is not considered, the basic equation of solute transport is shown in the following: where c i is the concentration of component i, mol/m 3 ; D i is the molecular diffusion, m 2 /s; R i is the reaction rate of component i, mol/(m 3 ·s).
Because the gel reaction of silica sol is the gradual connection of colloidal particles into a network structure, the gel reaction can be regarded as not affecting the molar concentration of silica.

Numerical Model and Boundary Conditions
4.1. Parameter Values. COMSOL software is used for numerical calculation. As shown in Figure 4, the model size is 100 × 100 m, the diameter of tunnel is 9.8 m, the diameter of grouting borehole is 5 cm, and the length of borehole is 3 m. The numerical simulation process can be divided into three steps: in situ stress initialization, tunnel excavation, and borehole grouting. The values of parameters in the numerical model are shown in Table 1.

Boundary Conditions for Solid Mechanics.
The upper boundary condition is given stress of 25 MPa, the lower boundary condition is fixed constraint, the left boundary condition is symmetrical boundary condition, and the right boundary condition is roller support. The boundary condition of the tunnel is simplified to free. After boring, the grouting pipe is put in, and the boundary condition is simplified to zero radial displacement.

Boundary Conditions for Darcy's Law.
The upper, lower, and right boundary conditions of the model are given hydraulic head of 50 m, and the left boundary condition is symmetrical. After tunnel excavation, timely grouting behind shield closes the cracks in surrounding rock, and the boundary condition for surrounding rock of tunnel is nonflowing. The boundary condition of borehole is given pressure of 2 MPa to simulate grouting pressure.

Boundary Conditions for Transport of Diluted Species in Porous
Media. The upper, lower, and right boundary conditions of the model are given pressure of 0 mol/m 3 , and the left boundary condition is symmetrical. The boundary condition of borehole is 2750 mol/m 3 , which simulates the grout concentration at the injection end. The viscosity, molar concentration, and permeability coefficient in the raw grout region remain basically unchanged at the same time. The viscosity and molar concentration in the transition region decrease with the distance from the grouting hole, while the permeability coefficient increases with the distance from the grouting hole. Moreover, when the gel time is not exceeded, the viscosity increases slowly and remains at a low level for a long time, which basically coincides with the indoor test.

Results and Discussion
The Darcy velocity at the same location decreases with the increase of time. The maximum Darcy velocity at the 15th minute is 3:21 × 10 −4 m/s, the maximum Darcy velocity at the 480th minute is 1:37 × 10 −4 m/s, and the maximum Darcy velocity at the 560th minute drops sharply to 1:32 × 10 −5 m/s. This is closely related to the late sharp increase in viscosity.
The normal grouting radius is higher than the axial grouting radius; the normal grouting radius, the axial grouting radius, and the transition region all increase with time, but the ratio of the transition region to the grouting radius remains at 36.22%.

Effect of Grouting Pressure on the Injection
Characteristics. Silica sol grout can be regarded as Newtonian fluid. Column theory and sphere theory for uniform permeability of Newtonian fluid in stratum are shown in Equations (9) and (10) [32]. According to the values of parameters in Tables 2 and 3, the theoretical grouting radius is calculated: where h 1 is the difference between grouting pressure head and groundwater head, cm; r 1 is the grouting radius, cm; β is the viscosity ratio of grout to water; and r 0 is the radius of grouting pipe, cm. Under the same initial permeability of 5 D, the numerical models with a grouting pressure of 2 MPa, 1.5 MPa, and 1 MPa are calculated, respectively. As shown in Figure 6, the numerical simulation and theoretical calculation of grout frontier under different grouting pressures are presented. The normal grouting radius is slightly higher than the axial grouting radius, the axial grouting radius is significantly higher than the column theoretical value, and the column theoretical value is higher than the sphere theoretical value. With the decrease of grouting pressure, the growth rates of normal grouting radius and axial grouting radius gradually decrease. When the grouting pressure is 2 MPa and 1.5 MPa, respectively, the numerical simulation value is larger than the calculation value of column theory, and the excess value is equal to the range of the transition region. The coincidence between the axial raw grout region and the column theoretical value is high, which indicates that column theory can explain the injection characteristics of raw grout region. The range of the transition region increases with time, which is attributed to the increasing effect of mechanical dispersion and molecular diffusion. The ratio of transition region to grouting radius fluctuates slightly with time under the same grouting pressure, and the average ratio increases slightly with the decrease of grouting pressure. The grouting design of the tunnel engineering should take into account not only the scope of the raw grout region but also the affluence factor of the transition region.

Effect of Initial Permeability on the Injection
Characteristics. Under the same grouting pressure of 2 MPa, the numerical models with an initial permeability of 5 D, 0.5 D, and 0.05 D are calculated, respectively. As shown in Figure 7, the numerical simulation and theoretical  The normal grouting radius is slightly higher than the axial grouting radius, and the axial grouting radius is significantly higher than the value of column theory. With the decrease of initial permeability, the growth rates of normal and axial grouting radius decrease gradually. Under all initial permeability conditions, the axial raw grout region has a high degree of coincidence with the column theory, which indicates that the column theory can explain the injection characteristics of the axial raw grout region. And under the initial permeability of 0.05 D, the range of the   10 Geofluids raw grout region is basically equal to the calculated values of column theory and sphere theory, which indicates that in the geotechnical medium media with low permeability, the multifield coupling mechanism of permeability grouting can explain the column theory and sphere theory well. The range of the transition region increases with time. Under the initial permeability of 5 D, the fluctuation of the ratio of the transition region to grouting radius is small, and the fluctuation increases with the decrease of the initial permeability, and the average ratio increases gradually with the decrease of  [17] Note: * Stocks-Einstein equation can be used to calculate molecular diffusion [31].   11 Geofluids the grouting pressure. The results show that for the grouting of low-permeability rock and soil medium, although the grouting radius is small, the range of the concentration gradient increases. Therefore, it is necessary to make full use of this phenomenon to develop a grouting material which is easy to diffuse and disperse. In addition, the maximum grouting radius is more than 16 m, so the "super grouting" technology can be developed without affecting the construction period. The fully enclosed grouting in full section of tunnel can be completed by using only two symmetrical grouting boreholes.

Conclusions
In order to explore the theoretical mechanism of nanosilica sol grouting technology for deep tunnels in unfavorable geological conditions, a dynamically changing model for nanosilica sol grouting in deep tunnels is established, taking into account the changing physical properties of grout and surrounding rock. Injection characteristics of nanosilica sol grouting in deep tunnels affected by grouting pressure and initial permeability are analyzed based on the project background of Xianglushan Tunnel of Central Yunnan Water Diversion Project. The rationality of this model can be verified by combining the column theory and sphere theory for Newtonian fluid. The main conclusions are as follows: (1) The grout range can be divided into the raw grout region and the transition region with the molar concentration of silica as the index. Viscosity, molar concentration, and permeability coefficient in the raw grout region remain basically unchanged at the same time; in the transition region, viscosity and molar concentration decrease with the distance from the grouting hole, while the permeability coefficient increases with the distance from the grouting hole (2) With the decrease of grouting pressure, the growth rates of normal grouting radius and axial grouting radius will gradually decrease; due to mechanical dispersion and molecular diffusion, the range of the transition region will gradually increase with time (3) The ratio of the transition region to grouting radius fluctuates slightly with time under the initial permeability of 5 D, and the fluctuation increases with the decrease of initial permeability, and the average ratio increases with the decrease of grouting pressure

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.