Abstract: In this research, the hydraulic characteristics of a horizontal disc, which is assumed as a rock joint, have been studied analytically and numerically with reference to radial flow occurred along a typical water pressure test that is also known as Lugeon test. The hydraulic characteristics such as Lugeon value, i.e., flowrate to hydrostatic pressure ratio, the joint hydraulic aperture, inertia conditions, turbulence and laminar conditions and the effect of relative roughness were investigated by use of well-known hydraulic relations coincided by numerical analysis. The numerical results were regarded as a verifier to evaluate the accuracy of empirical and analytical relations. The comparison between analytical and numerical results was based on Lu value which is recorded as the main result of Lugeon test. The results show that the simplest equation (the Darcy equation of radial flow) demonstrates the closer estimation to numerical one, in comparison with other analytical and empirical equations. Also, the Lugeon test must be considered as a single-well test and it could just identify a combination of some hydraulic parameters. A practical procedure is finally presented to calculate the possible hydraulic characteristics of tested section in the borehole according to Water Pressure test.
INTRODUCTION
Hydraulic characterization of rock joint is of importance in dam construction or tunneling especially for evaluation of seepage and groutability. The characterization mainly describes the flow fashion of joint as a conductive duct while the water is passed through it via hydraulic tests. The results of characterization could be reckoned as the introductory information for analyzing the hydraulic properties of joints network and finally rock mass in site.
Water Pressure Test (WPT or Lugeon test), which is a typical means to try out these characteristics, has been introduced by Morris Lugeon as a customary competitive task; nowadays is done throughout the site investigations. Through this well-known test, water with various pressures is injected to a packed section of borehole and flowrate is being recorded. Finally two decisive parameters flowrate, Q and effective pressure, P, are plotted and interpreted as the hydraulic behavior of tested section and the Lugeon value, Lu, calculated regarding the P-Q diagram.
Hydraulic characterization of rock discontinuities via WPT could be generally carried out after recording P-Q values, classification and interpretation of P-Q diagram and ascertaining Lugeon value. It must be considered that WPT is a single-well test and therefore could not exactly describe the hydraulic conditions of host rock, However, there are noticeable developments in analytical and empirical studies of the hydraulic characterization of rock joints using WPT. Using the Darcy assumptions, The basic equation of radial flow through a rock fracture for WPT situation, which is quite similar to the Thiem equation for steady state flow in a confined aquifer during a pump test (Swhartz and Zhang, 2003; Todd and Mays, 2005), would be extracted. Beyond the basic Darcy assumptions, a practical study on the hydraulic characteristics of a horizontal single joint regarding the WPT conditions are done by Elsworth and Doe (1986) based on Loius’ empirical results in relation to different conditions of flow. Some empirical modeling considering groutability of rock and hydrogeological properties of a rock fracture related to hydraulic tests is reported by Fransson (2001, 2002, 2004), Fransson and Gustafson (2000). Barton and Quadros (1997, 2003); also, Wittke (1990) has presented and discussed precisely about rock hydraulics on the subject of radial flow of injected water occurred during WPT and seepage flow through rock mass as well.
The situation regarded by Elsworth and Doe (1986) for developing the analytical formulation of radial flow was hard to achieve in a typical test. Also the formulation is not tested numerically for steady state. Fransson (2001, 2002, 2004) and Fransson and Gustafson (2000) have just considered laminar state of flow based on Thiem equation. Barton and Quadros (1997, 2003) have concentrated on the experimental results about joint roughness coefficient especially in laminar flow and did not verify the complex condition of turbulence situation. Some experimental formula also has been cited by Wittke (1990) but no numerical verification has been given.
In this research, the hydraulic characteristics of flow in a horizontal discontinuity are described with regard to Lu value as the main result of WPT. The flow is assumed to be radial toward a parallel horizontal disc. The horizontal smooth disc is considered as the simplest model that a joint would be had. The flow conditions in this disc are studied numerically and by use of analytical equations to characterize the hydraulic circumstances and also to verify the experimental equations. In this model, flow with a given pressure or a given flowrate is entered from cylindrical entrance that is taken into account as the joint orifice in the intersection area with borehole wall. Finally, the observed discrepancies between analytical and numerical results are discussed and a flowchart is proposed to characterize water pressure test results.
BASIC FORMULATIONS OF LOW-VELOCITY LAMINAR FLOW
The outward radial water flow between two parallel discs from a cylindrical inlet with radius ri toward a cylindrical outlet with the radius rf is considered initially. Assuming:
(1) |
and according to continuity equation as follows
(2) |
the Navier-Stokes equation is simplified as
(3) |
where, V is velocity magnitude, Vr is the radial velocity magnitude as a component of velocity vector f(z) is a simplified form of independent function describes velocity changes along z direction, r, z and θ are the components of cylindrical coordinate, vector nr is the radial unit vector, h is the hydrostatic pressure and g, ρw and μw are, respectively gravity acceleration, specific mass and dynamic viscosity of water.
The differential Eq. 3 consists of three terms Hydrostatic, dh/dr, which is equal to hydraulic gradient, i, Inertia, Vr2/gr, which is also deemed as kinematics term and Viscous, i.e., the right side of the Eq. 3. As the radius of flow increases, the radial velocity diminished according to Eq. 2 and the inertia term may be vanished. It may simplify the Eq. 3 which is finally solved as follows:
(4) |
(5) |
(6) |
(7) |
where, Vr* is the average velocity, b is the aperture of disc which is also defined as hydraulic aperture, hi and hf are the initial and final hydrostatic heads of flow, ri and rf are the initial and final flow radii, Q is the flowrate and i is hydraulic gradient.
Some parameters mentioned above are shown in Fig. 1.
| Fig. 1: | The schematic configuration of simplified rock joint
crossing the borehole used in analytical and numerical calculation.
(a) hydraulic aperture, (b) flow direction, (c) final radius, (d)
initial radius and (e) general trend of hydrostatic pressure or velocity
along radial flow direction |
The Eq. 5 is a highlighted equation at radial flow where both conditions of laminar and low-velocity are satisfied, which is on the base of Darcy assumptions. So, it would be introduced as Darcy equation of redial flow in this research. It is equivalent to Thiem equation which is considered for pumping water wells in the confined aquifer in steady state (Swhartz and Zhang, 2003; Todd and Mays, 2005).
Effective radius: The final radius of flow, rf, would be replaced by effective radius, re (in which the hydrostatic head drops to zero or to groundwater head) and the initial radius, ri, would be replaced by borehole radius, rw. Unfortunately, the effective radius may not be estimated because the WPT is a single-well test and two parameters flowrate and static pressure are merely determined all along the test. But, since the neperian logarithm of re/rw exists in such equations, the effect of uncertainty of re decreases. In the Table 1 some estimations of re are presented. The relations mentioned in Table 1 are extracted from the experimental equations proposed by some researchers to calculate the permeability. However, it is better that the unknown term of ln(re/rw) is expressed as boundary coefficient, ζ, to show the uncertainty. Since both rw and re are the properties of inlet and outlet boundaries, ζ is named as boundary coefficient. Also, some hydrogeological parameters such as joint transmissivity, T, would be described as T/ζ or Tζ (i.e., the normalized transmissivity).
Hydraulic aperture, permeability and transmissivity: Lugeon value, Lu, is a representative parameter with reference to permeability of tested section and is defined as follows:
(8) |
where, P is the effective hydrostatic pressure of water at joint entrance and L is the tested section length. In terms of Lugeon value, the Eq. 5 would be changed as follows:
| Table 1: | Some proposed alternative empirical equations to calculate the effective radius and ζ (After Ewert, 1985; Hamm et al., 2007) |
| *e is the neperian constant | |
(9) |
where, bζ is considered as normalized aperture with reference to the boundary coefficient. Note that all parameters except for Lu, which is presented in L min-1. MPa, must be entered in SI units in Eq. 9 and any equation that concerns Lu value here.
The equation of Darcy can be expressed in the form of V = K.i. Therefore, according to Eq. 4 and 6 the joint permeability is:
(15) |
where, Kζ is normalized permeability considering the boundary coefficient. If transmissivity is defined as T = K.b, it would also be calculated by use of the following equation:
(16) |
INERTIA AND TURBULENCE LIMITATIONS
The Eq. 3 has been solved with respects to some basic assumption including simple geometry (very low relative roughness of joint surfaces) and laminar flow with low velocity so that the inertia term is neglected. Actually, in a natural joint and actual water pressure test, the velocity of injected water at least around the joint entrance may be so noticeable that the inertia term would not be eliminated and the joint wall roughness may be so intensive that the flow behavior is changed especially near the wall. Another complicated situation that may occur specially at joint entrance is the turbulence condition depends upon both velocity intensity and the aperture magnitude. As a result, the Eq. 4 to 6 and subsequences do not satisfy the flow situation where the wall roughness and flow velocity are considerable.
Verification of inertia: Regarding Eq. 3, the ratio between inertia term and hydraulic gradient is interpreted as an index showing the importance of inertia term. Using some simplifications, this ratio, which is defined as the relative inertia intensity at the radius r, is given by:
(17) |
where, (Vr2/gr)* recognized as average inertia in direction of flow and almost-equal sign refers to the basic assumptions resulted in Eq. 5.
Elsworth and Doe (1986) have conducted some actual tests and have uttered that for values of ç greater than 0.5 the inertia term should not be neglected. By defining a threshold as critical relative intensity of inertia (ηc = 0.5), the Darcian radius, rD, is expressed by:
(18) |
If the Darcian radius is smaller than borehole radius (rD≤ rw) the inertia term is could be neglected. If rD is greater than rw, the inertia term must be taken into account in formulations. Maini (1971) has proposed an experimental equation where the inertia term is considerable.
Hydraulic aperture considering inertia term: The hydraulic aperture could be re-calculated using experimental equation of Maini (1971) by use of Lu values as follows:
b3-c1b-c2
= 0 |
(19) |
Where:
(20) |
and
(21) |
Equation 6 could be simply solved. Results have indicated that this calculated hydraulic aperture is a bit more than when the inertia term is neglected.
Reynolds number, critical flowrate and turbulence radius: The Reynolds number is an appropriate parameter to identify if turbulence condition exists considering both velocity (in terms of flow) and aperture (in terms of geometry). For a tubular flow such as flow in pipes the Reynolds number is given by:
(22) |
where, the hydraulic diameter, D, is the pipe diameter. But, for non-tubular flow geometry such as a duct, it is defined as:
(23) |
where, S is the saturated section of flow and W is the wet perimeter of duct.
The hydraulic diameter is around 2b for the radial flow in joint. So, the Reynolds number could be expressed by:
(24) |
It is interesting to note that Re is directly independent of both aperture and velocity. Consequently, the Reynolds number can be calculated straight from flowrates which are recorded during Lugeon test.
According to Munson et al. (2002), for a round duct assumed here, if Re≤ 2100, the flow is laminar. For 2100≤ Re≤ 4000, the flow is in transient condition and for 4000≤ Re, the flow is completely turbulence. However, Elthworh and Doe (1986) and Wittke (1990) have affirmed that critical Reynolds is equal to 2300. Prudently, to preclude turbulence effects, it is better to choose the lowest empirical critical Reynolds for calculations. Given critical Reynolds number, Recr = 2100, critical flowrate could be ascertained as follows:
(25) |
The critical flowrate may be interpreted as the upper limitation for laminar regime. The Lugeon number that will be applied in the equations must be chosen or extrapolated from P-Q diagram where Q is less than Qcr. Given Recr = 2100 in the Eq. 25 the turbulence radius could be calculated as follows:
(26) |
It seems that, as Reynolds number drops below the Recr farther than rT, the turbulence intensity is inconsequential and the laminar condition retains.
Head loss due to joint entrance: Since the water is entered from borehole surroundings into joint, a head loss is occurred due to the curvature of flowline and the local disturbance. Although, the analytical solution of this situation is more complicated, but Wittke (1990) has presented an experimental two-dimensional statement to find the length of disturbance zone, χ (where flow reaches to a completely parabolic velocity profile) and the hydrostatic pressure loss, ΔPχ. As regards WPT and the problem geometry, the equations revealed by Wittke (1990) could be rewritten as follows:
(27) |
and
(28) |
The borehole radius, rw and the hydrostatic pressure, P, should be corrected as P-ΔPχ and rw + χ if χ and ΔPχ are considerable in respect to rw and P.
Mechanical aperture: The mechanical aperture (or mean aperture) is the average distance between joint walls. This parameter is heavily dependent on how the openings vary along the flow paths. A quantity such as roughness may describe the changeability of fracture surfaces in a joint to approach the mechanical aperture.
There are some parameters describing the roughness condition of joint surface. The height of asperity amplitude along joint wall refers to the roughness height, ε and the relative roughness of joint wall is defined as ε/D or εr. However, measurement of these parameters is time-consuming and expensive in core drilling especially in a typical site investigation program. But, the usual field parameter of joint roughness coefficient, JRC, would be approximately recorded by inspecting the joint surface along cores. A beneficial empirical relation presented by Olsson and Barton (2000) (based on Barton, 1982) may be used to estimate the mechanical aperture, em, considering the hydraulic aperture and JRC as follows:
(29) |
This equation is only valid where bm/b ≥ 1. The mechanical aperture shows that there are many positions at joint surface which may have openings more than the mechanical aperture and less than it as well. For example, that could be deemed as the reason when a typical coarse grout would penetrate into a joint in spite of the fact that the equivalent hydraulic or mechanical aperture is estimated very low.
COMPREHENSIVE EQUATION OF RADIAL FLOW
In the form of power function, a relation between velocity and head gradient which also known as Missbach equation is expressed by:
(30) |
As Elsworth and Doe (1986) cited, the robustness of this law has been approved via some pervious simulations. Also some researches such as those of Wittke (1990), Fransson (2001) and Qian et al. (2005) have somewhat confirmed this law.
| Fig. 2: | Different possible flow regions in connection with Reynolds
number and relative roughness relating Table 1 (Elsworth
and Doe, 1986) |
| Table 2: | Determination of transmissivity, T, velocity exponent, m and boundary coefficient, ζ which are used in the solution of Missbach law equation regarding the flow region demonstrated in Fig. 1 |
According to the radial flow circumstances and Lu value, the general solution of Missbach law could be expressed by:
(31) |
It must be noted that in the Eq. 31 transmissivity is still equal to K.b. Considering Re and relative roughness, five flow regions would be defined, as demonstrated in Fig. 1. Three parameters T, m and ζ are described in respect of the flow regions in Table 2.
For example, in flow region LS (Fig. 2), Eq. 31 is converted to Eq. 9. For laminar flow in a rough joint, if Eq. 35 and 36 are substituted in Eq. 31, following relation would be defined between hydraulic aperture and Lu:
(32) |
where, κ = 1 + 8.8εr1.5. This relation shows that actually the Lu value just determines a combination of aperture, relative roughness and boundary coefficient in the form of b/(κ.ζ1/3) and non of them could be determine exactly by use of such a single well test.
| Fig. 3: | Hydraulic aperture versus relative roughness in different
Lu values and various effective radii according to Eq.
32. The Lugeon values are in lit/min.m.MPa as usual, (a) Lu =
1, (b) Lu = 5, (c) Lu = 10, (d) Lu = 30, (e) Lu = 60 and (f) Lu =
100 |
In Fig. 3a-f, variation of hydraulic aperture versus relative roughness is plotted in some different Lu values. Figure 3 shows that how many possible combinations of (b, εr, re) would be inferred for an exact Lu value just in the laminar condition.
If the aperture calculated by Eq. 32 and the aperture calculated by Eq. 9 are, respectively defined as equivalent rough aperture, br and equivalent hydraulic smooth aperture (abridged to hydraulic aperture) a relation such as Barton’s empirical relation can be define as follows:
(33) |
NUMERICAL AND ANALYTICAL RESULTS
To assess the accuracy of some useful aforementioned analytical and experimental relations, a series of numerical models have been conducted. The numerical modeling are done for several flat discs with a constant inner radius 3.5 cm (like a borehole radius), constant outer radius 1 m (like a one-meter section that re is set to its length), varying hydraulic aperture, b = 0.05, 0.1, 0.2, 0.5, 1 and 2 mm and varying flowrate Q = 1, 2, 3, 5, 10 L min-1 (below the critical Reynolds number) and 15, 20, 25, 30 and 35 L min-1 (above the critical Reynolds number). The general shape of models and used grid are demonstrated in Fig. 4.
The analytical relations are also calculated simultaneously for the several cases mentioned above. The solution method for numerical calculations has selected according to Finite Volume Method (FVM). The models are developed on the base of 3-D fully Cartesian equations Navier-Stokes; moreover according to critical Reynolds number, the turbulence equations of flow are also included where it is necessary. There are 200000 hexahedral cells that distributed into disc to calculate flow parameters. The minimum and maximum sizes of cells are 1.411E-11 and 1.797E-7 m3, respectively. As it is shown in Fig. 4, to achieve the more accurate results, the closer cells to inlet, the smaller volume they have, because the discrepancies are more noticeable nearby the inlet and lowers quickly toward outlet. The results are detailed in Table 3 and Fig. 5a-h.
Recording of joint roughness in core logging is not practical and convenient during a typical site investigation program. Some presented equations such as Eq. 32 show that the joint roughness could not be measured only by use of Lugeon test. But, with regard to Table 1, assuming an experimental value for re (or ζ), an equivalent hydraulic aperture for an assumed horizontal smooth flat joint would be determined for a given Lu value of a tested section. Actually, the tested section is simulated by a flat horizontal smooth disc including the mentioned equivalent hydraulic aperture. The calculated aperture is useful at least to estimate permeability and transmissivity (for example in seepage analysis) and also to make an initial estimate of rock groutability.
| Table 3: | Numerical modeling results vis-à-vis analytical and experimental equations |
| m1 and m2 are coefficients for - Q in L min-1 and P in MPa. r2 is the regression correlation coefficient; LuN is the average Lugeon value given by numerical calculations; LuA is the Lugeon value given by Eq. 9; LuCAI is the Lugeon value given by Eq. 9 in which re is replaced by rD; LuIN is the Lugeon value given by Eq. 19 | |
| Fig. 4: | General configuration of numerical model. (a) Upper wall, (b) Lower wall, (c) Inlet, (d) Outlet, (e) Interior grids and (f) General direction of flow |
DISCUSSION
Since the numerical results were conducted in complete form of Navier-Stokes equations, they are used to verify the efficacy of analytical equations in various states of flow and different aperture values. The error threshold has been selected around 1E-3 and the models, even for highest velocity and aperture, are so simple because of the general shape of model, uniform boundary condition and smooth wall set for solutions. All solutions reached to converged stable state with decreasing error. According to the numerical simulation, these results would be deemed as an adequate verifier to asses the accuracy of cited analytical and empirical equations at least at the smooth condition.
The analytical results which are calculated by use of Eq.
9 in the laminar state fairly correspond to numerical ones. In the
worst condition the relative error reaches to 6.7% (in low velocity laminar
state) and 21.3% (in laminar and considerable inertia state) and in the
practical range of Lu value (between 1 and 100) the error is below 4.6%.
When the inertia is considerable (i.e., the rD≥ re),
the results indicated that Eq. 9 could be still useful
if it is modified by substituting rD for re in boundary
coefficient. It can reduce the error from maximum 21.3% to maximum 8.5%.
The empirical Eq. 19 did not match with the numerical results and showed a noticeable relative error from minimum 19% to maximum 43%. It seems that this equation has been achieved for some special experiments and the parameters have been individually sat according to series of experiments conducted by Maini (1971).
| Fig. 5: | The Lugeon value in the transient and turbulent conditions
(Re ≥ 2100) (a) bh = 0.01 mm, (b) bh = 0.05
mm, (c) bh = 0.1 mm, (d) bh = 0.2 mm, (e) bh = 0.22 mm, (f) bh = 0.5 mm, (g) bh = 1 mm and
(h) bh = 2 mm LuN is the Lugeon value given
by numerical modeling; LuCAT is the Lugeon value given
by Eq. 9 in which rw is replaced by rT;
LuA is the Lugeon value given by Eq. 9;
LuTR is the Lugeon value given by Eq. 31 considering Eq. 37 and 38; LuCTR is the Lugeon value given by Eq. 31 considering Eq. 37 and 38 in which re is replaced by rT |
The equation extracted for analytical solution in smooth transient-turbulent condition Eq. 31 regarding NS in Fig. 2 and Table 2 has also not indicated an appropriate correlation with numerical results (Fig. 5). It is interesting to note that, the Eq. 9 has lesser error than Eq. 31 with reference to numerical results. The Eq. 9 has used in two forms, original and modified. In modified form, like laminar high velocity state, re is replaced by rT, but it did not diminish the error.
The complicated shape of flow in turbulent condition and several assumptions imposed to achieve an analytical solution could be considered as the main reasons for these discrepancies between analytical and numerical results. Also, it must be noted that turbulence condition adds more parameters to numerical solution, which causes to insert at least 4 differential equations in the simplest form, these equations are not considered in analytical-empirical solutions.
Another result of numerical solutions is about the changes of hydrostatic pressure along the flow line in turbulent condition. According to Eq. 31 the hydrostatic pressure must be dropped along the general direction of flow toward outlet. But, as demonstrated in Fig. 6 the numerical analyses have indicated that at first, pressure increases to a maximum level and then it reduces as expected. This phenomenon is caused by an initial increase in turbulence velocity and consequently total velocity in the inlet that intensifies the dynamic pressure and diminishes the static pressure. But, at the middle part of the disc, it seems the turbulent velocity disappears and accordingly the falling trend of the hydrostatic pressure retains like laminar flow.
As the velocity increases, a discrepancy between numerical and analytical results is also stepped up. It can be considered as a significant reason to choose the lowest flowrate and pressure value among several (P, Q) pairs of a typical water pressure test for calculating the Lu value and consequently the hydraulic characterization.
| Table 4: | Obtaining some hydraulic characteristics with respect to Lugeon test |
| Fig. 6: | The changes of hydrostatic pressure along radial flow path toward outlet in turbulent condition. As expected, analytical Eq. 9, PA and also empirical Eq. 31, PCTR, in NS flow state show a falling trend but numerical results show initial increases of pressure, Q = 35 L min-1, Re = 5280, (a) b = 0.001 m and (b) 0.002 m |
Ewert (1997) has described a Lu value that calculated based on lowest (P, Q) point as the true permeability of rock mass, i.e., without being influenced by flow effects or changing in geometry of conductive spaces in rock joints. Actually, to achieve the best estimation of hydraulic aperture and other hydraulic characteristics, at least one of the WPT steps must be set in a way that the laminar and low inertia conditions would be considered. A summarized procedure is presented in Table 4 to show how the hydraulic characteristics of a borehole section would be determined by use of Lugeon test.
CONCLUSIONS
According to numerical and analytical calculations presented and discussed above, the following significant conclusions can be drawn.
As the Lugeon test is a single well test, the accuracy of hydraulic parameters estimated by use of this test is completely in doubt. For example, in a simple situation including a horizontal rough joint, as regards Fig. 3, merely a combination of three parameters, i.e., hydraulic equivalent aperture, relative roughness and effective radius, would be determined around borehole section test. In a typical tested section, some particular features such as joint network, joint filling material, contact areas in joints and probability of erosion, clogging, swelling or hydraulic fracturing are existed. Also, the flow is due to turbulence eddies and pressure losses caused by test equipments and geometry of rock conduits. According to those complicated situation, an applicable values describing aperture or roughness could not be identified. So, it is better to define some corresponding parameter such as equivalent hydraulic aperture or relative roughness to describe a general condition governing over rock nearby borehole.
There are several equations introduced here that are set empirically or analytically to portray the radial flow characteristics. Where the numerical method used as a verifier for these equations, the results show that despite of different assumption of analytical Eq. 9, only this simple equation generally coincide with the numerical results. In laminar low velocity conditions, this equation is perfectly matched to numerical results. In laminar condition when the inertia term is considerable, this equation could be corrected by a substitution Darcian radius, rD, for well radius, rw. Additionally, it also shows most well-matched result to numerical ones in comparison with other analytical relations extracted for turbulence conditions.
The two parameters rD and rT are also verified using numerical calculations. The results approved that critical relative intensity of inertia (ηc = 0.5) are choose appropriately because after correction of initial radius in Darcy equation the error decreases significantly. For verifying rT, The hydrostatic pressure is plotted along flow direction for turbulence conditions as shown in Fig. 6. If the maximum point of pressure curve is accepted as the point in which the turbulence effects disappear adequately after that (because of lowering the turbulence velocity), the turbulence radius did not match well enough to that point. It indicates that the turbulence conditions may differ and be complicated in radial flow in comparison with pipe or rectangular duct flow. Actually, the fluid has a potential of generating large eddies in horizontal plane at least near to inlet. But in the pipe or duct flow the fluid is limited to closer walls. So, the Reynolds number may not describe complete condition of turbulence situation if be calculated based on experimental Eq. 23 and 24 in radial flow.