Back Analysis of Geomechanical Parameters of Rock Masses Based on Seepage-Stress Coupled Analysis

Abstract

Given that rock masses are complex, the geomechanical parameters of rock masses are hard to determine in underground engineering. In this paper, the inverse model and method are established to identify the parameters based on the coupled stress and fluid flow model combined with the finite element method and adaptive genetic algorithm. Moreover, the model and method are applied in the Lianghekou highway tunnel, and the initial permeability coefficients of the stratum and the lateral pressure coefficients of the initial ground stress are identified by the back analysis with relative errors of the measured and fitted values at measuring points below 5%. Results show that the inverse model and method are effective and sound.

Xianghui Deng 1 and Dongyang Yuan 1 and Dongsheng Yang 1 and Changsheng Zhang 1

Academic Editor:Roberto Fedele

School of Civil and Architecture Engineering, Xi'an Technological University, Xi'an 710032, China

Received 8 October 2016; Revised 13 December 2016; Accepted 23 January 2017; 19 February 2017

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Rock masses exist in certain geological environments. Groundwater and ground stress are the most important factors that influence the geological environment. Water flow in rock mass changes the initial ground stress status of rock masses, and the change of their stress status influences the characteristics of the fluid flow in the rock masses. The fluid flow and stress affect each other and cause coupling, regardless of whether the fluid flow changes first or vice versa [1].

The phenomenon of coupled stress and fluid flow has been widely studied in academic and engineering circles, and scholars have conducted extensive research on the seepage-stress coupled model. In 1969, Snow first proposed the mathematical equations that show the effect of normal stress on the permeability coefficient by conducting a single-fracture permeability test [2]. In 1974, after conducting a pumping test and theory analyses, Louis proposed that the seepage discharge of fractured rock masses decreases with the increase of the normal stress, and he proposed a corresponding empirical formula [3]. In 1975, Nelson presented the empirical formula of permeability coefficients based on the Navajo sandstone sample seepage-stress test [4]. Nelson's empirical function first displayed the influence of effective stress on the rock mass permeability. Through laboratory tests, Y. Z. Zhang and J. C. Zhang [5] concluded that if the initial fracture width was small, then the seepage discharge and the stress of the fractured rock mass would not have a negative exponential function relationship with each other but a biquadrate relationship. They also suggested that the seepage discharge of the fractured rock mass declined with the increase of the compressive stress and that the seepage discharge increased with the unidirectional compressive stress parallel to the fracture surface. By the investigation and theoretical analysis, Su and coauthors proposed a negative correlation between fracture permeability and fracture normal effective stress [6]. Through the single-fracture 3D stress test, Chang et al. discovered that the fracture fluid flow was influenced by the normal stress and was strongly affected by the fracture lateral deformation caused by the fracture lateral stress [7]. Furthermore, the negative exponential formula was established between the fluid flow and the normal stress. Many researchers have also realized similar relationships in subsequent research [8, 9]. In general, seepage-stress coupling is reflected in the following aspects: on the one hand, effective stress (normal and lateral) controls the fracture width and other geometric shapes of rock masses or the porosity of porous mediums. The effective stress determines the characteristics of fluid flow, indicating the effects of the stress field on the seepage field [10, 11]. On the other hand, the groundwater influences or changes the rock mass structure by imposing hydrostatic and hydrodynamic pressure, thereby changing the stress status of the rock mass that shows the effect of seepage field on the stress field. The fluid flow and stress interaction enable the rock mass to maintain a dynamic balance [12].

The analysis of the hydromechanical behavior of rock masses remains an important topic in rock mechanics. It is a critical phenomenon in ongoing challenging issues such as tunneling under high groundwater pressures and the extraction of hydrocarbons from deep and pressurized petroleum reservoirs. Despite continuing and extensive efforts, such analysis continues to be difficult. First, establishing a reasonable hydromechanical constitutive model is difficult. The major difficulty in modeling the fluid flow in fractured rocks involves handling the solid-fluid interaction. The equivalent continuum approach and the discrete fracture network approach have been developed based on the mechanical and hydraulic natures of the rock mass for such modeling [13, 14]. The continuum and discrete approaches have been combined to propose the dual permeability model [15]. In this model, flow in natural pores and cracks are governed by different equations, which may or may not be coupled [16]. Many hydromechanical models have been proposed to overcome these problems. For example, Shao et al. established a coupled constitutive model for anisotropic damage and permeability variation in brittle rocks under deviatoric compressive stresses. The formulation of the model is based on experimental evidences, and the main physical mechanisms involved in the scale of microcracks are considered [17]. From a phenomenological point of view, on the one hand, a microscopic approach is often used to analyze the permeability evolution by the fluid flow through cracks [18]. On the other hand, a macroscopic approach is appropriate for studying the mechanical characteristics of materials, such as the stress-strain relationship after damage. Therefore, Pereira and Arson established the double-porosity model based on a relationship between the microscopic and macroscopic damage tensors, which can simulate the flow through the network of cracks/porosity and evaluate the equivalent permeability [19]. Moreover, Pereira and Arson believed that the pore size distribution (PSD) of the material is coupled to the mechanical behavior of the rock and models the influence of deformation and damage on the permeability and retention properties of cracked porous media [20, 21]. De Bellis et al. referred to the above research to simplify the damaging porous material model through consistent linearization [22]. Their model consists of nested families of equispaced frictional and cohesive faults in an otherwise elastic matrix material. The linear kinematic model preserves the main microstructural features of the finite kinematic one but offers superior computational performance. Meanwhile, determining the parameters of the hydromechanical model is very difficult for actual engineering projects. Three methods are commonly used, that is, theoretical analysis, field measurement, and the back analysis method [23, 24]. First, assumptions are generally made for the derivation, which vary greatly from the actual condition, causing difficulties in applying the derived formula. Test methods consist of laboratory and in situ tests. Laboratory tests usually have an obvious "size effect," and the accuracy of parameters cannot be guaranteed. Meanwhile, in situ tests have a limited measurement range, and the measurement result only indicates the characteristics of rock masses near the sampling point. Several uncertainties occur in identifying the parameters caused by the disturbance of sampling. In addition, in situ tests have disadvantages such as data divergence, less representation, and high costs. Comparatively, the back analysis method is based on the measured physical information (the displacement, strain, water level, and other factors) which reflects the systematic mechanical behaviors. The inversion model can be used to obtain initial parameters of the surrounding rock, and even the inverse model is analyzed occasionally. Now, the back analysis method is widely used because it is a relatively easy and cost-effective technique.

According to the mechanical behavior, the seepage-stress coupled model can reflect the real stress and fluid flow characteristics of rock masses. At present, the back analysis mostly concentrates on the inversion of the model parameters in the single field (the single field is the uncoupled problem, such as the seepage field, stress field, and temperature field), considering that the coupled interaction is obviously different from that of the single field. If the result of a single field back analysis is directly applied to the coupled stress and fluid flow model, then large errors may be generated. Therefore, identifying the parameters with the coupled model is necessary.

Until now, identifying parameters has been the problem when the seepage-stress coupling is considered. Three research ideas have been proposed, one of which is the parameter inverse of the equivalent continuous model for coupled stress and fluid flow analysis, that is, conducting a back analysis on the parameters of the model based on various types of monitoring data and coupling the forward analysis method in a continuous medium or a region considered as a continuous medium [25]. The credibility of the inversion result can be improved significantly compared with that in a single field. The second idea involves considering the rock mass as a discrete fracture medium when the back analysis is conducted. The fractured rock mass or a certain region is considered as a discrete medium, and the discrete fracture model is established in the fluid flow analysis. Meanwhile, the field monitoring data are employed for fitting, and the minimum function is taken as the objective function to identify the parameters of the coupled model [26]. Given the distribution complexity and randomness of the fractures in the rock mass, conducting numerical simulations for the discrete fracture distribution in actual engineering is difficult [27, 28]. Therefore, this idea is rarely applied in actual engineering projects. The third idea is the dual permeability coupled model [29, 30]. The actual monitoring information is applied to identify parameters, which is similar to the process of the two aforementioned ideas. According to this idea, the continuous medium model is adopted when the region has less fracture and poor permeability. The discrete fracture model is adopted when the region has a large fracture or good permeability, which largely depends on the fractures. The third idea is theoretically reasonable. However, the numerical model is complex, and numerous factors have to be considered. Therefore, applying it in practical engineering is complicated.

In this paper, the inverse model is established based on the equivalent continuum coupled stress and fluid flow model and is combined with the finite element method (FEM) and adaptive genetic algorithm (AGA). Meanwhile, the optimal parameter combination of the coupled model is identified.

2. Equivalent Continuum Model for Coupled Stress and Fluid Flow

2.1. Equivalent Continuum Mathematical Model for Coupled Stress and Fluid Flow

Generally, the water quantity is constant when water flows in the rock mass. According to the water balance principle, the mathematical equation of water balance can be established as follows: [figure omitted; refer to PDF] where _vx , _vy , and _vz are the conductivity velocity of rock mass in X, Y, and Z directions; ρ is the water density; α is the compression coefficient of rock mass; α is the compression coefficient of water; n is the porosity ratio of rock mass; h is the hydraulic head; t is time. The formula of Darcy's law is shown as follows: [figure omitted; refer to PDF]

Combining (1) with (2), (3) is obtained as follows: [figure omitted; refer to PDF] where _Ss is the unit storage volume; _kx , _ky , and _kz are the hydraulic conductivity of rock mass in X, Y, and Z directions; Q is the source sink term; _(h)Γ1 is the hydraulic head boundary condition; and _{kn((∂h/∂n))Γ2} is the discharge boundary condition.

According to the properties and seepage characteristic of the fractured rock mass, the seepage-stress coupled model can be classified into the equivalent continuous coupled model, the discrete fracture coupled model, and the dual permeability coupled model [31]. The aforementioned analysis shows that the discrete fracture coupled model and dual permeability coupled model are difficult to apply in the actual engineering. The representative element volume (REV) is relatively small for some rock masses. Therefore, the equivalent continuous coupled model is used for the seepage analysis.

The seepage control equation can be deduced from formula (1). The seepage control in Equation [32] is as follows: [figure omitted; refer to PDF] where P is the seepage water pressure, including hydrostatic pressure and hydrodynamic pressure; [K] is the total seepage matrix; {Q} is the source sink term; [S] is the unsteady seepage storage matrix; and (∂P/∂t) is the array that shows the change rate of the hydraulic head with time.

In an equivalent continuous medium, the FEM discretized process of the stress field is derived in [32]. The stress control equation is shown as follows: [figure omitted; refer to PDF] where (σ) is the total stress array; (F) is the known load array; {U} is the displacement array; (_Kn ) is the stiffness matrix; [D] is the elastic matrix; and (B) is the geometric matrix.

The linear interpolation functions are used when the FEM method is adapted to solve coupled problem. By combining the seepage equations of equivalent continuous medium, stress field equations, and seepage-stress coupling formula, the mathematical model for the coupled stress and fluid flow analysis [33] is as follows: [figure omitted; refer to PDF] where K=f(_σe ,P) is the empirical formula of the coupled stress and fluid flow, _σe is the effective normal stress, and K is the permeability coefficient tensor. The permeability coefficient tensor of the equivalent continuum coupled model was listed in [34]. Details are as follows: [figure omitted; refer to PDF] where _Kij is the component of the permeability coefficient tensor; _Kx , _Ky , and _Kz are the coefficient tensors of the X, Y, and Z directions. The permeability coefficients _Kx , _Ky , and _Kz are used in seepage analysis of the equivalent continuous model. The permeability coefficient K is closely related to the effective normal stress and the seepage pressure according to the current research. In fact, some scholars think the permeability coefficient is positively correlated with the water pressure but negatively related to the stress [30, 34]. The initial permeability coefficients of the X, Y, and Z directions (_K0x , _K0y , and _K0z ) will be identified in this paper.

2.2. Numerical Method of Equivalent Continuum Coupled Model

In the equivalent continuum coupled model, the iteration method is used for the analysis. Before the analysis process, the numerical analysis programs of the seepage field and stress field are set up, respectively. The simulation result of the single field is employed as the boundary conditions of each other and the external load. When the calculation precision meets the convergence condition, the simulation result is obtained using an iterative calculation. The basic steps are as follows:

(1) The initial stress field of rock masses _{^σ(0) ij} is calculated.

(2) The initial permeability tensor _{^K(0) ij} of different strata in the calculation area is computed according to the calculated initial stress field _{^σ(0) ij} .

(3) The seepage field is estimated according to the boundary conditions and the initial permeability tensor.

(4) The hydrodynamic pressure and hydrostatic pressure are determined according to the seepage field calculation results. Then other load increments are considered, the stress increment Δ_σij is analyzed, and the stress field _{^σ(1) ij} of this moment is obtained.

(5) The seepage tensor _{^K(1) ij} is assessed in the new stress field according to the stress field and the empirical relationship of the two fields.

(6) Steps (3) to (5) are repeated until the seepage field and stress field calculation errors in the adjacent time satisfy the convergence precision.

3. Parameter Optimization of the Inversion Method

3.1. Main Parameters to Be Inversed

Many parameters are required when conducting the coupled analysis, which include the parameters of the seepage field, the parameters of the stress field, and the coupled coefficients of the two fields. In this process, parameters that need to be determined are classified into four types according to their properties: (1) physical and mechanical parameters of the stress field, namely, the rock mass gravity γ, elasticity modulus E, Poisson ratio μ, cohesion C, and internal friction angle [straight phi]; (2) parameters of the seepage field, including the permeability coefficient of different strata or the facture aperture, and the normal stiffness of the facture rock masses; (3) parameters related to the in situ stress field, for example, the lateral pressure coefficient of initial ground stress; and (4) the coefficients of the coupled empirical relationship. The details are shown in Table 1.

Table 1: Parameters needed in the equivalent continuous coupled model.

Parameters of the seepage field	Parameters of the stress field	Coupled model parameters	Physicomechanical parameters of the stratum
Permeability coefficients of different strata Kx , Ky , and Kz	Displacements of structures or surrounding rocks ui , stresses of structures or surrounding rocks σi , and the lateral pressure coefficients of ground stresses λx and λy	Coefficients of the coupled empirical relationship	Gravity γ	Modulus of elasticity E	Poisson ratio μ	Cohesion C	Internal friction angle [straight phi]

Several parameters are determined easily by the laboratory tests, such as the rock mass gravity γ, elasticity modulus E, Poisson ratio μ, cohesion C, and internal friction angle [straight phi]. Other parameters (displacements of structures or surrounding rocks ^ui , stresses of structures or surrounding rocks ^σi ) are also identified by the field tests. In fact, determining the permeability coefficients of different strata, the coefficients of the coupled empirical relationship, and the lateral pressure coefficients of ground stress by the tests is difficult. Therefore, these parameters are usually selected as the inversed parameter.

3.2. Objective Function of Back Analysis

Considering that the hydraulic heads and displacement information belong to different series, in the calculation area, the minimum error F between the measured value and computation value in the measuring points is taken as the objective function. Therefore, the objective function of back analysis is constructed as follows: [figure omitted; refer to PDF] where _hi is the calculated hydraulic head value at the measuring point i; ^{hi[variant prime]} is the measured hydraulic head value at the measuring point i; N and T are the numbers of hydraulic head and displacement measuring points, respectively; _ui is the calculated displacement value at the measuring point i; and ^{ui[variant prime]} is the measured displacement value at the measuring point i. In the aforementioned formula, the relative values of the hydraulic head and displacement are adopted to make the objective function a dimensionless numerical function. Thus, the dimensionless numerical function can avoid other problems caused by dimensional differences and facilitate the algorithm convergence judgment during the optimization.

3.3. Adaptive Genetic Algorithms

In the traditional genetic algorithm, crossover rate _Pc and mutation rate _Pm are the key factors influencing the behavior and performance of genetic algorithm. However, these factors are predetermined in the traditional genetic algorithm, and _Pc and _Pm should be determined through repeated tests specific to different optimization problems so that they can be trapped by local optimal solution. To solve the problem, an AGA, in which _Pc and _Pm change automatically with the adaptability, is presented in this paper [35]. In an AGA, the algorithm increases the values of crossover rate and mutation rate automatically according to the individual fitness. When the individual fitness is higher than the average fitness corresponding to lower _Pc and _Pm , the next generation can be protected by the solution. Instead, when the individual fitness is lower than average, the AGA adopts the higher crossover rate and mutation rate. Thus, the optimal solution obtained can prevent the algorithm from falling into the local optimum solution.

3.4. Inversion Analysis Procedures

The inversion analysis approach and procedures can be established as follows based on the objective function, the type of parameters to be inversed, the AGA, and the coupled characteristics between the stress and fluid flow of rock masses:

(1) The parameters to be inversed are selected according to actual engineering conditions.

(2) The variation range of the parameters to be inversed is defined according to test results or the engineering experience.

(3) The AGA is adopted to transform the parameters into coded strings. Meanwhile, n initial groups are generated at random within the preset parameter scope.

(4) The initial solvable group is encoded into 3D back analysis programs of the coupled stress and fluid flow, and the objective function value F is computed.

(5) The corresponding objective function value F is estimated with the preset convergence condition.

(6) The optional, cross, and mutation operations are performed, and a new generation of solvable groups is formed according to the operational approach of the AGA.

(7) Steps (3) to (5) are repeated until the objective function F meets the error precision.

(8) The optimized objective function value and optimized inversion parameters are obtained.

The procedure for the inversion of the parameters is shown in Figure 1. The back analysis program of the coupled stress and fluid flow (3D-BackCSF.For) is compiled with Fortran according to Figure 1.

Figure 1: Flowchart of back analysis.

[figure omitted; refer to PDF]

4. Case Study

4.1. Engineering Introduction

Highway #1 of a hydropower station is situated in the Yajiang County belonging to a deep-cutting high mountainous area in the West Sichuan Plateau in Sichuan Province, east of the Tibetan Plateau. The overall terrain is high in the east and west and low in the middle. A highway tunnel is the main channel connecting both banks of the Yalong River.

Highway tunnel #1 for hydropower traffic engineering is 5855 m through the Xiala Mountain. The tunnel pile number is K7+480-K13+335. Rock masses of the tunnel area of the highway tunnel #1 are mainly distributed in the middle of the Lianghekou group (T3In2). The rocks are dark-gray to gray-black sandy carbonaceous slate and metamorphic quartz sandstone, which are developed in layers, and the joints are undeveloped. The main mineral composition of the rocks is feldspar and quartz stone. No large fault and fault structure exist in the engineering area, where the rock masses mainly contain structural joints. Groundwater is classified into quaternary loose accumulation layer pore water and bedrock fissure water. The rechargeable sources of the quaternary loose accumulation layer pore water mainly comprise surface runoff and atmospheric precipitation. By contrast, the recharge source of bedrock fissure water is far less diversified. Most of the water comprises atmospheric precipitation and a small amount of melted snow water from high mountains. The Lianghekou Tunnel has a simple monoclinal structure composed of hard fine-grained quartz and sandstone. The rock masses are broken without bad geological and special strata. The surrounding rock of the tunnel comprises Lianghekou group weak weathering fine-grained quartz sandstone, which is mainly classified as grade IV.

The pile number of field test section in the tunnel is K8+450-K8+600, whereas the average buried depth ranges from 595 m to 603 m. Rich groundwater is contained in the rock formation because of abundant rainfall in the area. According to the engineering geological investigation report, the buried depth of groundwater in the test section is approximately 127-135 m, and the groundwater is located in the weakly weathered layer. The rock masses in the area of highway tunnel #1 are mainly distributed in the middle of the Lianghekou group (T^3In2 ), which are dark-gray to gray-black sandy carbonaceous rocks. The grade of the surrounding rock in the test section tunnel is IV. The tunnel was built by the benching tunneling method. The preliminary support parameters are as follows: the diameter of the anchor is [varphi]25, and the length is 4.0 m. The longitudinal seam and circumferential spacing are 1.2 m × l.2 m. The anchor wall is constructed by [varphi]8 (20 cm × 20 cm) reinforcing mesh. C25 shotcrete with a thickness of 20 cm is used in the preliminary support structure. The waterproof board is 1.2 mm thick. The steel frame implements are I16 H and @80. The mechanical parameters of the shotcrete and grouting anchor layer are listed in Table 2. The construction method is shown in Figure 2.

Table 2: Physicomechanical parameters of shotcrete and grouting anchor layer.

Area	Gravity γ (KN/m3 )	Modulus of elasticity E (GPa)	Poisson ratio μ	Cohesion C (MPa)	Internal friction angle [straight phi] (°)
C25 shotcrete layer	22	24	0.23	1.6	34
Grouting anchor layer	23	22	0.25	2.0	36

Figure 2: Benching construction method.

[figure omitted; refer to PDF]

4.2. Parameters to Be Inversed

The effects of high water pressure on tunnel stability should be considered because the groundwater level is high. The buried depth is high and the ground stress is huge because the area is located in the deep-cutting high-mountain area to the east of the Tibetan Plateau. Therefore, the effects of the initial ground stress on the tunnel stability should also be considered. According to the aforementioned analysis, the coupled function should be considered in the tunnel.

In this project, the initial ground stress is high because the buried depth is about 600 m. Therefore, determining the lateral pressure coefficients _λx and _λy of the initial ground field in the tunnel area is important. Meanwhile, the effects of groundwater on the surrounding rock stress are the key to analyzing the seepage-stress coupling. The initial permeability coefficient _K0i is selected as the inversed parameter. Given that the groundwater is located in the weakly weathered sandstone layer, the permeability characteristics of the weakly weathered sandstone layer, metamorphic sandstone layer, and metamorphic quartz sandstone layer are relatively close on the basis of the engineering experience. Therefore, the initial permeability coefficients of the three strata are assumed to be identical in the process of the inverse analysis. The modulus of elasticity E, Poisson ratio μ, and frication angle [straight phi] have already been measured by the laboratory test and will not be inversed this time.

4.3. Numerical Model and Boundary Condition

4.3.1. Finite Element Model

Tunnel pile number K8+500-K8+550 is employed as the test segment. A geological survey shows that the average buried depth of this tunnel in vertical direction is approximately 600 m, and the average buried depth of the groundwater is approximately 132 m from the tunnel vault. The strata are generally divided into five layers from top to bottom: the strongly weathered sandstone layer, the weakly weathered sandstone layer, the metamorphic sandstone layer, the metamorphic quartz sandstone layer, and the weakly weathered fine quartz sandstone layer. Table 3 shows the mechanical parameters of each geological layer.

Table 3: Physicomechanical parameters of the stratum.

Stratum	Gravity γ (kN/m3 )	Modulus of elasticity E (GPa)	Poisson ratio μ	Cohesion C (MPa)	Internal friction angle [straight phi] (°)
Strongly weathered sandstone	22	5.1	0.33	0.9	20
Weakly weathered sandstone	22	12.5	0.30	1.2	23
Metamorphic sandstone	22	15.2	0.28	1.5	25
Metamorphic quartz sandstone	24	19.0	0.29	1.6	32
Weakly weathered fine quartz sandstone	24	21.0	0.28	1.5	35

A total of 83,176 nodes and 462,398 tetrahedron elements exist in the calculation area. The tetrahedron elements include 1,816 nodes for the surrounding rock of grouting anchor layer, which are divided into 6,237 elements, and 1,632 elements for the C25 sprayed concrete layer, which are distributed into 3,200 elements. Figure 3(b) shows the FEM.

Figure 3: Numerical models.

(a) Physical model

[figure omitted; refer to PDF]

(b) Finite element model

[figure omitted; refer to PDF]

In the whole coordinate system, the coordinate origin is O, which is located 100 m below the tunnel vault. The width of 50 m is chosen along the left and right of the calculation area. The actual calculation area is 244 m × 100 m × 50 m. Direction X is perpendicular to the tunnel axis. Direction Z is the vertical direction, and forward denotes the positive direction. Direction Y is consistent with the direction of the tunnel axis (see Figure 3(b)). The engineering geological condition is divided into four layers in the physical model from top to bottom: the weakly weathered sandstone layer, the metamorphic sandstone layer, the metamorphic quartz sandstone layer, and the weakly weathered fine quartz sandstone layer. The different strata thicknesses are as follows: 60 m for the weakly weathered sandstone layer, 50 m for the metamorphic sandstone layer, 64 m for the metamorphic quartz sandstone layer, and 70 m for the weakly weathered fine quartz sandstone layer. The thickness of the grouting anchor layer is 4 m. The calculation area is shown in Figure 3(a).

4.3.2. Boundary Condition

(1) The geological investigation shows that, in the area of tunnel pile number K8+500-K8+550, the groundwater level measured above the surrounding rocks is approximately 132 m, and the groundwater is located in the weakly weathered sandstone layer (see Figure 3(a)). According to the actual construction condition, a large amount of groundwater should have been infiltrated. However, given that the new Austrian method is adopted, the tunnel grouting sealing and preliminary support structure are performed timely. Therefore, the change of groundwater level will be ignored, and the buried depth of the groundwater is assumed to be 132 m.

(2) In the seepage field, the free surface above the tunnel is 132 m, the computation region is assumed to be a plane, the profile surface inside the tunnel is a seepage-free surface, and the pore water pressure is 0. The OXY at the bottom, left, and right sides are all impervious boundaries.

(3) The OXY at the bottom is far from the tunnel excavation face. Thus, the displacement of the OXY plane at these three directions is assumed to be 0. The displacement at direction X in the YZ plane at the left and right sides is also supposed to be 0. Meanwhile, if K8+500 is used as the start point, then the Y=0 m profile along the axis direction is the common plane of the front and rear construction sections. The deformation is stable because the tunnel has already been excavated and the preliminary support structure was made in the previous period. Therefore, the displacement of this section at direction Y is assumed to be 0. For the Y=50 m section along the axis direction, the rear surrounding rock is not excavated and the displacement of the section at direction Y is also assumed to be 0.

(4) The initial ground stress at direction Z is generated by the dead weight of rock masses, and the initial ground stresses at directions X and Y are related to buried depth H and the tectonic stress. The main objective is to identify the lateral pressure coefficients. To simplify the analysis, we assume that the lateral pressure coefficients of different buried depths were the same in this model. Initial ground stresses at directions X and Y are assumed to be _λx γH and _λy γH, respectively.

4.3.3. Objective Function

Considering that the hydraulic head is constant, measuring displacements are only adopted as compared information. Therefore, the objective function of back analysis is constructed as follows: [figure omitted; refer to PDF]

4.3.4. Water Pressure Analysis

When the groundwater effect is analyzed, the hydrostatic pressure and seepage hydrodynamic pressure should be considered. The seepage hydrostatic pressure p is [figure omitted; refer to PDF]

The seepage hydrodynamic pressure f[arrow right] is [figure omitted; refer to PDF] where _fx , _fy , and _fz are the components of the seepage body forces in directions x, y, and z, respectively.

In the equivalent continuous medium coupling model, the relational expression between the seepage field and stress field is considerably important. In this analysis, the relational expression between the stress and permeability coefficient is used as the relation equation, and the parameters are inversed by using the equivalent continuous coupled model. The relational expression of the coupled stress and fluid flow adopts a Louis empirical formula, as shown in the following equation: [figure omitted; refer to PDF] where _K0i are the initial permeability coefficients of the X, Y, and Z directions; μ is the Poisson ratio; ρ is the rock mass density; _σe is the effective normal stress and its direction is perpendicular to the main permeability direction; P is the seepage pressure; α is the coupling coefficient between the seepage and stress fields; and 1.0 is the calculated value.

4.4. Tunnel Deformation Monitoring Scheme

Tunnel vault settlement and clearance convergence measurements are conducted in test section K8+500-K8+550. A vault settlement point and a clearance convergence line are installed in one section. The vault settlement point is A, whereas the clearance convergence measuring point is BC, as shown in Figures 4 and 5.

Figure 4: Layout of deformation monitoring points.

(a) Measuring point of vault settlement

[figure omitted; refer to PDF]

(b) Measuring point of clearance convergence

[figure omitted; refer to PDF]

Figure 5: Measuring point and measuring line arrangement.

[figure omitted; refer to PDF]

The vault settlement is measured using a PENTAX R-322 total station with a precision of 1 mm + 1 ppmm, whereas the clearance convergence is measured using a JTM-J7100 steel rule convergence gauge with a precision of 0.01 mm. The deformation of six sections from K8+500 to K8+550 is monitored, and the results are shown in Table 4.

Table 4: Measured and fitted values of measuring points.

Pile number section	Measured vault settlement value (mm)	Fitted vault settlement value (mm)	Measured clearance coverage value (mm)	Fitted clearance convergence value (mm)
K8+500	50.1	51.1	44.4	45.9
K8+510	50.5	51.6	42.1	46.9
K8+520	50.4	52.2	47.9	47.7
K8+530	48.9	52.7	50.7	48.3
K8+540	54.8	53.1	44.8	46.1
K8+550	52.4	53.5	43.2	45.6

4.5. Back Analysis of Parameters

The 3D-BackCSS back analysis procedure is compiled and used in this paper. The pore water pressure, normal stress _σx in direction X, and normal stress _σz in direction Z of each tunnel section are obtained through the coupled analysis, as shown in Figures 6 and 7. Although the analysis results of the five sections are different, the results of the section of each pile number have a similar tendency. Therefore, only sections K8+510 and K8+550 are used to reveal the law.

Figure 6: Contour map of pore water pressure/MPa.

(a) Pore water pressure of contour map of section K8+510/MPa

[figure omitted; refer to PDF]

(b) Pore water pressure of contour map of section K8+550/MPa

[figure omitted; refer to PDF]

Figure 7: Contour map of the stresses.

(a) Contour map _σZ of section K8+510/MPa

[figure omitted; refer to PDF]

(b) Contour map _σZ of section K8+550/MPa

[figure omitted; refer to PDF]

(c) Contour map _σx of section K8+510/MPa

[figure omitted; refer to PDF]

(d) Contour map _σx of section K8+550/MPa

[figure omitted; refer to PDF]

Figure 6 shows that the internal profile of the tunnel is a free surface and that the pore water pressure is 0 MPa. In the tunnel profile line with a 30 m × 56 m (height × width) area, the hydraulic head is regular in the circular ring. Meanwhile, the maximum hydraulic head occurs at the bottom of the calculation area. This trend is consistent with theory and engineering practices, and the size of pore pressure is reasonable.

The maximum normal stress _σZ of the two sections in direction Z occurs near the vault of the tunnel, as shown in Figure 7. The vault and inverted arch with the largest deformation have the minimum stress, and _σZ is approximately 0.5-0.6 MPa. The location with the largest deformation has a smaller stress value because part of the stress is released. In a certain range of the tunnel portal, especially near the bottom arch, stress concentration is likely to occur as a result of large deformation constraint. Therefore, the stress is large, which is consistent with the actual condition.

After 180 iterations of the optimization calculation, if the value of objective function is 0.0318, then the measured displacement and fitting displacement deviation of the measuring point reach the minimum. Given that the model is established with a large height and width, all calculation results of this project indicate that the range is 30 m × 56 m (height × width) for better analysis and comparison of the numerical calculation results. The results are presented in Figure 8. The results show that the maximum settlement of the tunnel occurs at the vault and that the displacement value is approximately 50 mm. The maximum displacement in direction X occurs in the part of the tunnel with the maximum width, and the maximum horizontal displacement reaches 48.6 mm.

Figure 8: Contour map of the displacements.

(a) Displacement contour map of section K8+510 in direction Z/m

[figure omitted; refer to PDF]

(b) Displacement contour map of section K8+550 in direction Z/m

[figure omitted; refer to PDF]

(c) Displacement contour map of section K8+510 in direction X/m

[figure omitted; refer to PDF]

(d) Displacement contour map of section K8+550 in direction X/m

[figure omitted; refer to PDF]

Deviation between the measured and fitted value can be determined by comparing the measured vault settlement value and clearance convergence value with the fitted vault settlement and clearance convergence value in the same section. The details are shown in Table 4. As listed in Table 4, section K8+530 has the largest error between the measured and fitted vault settlement values, with an error of 3.8 mm and a relative error rate of 2.2%. Section K8+510 has the smallest error between the measured and calculated vault settlement values, with an error of 1.1 mm and a relative error rate of 2.2%. Section K8+510 has the largest clearance convergence displacement error of 4.8 mm and a relative error rate of 11.4%, whereas section K8+520 has the smallest clearance convergence displacement error of 0.2 mm and relative error rate of 0.4%. In general, the average error is below 5%, and the result is rational.

After optimization analysis, the error between measured displacement and fitted displacement of the measuring point reaches the minimum when the value of objective function is 0.0318. The inverse parameters are identified, and the result is shown in Table 5. After comparing the measured and fitted displacement values of the five pile sections, aside from individual values, the law is obvious: the fitted values of the tunnel vault and fitted value of clearance convergence value are larger than the measured values. In fact, this result reflects a problem in the displacement measurement process. During the tunnel excavation process, a time difference exists between the completion of excavation and the first measurement. Generally, the time difference exceeds 6 hours. The displacement will occur in this time. In addition, the excavation causes the deformation of surrounding rocks before the excavation of the tunnel face. Therefore, the fitted values are always larger than the measured values. From the measuring process, the tunnel has a rapid deformation rate in the early phase. During this period, the part of the deformation occurs in a short period. However, the deformation monitoring lags behind the actual deformation. Therefore, the actual deformation is larger than the measured deformation, but estimating the specific value is challenging. The calculation and analysis are performed according to the surrounding rock condition, relative mechanical parameters, and stress condition, with the fitted tunnel deformation certainly larger than the measured value. Thus, the parameters obtained from back analysis are greater than the measured value.

Table 5: Optional results of the back analysis.

Lateral pressure coefficient λx	Lateral pressure coefficient λy	Initial permeability coefficient
K 0 x (m/s)	K 0 y (m/s)	K 0 z (m/s)
1.15	0.76	7.8 × 10-7	1.3 × 10-7	5.3 × 10-6

5. Conclusions

The equivalent continuous coupled model was adopted in this study by considering the hydrostatic pressure and seepage pressure. The error between the measured and fitted values of the vault settlement and clearance convergence was adopted as the objective function. Meanwhile, the back analysis was conducted and the initial permeability coefficients and lateral coefficients of initial ground stress are identified by combining the FEM with the AGA. The main conclusions are as follows:

(1) The inversion model and method for the coupled stress and fluid flow analysis were established by combining FEM with AGA. The initial permeability coefficients of the stratum and the lateral pressure coefficients of the initial ground stress were identified by the back analysis. This algorithm overcomes the limitations of the traditional optimization algorithms, in which the inversion result significantly depends on the initial value of _Pc and _Pm and tends to fall into the local optimum when the inversion analysis is performed.

(2) The inverse model was applied in Lianghekou highway tunnel. According to the measured results of tunnel test sections, the relative error of vault settlement in five sections ranged from 2.2% to 7.7%, whereas the relative error of clearance convergence ranged from 0.4% to 11.4%. In general, the average error was below 5%, and the error was small. The results show that the established inversion analysis method and model are effective.

(3) After a 180-iteration optimization calculation, the optimal parameter combination of the initial permeability coefficients and the lateral pressure coefficients was obtained when the objective function value F was 0.0318. According to the actual conditions, the inversion results are reasonable. The proposed method takes advantage of the AGA and has both the speed and precision of the back analysis.

Acknowledgments

This research was financially supported by the National Natural Science Foundation of China (Grant no. 51408054 sponsored), the Scientific Research Foundation (15JK1337) by the Education Department of Shaanxi Province, the Scientific Research Program (KLTLR-Y14-15) for Technology of Highway Construction and Maintenance Technology of National Transportation Industry Key Laboratory, and the Research Program (XAGDXJJ16003) sponsored by Xi'an Technological University.

References

[1] Z. Q. Liu, Q. Zhang, "A review on the state of art of the saturated seepage-stress coupling models in rock mass,", Advances in Mechanics , vol. 38, no. 5, pp. 585-600, 2008.

[2] D. T. Snow, "Anisotropie permeability of fractured media,", Water Resources Research , vol. 5, no. 6, pp. 1273-1289, 1969.

[3] C. Louis, Rock Hydraulics in Rock Mechanics , Springer, New York, NY, USA, 1974.

[4] R. Nelson, Fracture Permeability in Porous Reservoirs: Experimental and Field Approach , Texas A&M University, College Station, Tex, USA, 1975.

[5] Y. Z. Zhang, J. C. Zhang, "Experimental study of the seepage flow-stress coupling in fractured rock masses,", Rock and Soil Mechanics , vol. 18, no. 4, pp. 59-62, 1997.

[6] B. Y. Su, M. L. Zhan, Y. Wang, "Experimental study of the characteristic of seepage flow-stress coupling in fracture,", Rock and Soil Mechanics , vol. 19, no. 4, pp. 73-76, 1997.

[7] Z. Chang, Y. Zhao, Y. Hu, D. Yang, "Theoretic and experimental studies on seepage law of single fracture under 3D stresses,", Chinese Journal of Rock Mechanics and Engineering , vol. 23, no. 4, pp. 620-624, 2004.

[8] Z. X. Chang, Y. S. Zhao, Y. Q. Hu, D. Yang, "Theoretic and experimental studies of the coupling of seepage flow and 3D stresses in fractured rock masses,", supplement 1, Chinese Journal of Rock Mechanics and Engineering , vol. 23, pp. 4907-4911, 2004.

[9] S. E. Minkoff, C. M. Stone, S. Bryant, M. Peszynska, M. F. Wheeler, "Coupled fluid flow and geomechanical deformation modeling,", Journal of Petroleum Science and Engineering , vol. 38, no. 1-2, pp. 37-56, 2003.

[10] A. Kobayashi, T. Fujita, M. Chijimatsu, "Continuous approach for coupled mechanical and hydraulic behavior of a fractured rock mass during hypothetical shaft sinking at Sellafield, UK,", International Journal of Rock Mechanics and Mining Sciences , vol. 38, no. 1, pp. 45-57, 2001.

[11] G. Cammarata, C. Fidelibus, M. Cravero, G. Barla, "The hydro-mechanically coupled response of rock fractures,", Rock Mechanics and Rock Engineering , vol. 40, no. 1, pp. 41-61, 2007.

[12] Y. W. Tsang, C. F. Tsang, "Channel model of flow through fractured media,", Water Resources Research , vol. 23, no. 3, pp. 467-479, 1987.

[13] National Research Council, Rock Fractures and Fluid Flow , National Academy Press, Washington, DC, USA, 1996.

[14] L. Jing, Y. Ma, Z. Fang, "Modeling of fluid flow and solid deformation for fractured rocks with discontinuous deformation analysis (DDA) method,", International Journal of Rock Mechanics and Mining Sciences , vol. 38, no. 3, pp. 343-355, 2001.

[15] L. Jing, J. A. Hudson, "Numerical methods in rock mechanics,", International Journal of Rock Mechanics and Mining Sciences , vol. 39, no. 4, pp. 409-427, 2002.

[16] T. Vogel, H. H. Gerke, R. Zhang, M. T. Van Genuchten, "Modeling flow and transport in a two-dimensional dual-permeability system with spatially variable hydraulic properties,", Journal of Hydrology , vol. 238, no. 1-2, pp. 78-89, 2000.

[17] J. F. Shao, H. Zhou, K. T. Chau, "Coupling between anisotropic damage and permeability variation in brittle rocks,", International Journal for Numerical and Analytical Methods in Geomechanics , vol. 29, no. 12, pp. 1231-1247, 2005.

[18] S. C. Yuan, J. P. Harrison, "Development of a hydro-mechanical local degradation approach and its application to modelling fluid flow during progressive fracturing of heterogeneous rocks,", International Journal of Rock Mechanics and Mining Sciences , vol. 42, no. 7-8, pp. 961-984, 2005.

[19] K. Maleki, A. Pouya, "Numerical simulation of damage--permeability relationship in brittle geomaterials,", Computers and Geotechnics , vol. 37, no. 5, pp. 619-628, 2010.

[20] J.-M. Pereira, C. Arson, "Retention and permeability properties of damaged porous rocks,", Computers and Geotechnics , vol. 48, pp. 272-282, 2013.

[21] C. Arson, J.-M. Pereira, "Influence of damage on pore size distribution and permeability of rocks,", International Journal for Numerical and Analytical Methods in Geomechanics , vol. 37, no. 8, pp. 810-831, 2013.

[22] M. L. De Bellis, G. Della Vecchia, M. Ortiz, A. Pandolfi, "A linearized porous brittle damage material model with distributed frictional-cohesive faults,", Engineering Geology , vol. 215, pp. 10-24, 2016.

[23] G. Gioda, G. Maier, "Direct search solution of an inverse problem in elastoplasticity: identification of cohesion, friction angle and in situ stress by pressure tunnel tests,", International Journal for Numerical Methods in Engineering , vol. 15, no. 12, pp. 1823-1848, 1980.

[24] S. Sakurai, K. Takeuchi, "Back analysis of measured displacements of tunnels,", Rock Mechanics and Rock Engineering , vol. 16, no. 3, pp. 173-180, 1983.

[25] J. Sheng, B. Su, Y. Wang, M. Zhan, "Coupling analysis of elasto-plastic stress and fluid flow in jointed rock masses,", Chinese Journal of Rock Mechanics and Engineering , vol. 19, no. 3, pp. 304-309, 2000.

[26] Y. Wang, J. Liu, "Parameter inversion for fully coupled problem of steady fluid flow and stress in fractured rock masses based on sensitivity analysis,", Rock and Soil Mechanics , vol. 30, no. 2, pp. 311-317, 2000.

[27] L. Guo, X. Li, Y. Zhou, Y. Zhang, "Generation and verification of three-dimensional network of fractured rock masses stochastic discontinuities based on digitalization,", Environmental Earth Sciences , vol. 73, no. 11, article no. 29, pp. 7075-7088, 2015.

[28] Y. Q. Wu, "Multi-level fracture network model and FE solution for ground water flow in rock mass,", Journal of Hydraulic Research , vol. 43, no. 2, pp. 202-207, 2007.

[29] J. Zhang, J. Wang, "Coupled behavior of stress and permeability and its engineering applications,", Chinese Journal of Rock Mechanics and Engineering , vol. 25, no. 10, pp. 1981-1989, 2006.

[30] T. H. Yang, P. Jia, W. H. Shi, P. T. Wang, H. L. Liu, Q. L. Yu, "Seepage-stress coupled analysis on anisotropic characteristics of the fractured rock mass around roadway,", Tunnelling and Underground Space Technology , vol. 43, pp. 11-19, 2014.

[31] J. Chai, Y. Wu, S. Li, "Analysis of coupled seepage and stress fields in rock mass around the Xiaowan arch dam,", Communications in Numerical Methods in Engineering , vol. 20, no. 8, pp. 607-617, 2004.

[32] Y. Q. Wu, Z. Y. Zhang 1995 (Chinese), An Introduction to Rock Mass Hydraulics , Southwest Jiaotong University, Chengdu, China

[33] F. Cappa, Y. Guglielmi, P. Fénart, V. Merrien-Soukatchoff, A. Thoraval, "Hydro-mechanical interactions in a fractured carbonate reservoir inferred from hydraulic and mechanical measurements,", International Journal of Rock Mechanics and Mining Sciences , vol. 42, no. 2, pp. 287-306, 2005.

[34] M. Wu, P. Chen, J. Chai, "Coupled unsteady 3D seepage and stress fields in fracture rock mass,", Journal of Hydroelectric Engineering , vol. 29, no. 6, pp. 199-204, 2010.

[35] X. Deng, H. Fang, J. Chen, J. Kou, "Identification of hydraulic conductivity in aquifer for coupled FEM and adaptive genetic algorithm,", Mathematical Problems in Engineering , vol. 2015, 2015.