An in situ thermo–hydraulic experiment in a saturated granite II: analysis and parameter estimation

ARTICLE IN PRESS

International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1395–1411 www.elsevier.com/locate/ijrmms

An in situ thermo–hydraulic experiment in a saturated granite II: analysis and parameter estimation E. Detournaya,, T. Senjuntichaib, I. Berchenkoc a

Department of Civil Engineering, University of Minnesota, 500 Pillsbury Drive S.E, Minneapolis, MN 55455, USA b Chulalongkorn University, Bangkok, Thailand c Shell Exploration and Production, The Netherlands Accepted 9 September 2004

Abstract This paper is Part II of a series of two papers concerned with an in situ Thermo–Hydraulic Experiment carried out at the Underground Research Laboratory of Atomic Energy of Canada Limited. The Thermo–Hydraulic Experiment was designed to determine the thermoporoelastic parameters that control the magnitude of the pore pressure induced by thermal loading. The experimental set-up involved a heater installed in a sub-horizontal borehole drilled from an underground gallery, and piezometers and thermistors located at different distances from the heater in auxiliary boreholes drilled from an adjacent gallery. Several water injection and heater tests were conducted during this experiment. Part II focuses on the interpretation of the experimental data to estimate some thermoporoelastic constants of the Lac du Bonnet Granite. The interpretation is based on matching the pore pressure and temperature data with theoretical predictions based on singular thermoporoelastic solutions. To estimate the parameters, the distances from the center of the test zone to all piezometers and thermistors have to be known. However, sufficiently accurate distances for back-analyses of the experiment could not be obtained due to measurement uncertainties. A two-step scheme is described in this paper to minimize the uncertainty in the distances. In a first step, parameters lumping distance and material constants were estimated by matching a theoretical prediction to experimental data. In a second step, the material constants were extracted from the lumped coefficients using a least square technique involving the tentative distances deduced from the borehole survey. r 2004 Elsevier Ltd. All rights reserved.

1. Introduction This paper is Part II of a series of two papers concerned with an in situ Thermo–Hydraulic Experiment (abbreviated to THE) carried out at the Underground Research Laboratory (URL) of Atomic Energy of Canada Limited. This field experiment, which included a series of water injection and heater tests, is described in Part I [1]. The main objective of the THE was to determine the values of thermoporoelastic parameters through back-analyses of the results from the in situ experiment. Corresponding author. Tel.: 612 625 3043; fax: 612 626 7750.

E-mail address: [email protected] (E. Detournay). 1365-1609/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijrmms.2004.09.007

Interpretation of THE to determine in situ values of the hydro–thermal coupling parameter g as well as the permeability k; the thermal diffusivity cn and the hydraulic diffusivity c of the Lac du Bonnet granite is based on matching the pore pressure and temperature data with the theoretical predictions based on singular thermoporoelastic solutions. Various types of source model are employed in an attempt to account for the geometry of the test zone. Due to some uncertainty involving the distances from the center of the test zone to the piezometer and thermistor sensors, a two-step scheme is developed to estimate parameters from THE. First, parameters lumping distance and material constants are estimated by matching theoretical predictions from selected source models to the experimental data.

ARTICLE IN PRESS 1396

E. Detournay et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1395–1411

Second, the thermoporoelastic constants are extracted from the lumped coefficients using a least square technique involving the tentative distances deduced from the borehole survey. Two distinct approaches are considered in the estimation of the lumped parameters in the first step. For the first method, the parameters are determined by matching only the time and the amplitude of the peak between the experimental responses and theoretical predictions. For the second one, the parameters are estimated by employing a nonlinear least square method [2,3] to match the full history of experimental data with predictions from the source solutions.

2. Thermoporoelasticity Our analysis of the THE experiment is carried out within the framework of the linear theory of thermoporoelasticity. Our favored formulation of this theory [4,5] is an extension of the classic Biot theory [6,7]. In this formulation, the porous medium is viewed as an open system [5]; i.e. an element of the porous material, identified by an element of the solid skeleton, can gain or lose fluid mass as it deforms. The essential thermo–hydraulic–mechanical (T–H–M) coupled effects are embodied in the volumetric response equations (the deviatoric response of an isotropic material is purely elastic) which can be written as s ¼ K u e  bMz  bu K u y;

p ¼ Mðz  be þ bm yÞ:

kr2 p ¼

1 @p @e @y þ b  bm ; M @t @t @t

cn r 2 y ¼

(2) @y : @t

@p @y g ; (4) @t @t where c ¼ k=S is the hydraulic diffusivity and g ¼ b0 =S is a hydro–thermal coupling coefficient. Both the storage coefficient S and b0 (a coupling coefficient of the same nature as bm ) are defined under oedometric condition (uniaxial strain) cr2 p ¼

S¼

1 3b2 þ ; M 3K þ 4G

b0 ¼ bm 

3bbK : 3K þ 4G

(5)

The thermal and hydraulic diffusivities, cn ¼ kn =C and c, the mobility k; and the hydro–thermal coupling coefficient g are the main parameters controlling the temperature and the pore pressure induced by heating of the rock.

3. Methodology of data interpretation

(1)

In the above equations, s ¼ skk =3 is the mean stress, p is the pore pressure, e is the volumetric strain, z is the variation of fluid content [6], and y is the temperature. The fields s; p, and y are variations with respect to initial equilibrated values and all the fields are initially zero. The constitutive equations can be combined with the balance laws (momentum, mass, energy) and conduction laws (Darcy and Fourier) to yield field equations for the displacement ui ; pore pressure p, and temperature y: For quasi-static processes of concern here, these equations can be expressed as 1 Gr2 ui þ ðG þ 3KÞe;i ¼ bp;i þ bKy;i ; 3

coefficient bu ; hydro–thermal coupling coefficient bm ); (v) transport (mobility coefficient k). Not all of these coefficients are independent, however, as the following relationships can be established between them: K u ¼ K þ b2 M; bM ¼ BK u ; bu ¼ bK=K u þ Bbm : A linear isotropic rock is thus characterized by nine independent thermoporoelastic parameters. For the particular case of an irrotational displacement field induced by thermal loading in an infinite domain, the pore pressure diffusion equation (3) reduces to [1]

(3)

The material constants appearing in these equations can be conveniently organized into five groups: (i) mechanical (shear modulus G, drained bulk modulus K); (ii) thermal (drained thermal expansivity b; volumetric heat capacity C, thermal conductivity kn ); (iii) hydromechanical (Biot coefficient b, Biot modulus M, Skempton coefficient B, undrained bulk modulus K u ); (iv) thermohydromechanical (undrained thermal

We now turn our attention to the determination of several thermoporoelastic parameters of the Lac du Bonnet granite (c, k; cn ; L0 =kn ; g) from data collected during the in situ THE conducted at the Underground Research Laboratory. This experiment, which is described in the companion paper [1], consisted of an injection test and three heater tests. The approach used to estimate the thermo-hydraulic parameters makes use of closed form thermoporoelastic heat and fluid sources and proceeds in two steps. First, it relies on the determination of two lumped quantities at each sensor (and for each test), by matching the theoretical response with data and in a second step the material parameters are deduced from the lumped quantities using the information about the distance between the source and the measurement points. As explained below, this approach attempts to minimize the consequence of the inaccuracies in the exact location of the pore pressure and temperature sensors. Various types of source solutions are actually used to interpret the field data. For the injection test, we make use of fluid sources characterized by an instantaneous injection of a (known) finite volume of fluid O0 ; in these source solutions, the rate of fluid volume injected varies according to a Dirac delta function. Both a point source

ARTICLE IN PRESS E. Detournay et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1395–1411

and line source are used to model the injection test; the line source with a length L corresponding to the length of the test zone is an attempt to crudely account for the geometry of the test zone. Interpretation of the heater tests relies on thermal point source and line source to model the temperature, and a combination of a step heat source and a time-varying fluid sink point source and a spherical source to interpret the thermally induced pore pressure. All the relevant source functions are listed in Appendix A. Those source functions are functions of R or r and z depending on the types of model we consider. These variables represent the distances from the center of the test zone to the piezometer and thermistor sensors and play an important role in the interpretation of the field data. As discussed in the companion paper [1], some uncertainty in the location of the center of the test zone and the sensors meant that the distances compiled from the borehole survey can only be considered as tentative values. Since uncertainty in these distances can result in a significant error in the estimated values of parameters, the back-analysis should involve a procedure to improve their estimation. Consideration of these distances as parameters to be estimated together with the material parameters has a few shortcomings. First, it involves a lot more parameters to be estimated, e.g. 12 more parameters (two unknown distances, r and z; for six thermistors) would have to be considered in the interpretation of the temperature data with the heater modeled as a line source. More parameters usually lead to more uncertainty in the values of the estimated parameters. In addition, in some cases, a unique estimate of the distances cannot be obtained together with the material parameters. For example, in the case of a point source model, the distances normally appear with the material parameters; as a result, the distances and the material parameters cannot be simultaneously and uniquely determined. A two-step scheme has then been devised to extract the material parameters from the pore pressure and temperature data at the various sensors, in an attempt to minimize the effect of the uncertainty in the distances between the sensors and the test zone. In a first step, parameters lumping distance and material constants are estimated by matching a theoretical response to an experimental one at a particular sensor. In fact, each measured quantity m (pore pressure or temperature) can generically be expressed using a source solution as m ¼ Pm1 f ðt; P2 Þ;

(6)

where f is a known function, Pi ¼ ani RZ ; i ¼ 1; 2 is a lumped parameter that depends on another constant ai (usually a material constant) and on the distance R between the source (or the center of the source) and the

1397

point of measurement (i.e. the sensor), and m; n; and Z are power law exponents. (Actually, m and n only take the value 1 or 1/2 and Z is simply 1 or 1:) By matching the theoretical response (6) with the measured response fmðkÞ j ; tj g at a particular sensor k, the two lumped parameters PðkÞ i ; i ¼ 1; 2 can be estimated for this sensor. Two methods are actually employed to determine the two lumped coefficients from the temperature and pore pressure data. In the first method, the lumped coefficients are estimated by matching the time and amplitude of the peak of the pore pressure and temperature responses with theoretical predictions. The second method is based on matching the full history of the experimental data with a theoretical response. The concept of nonlinear parameter estimation involving a least square method is employed to find the unknown parameters giving the least error in the objective function. The objective function is the summation of the squared difference between predicted and measured quantities. In the second step, the material constants are extracted from the lumped coefficients using a least square technique involving the tentative distances.

4. Interpretation of THE data using source solutions 4.1. Pore pressure in injection test The solutions of an instantaneous fluid source (point and line) are used to interpret the pore pressure induced by injection of a volume O of water in the rock mass. As discussed above, the pore pressure measurements are first interpreted in terms of two lump quantities U1 and U2 which depends on the distance R between the source (or the center of the source) and the point where the pore pressure is measured. These two lumped quantities are defined as pffiffiffi c and U2 ¼ kR: U1 ¼ (7) R Consider first the pore pressure induced by an instantaneous point source (see Appendix A). In terms of the lumped quantities U1 and U2 ; the pore pressure change can be expressed as ! O0 1 pffiffiffiffiffiffiffiffiffi exp  2 : p¼ (8) 4U1 t 8U1 U2 p3 t3 The pore pressure is characterized by a peak amplitude pp reached at time tp : The parameters U1 and U2 can be determined from the knowledge of pp and tp as sffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffi 1 O0 6 U1 ¼ ; U2 ¼ : (9) 6tp 8pp tp p3 e3

ARTICLE IN PRESS E. Detournay et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1395–1411

Consider next the line source solution. The pore pressure equation can now be rewritten as !  "  O0 cos2 y sin y þ L pffiffiffiffiffi exp  p¼ Erf 16pU2 Lt 4U12 t U1 4t   # sin y  L pffiffiffiffiffi  Erf ; ð10Þ U1 4t where Erf is the error function and

(12) where Z Kðx; ZÞ ¼ pffiffiffi þ ð1  x2 ÞErfðZÞ: p expðZ2 Þ Thereafter, U2 can be deduced from (10) by setting t ¼ tp and p ¼ pp : Alternatively, U1 and U2 can be estimated from curve matching using either the point source or the line source solution, see Appendix B. Determination of U1 and U2 using either the time and amplitude of the pore pressure peak or the matching of the full response is summarized in Appendix C. Figs. 1 and 2 show a comparison between experimental data and theoretical solutions using the point source and line source models, respectively. A very good agreement between experimental and theoretical responses is noted in both figures, especially for PZ6 for

0.08 ED PM LS

p (MPa)

PZ1

0.06

0.04

PZ6 0.00 0

0.00 5

10

15

10

15

20

20

25

30

35

t (day)

Fig. 1. Injection test: pore pressure at PZ1 and PZ6 (point source model).

25

30

35

Fig. 2. Injection test: pore pressure at PZ1 and PZ6 (line source model).

both peak and curve matching. The difference between the predicted and observed responses using the point source model shown in Fig. 1 for PZ1 is most likely associated with the influence of the actual geometry of the test zone, as PZ1 is the closest piezometer to the main borehole. 4.2. Temperature induced by heating The temperature data at the various thermistors can be compared to the temperature field corresponding to a continuous point source and a continuous line source to estimate the following lumped quantities X1 and X2 defined as pffiffiffiffiffi cn L0 and X2 ¼ ; (13) X1 ¼ R kn R where R is the distance between a field point and the center of the source and L0 is the heat rate applied to the rock. As discussed in the companion paper [1], L0 is less than the power of the heater due mainly to losses associated with heat conduction through the borehole casing and the packer system. Consider first the point heat source (see Appendix A). The temperature response can be rewritten in terms of X1 and X2 as X2 ½Os ðt; X1 Þ  Os ðt  td ; X1 Þ; 4p where   1 pffiffiffiffiffi : Os ðt; X1 Þ ¼ Erfc X1 4t

PZ6

0

5

t (day)

y¼

0.04

0.02

ED PM LS

PZ1

0.02

r z L (11) cos y ¼ ; sin y ¼ ; L ¼ : R R R Note that y and L are computed from the tentative values of r; z and R and the half-length of the test zone L. The lumped coefficient U1 can be estimated from the peak time tp using ! ! cos y sin y þ L cos y sin y  L pffiffiffiffiffiffi ; pffiffiffiffiffiffi  K pffiffiffiffiffiffi ; pffiffiffiffiffiffi ¼ 0; K U1 4tp U1 4tp U1 4tp U1 4tp

0.06

0.08

p(MPa)

1398

(14)

The temperature response is characterized by a peak amplitude yp at a time tp which is larger than the duration of heating td : The two quantities X1 and X2 can thus be estimated from tp and yp ; respectively. By setting the time derivative of (14) to zero, X1 can be expressed in

ARTICLE IN PRESS E. Detournay et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1395–1411

terms of tp and td as

4 PZ1-T

where y and L are given by (11). Then, by the same technique, X1 can be determined from (17)

θ (oC)

3

Once X1 has been estimated, the parameter X2 can be computed from (14) using the peak amplitude yp and tp : It can be shown, however that the ratio tp =td tends towards 1 with increasing td : Thus for ‘‘large’’ td ; any uncertainty in the time of the peak tp from the temperature response will result in a significant error in the estimation of X1 : The time of the peak tp is distinct enough from td at the various thermistors in the first and second heater test, but not in the third test where td 50 days. Thus the technique of peak matching to estimate X1 and X2 is only used for the first two heater tests. For all three tests, X1 and X2 can also be estimated by curve matching, as discussed in Appendix B. Consider next the line source model. The induced temperature is also given by (14) with Os replaced by O defined as Z L 1 1 OðtÞ ¼ s qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2L L ðsin y  xÞ2 þ cos2 y sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! ðsin y  xÞ2 þ cos2 y

Erfc dx; ð16Þ 4X21 t

F ðtp Þ  F ðtp  td Þ ¼ 0;

2

1

PZ6-T

0 0



 cos y sin y þ L pffiffi ; pffiffi ; t 2X1 t 2X1 t   cos y sin y  L pffiffi ; pffiffi ; t G 2X1 t 2X1 t

F ðtÞ ¼ G

with ErfðyÞ : t expðx2 Þ

Thereafter, the parameter X2 can be determined from (14), with Os ðtÞ replaced by OðtÞ; using yp and tp : Alternatively, X1 and X2 can be estimated by curve matching. Details of the estimation of X1 and X2 can be found in Appendix C. Figs. 3 and 4 present, for the second heater test, a comparison between the measured and predicted temperature at PZ1-T and PZ6-T, based on the point and line source models, respectively. It can be seen from these figures that a very good agreement between experimental and theoretical responses are observed for both peak matching and least square fitting. As expected, a better fit in the figures is once again found at the farthest thermistor (PZ6-T).

3

6

9

t (day)

Fig. 3. Second heater test: temperature at PZ1-T and PZ6-T (point source model).

4 ED PM LS

PZ1-T 3

2

PZ6-T

1

0 0

where

Gðx; y; tÞ ¼

ED PM LS

(15)

θ (oC)

1 X1 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 6tp ðtp =td  1Þ½ln tp =ðtp  td Þ

1399

3

6

9

t (day)

Fig. 4. Second heater test: temperature at PZ1-T and PZ6-T (line source model).

Comparison between experimental data and theoretical predictions of the temperature at PZ1-T and PZ6-T is shown in Fig. 5 for the third heater test for both point and line source models using a nonlinear least square fit. Note the excellent agreement between the measured and computed responses for both PZ1-T and PZ6-T. 4.3. Pore pressure induced by heating The pore pressure field induced by heating can be modeled using either the combined solution of a fluid and heat point source, or the solution of a spherical cavity with a constant heat flux and constant pore pressure imposed at the boundary (see Appendix A). The lumped quantities to be determined from the pore pressure response at a piezometer (either by matching or by using the peak) are G1 and G2 defined as g a (18) G1 ¼ ; G2 ¼ ; R R

ARTICLE IN PRESS E. Detournay et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1395–1411

12

θ (oC)

ED PS LS

PZ1-T

8 PZ6-T

4

0 0

25

50

75

100

t (day)

Fig. 5. Third heater test: temperature at PZ1-T and PZ6-T.

where a denotes the effective drainage radius. This parameter is the distance from the center of the source at which the pore pressure is assumed to be held constant. In the experiment this occurred at the borehole boundary. However, for an effective spherical cavity the value of a was considered an unknown parameter determined through back-analysis. Consider first the combined heat and fluid source singularity. The pore pressure can be rewritten as p¼

L0 G1 ½P1 ðt; G2 Þ  P1 ðt  td ; G2 Þ; 4pkn ð1  o2 Þ

where



"  2 # 1 1 pffiffiffiffiffiffiffi exp  pffiffi MðtÞ ¼ 2X1 t 2X1 pt3   1  G2 G2 pffiffiffiffiffiffiffi þ pffiffiffiffiffiffiffi  2oX1 pt3 2X1 pt3 "  2 # 1  G2 G2 pffiffi þ pffiffi

exp  : 2oX1 t 2X1 t

 1 pffiffi P1 ðt; G2 Þ ¼ HðtÞ Erfc 2X1 t   1  G2 G2 pffiffi þ pffiffi Erfc 2oX1 t 2X1 t

and X1 is the thermal lumped parameter defined by (13). The pore pressure p at a fixed distance R is once again characterized by a peak amplitude pp at time tp : Depending on the duration of the thermal load td and the distance R, the time of the peak tp could be smaller or larger than td : If tp otd the peak is said to be natural and both tp and pp are obviously independent of td : (In heater test #3, the pore pressure peaks naturally at piezometers PZ1, PZ3, and PZ5.) The peak time tp is thus given explicitly by   ð1  G2 þ oG2 Þ2  o2 1 ; tp otd : tp ¼ (20) 4o2 X21 lnð1  G2 þ oG2 Þ  ln o The above equation can be solved for G2 using the peak time tp and assuming that o and X1 have been determined before. If tp 4td (as in the heater tests #1 and #2), G2 is deduced from the time of the peak by solving (21)

ð22Þ

The lumped parameter G1 is then computed from (19) identified at the pore pressure peak. Alternatively, G1 and G2 can be estimated by matching the theoretical pore pressure response (19) with the experimental data. For the spherical cavity model, the change in pore pressure is also given by (19) but with P1 ðt; G2 Þ being replaced by P2 ðt; G2 Þ 

  1  G2 ; G2 ; P2 ðt; G2 Þ ¼ HðtÞ Cðt; 1  G2 ; G2 Þ  C t; o (23) where    x x X21 t pffiffi  exp þ 2 G2 G2 2X1 t pffiffi  x X1 t pffiffi þ

Erfc : G2 2X1 t

 Cðt; x; G2 Þ ¼ Erfc

(19)



Mðtp Þ  Mðtp  td Þ ¼ 0;

where

ð24Þ

Comparison between measured and predicted temperature at PZ1 and PZ6 for the second and third heater test are shown in Figs. 6–9. A very good agreement between the experimental and predicted responses is clearly observed in these figures. Note that the theoretical predictions obtained from the spherical source model yield a better fit than those from the point source one.

ED PM LS

PZ1 2

p (MPa)

1400

1 PZ6

0

0

3

6

9

12

15

t (day)

Fig. 6. Second heater test: pore pressure at PZ1 and PZ6 (point source model).

ARTICLE IN PRESS E. Detournay et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1395–1411

1401

5. Parameter estimation PZ1

ED PM LS

p (MPa)

2

1 PZ6

0

0

3

6

9

12

15

t(day)

Fig. 7. Second heater test: pore pressure at PZ1 and PZ6 (spherical source model).

3 PZ1

ED PM LS

2

p (MPa)

PZ6 1

0

-1

-2 0

25

50

75

100

t (day)

Fig. 8. Third heater test: pore pressure at PZ1 and PZ6 (point source model).

3 PZ1

ED PM LS

2 PZ6 p (MPa)

1

0

-1

-2 0

25

50

75

100

t (day)

Fig. 9. Third heater test: pore pressure at PZ1 and PZ6 (spherical source model).

The hydraulic diffusivity c, the mobility coefficient k; the thermal diffusivity cn the hydro–thermal coupling coefficient g; as well as L0 =kn and the effective drainage radius a can finally be assessed using a least square fit involving the lumped parameters U1 ;U2 ðc; k), X1 ; X2 (cn ; L0 =kn ) and G1 ; G2 ðg; aÞ and the tentative distances between the sensors and the test zone based on field measurements. The detailed results for all tests are presented in Appendix C. Hydraulic diffusivity and mobility. Interpretation of the injection test data yields an estimated hydraulic diffusivity c varying between 4:6 108 and 1:8 107 m2 =s depending on the model used and the method of parameter estimation and the estimated mobility k varying in the range 6 1015 –9:3 1015 m2 =kPa s (see Table 4). The estimated distances for R1 and R6 obtained from the point source model seem to be reasonable when compared to the distances obtained from the borehole survey [1]. However, there is a discrepancy between the distances estimated with the line source model and those determined from the borehole survey. Thermal diffusivity. The estimated values of cn and L0 =kn as well as the computed value of the distances between the center of the test zone and six thermistors can be found in Tables 8–10 for the first, second and third heater tests, respectively. Estimate of cn ranges from 1:0 106 to 1:9 106 m2 =s; while L0 =kn is estimated to be in the interval between 87 and 105 1C m. There is therefore a very good agreement between the various values of cn and L0 =kn obtained with different types of matching and source models. Using kn ’ 3:5 W=m C as estimated from other in situ experiments at the URL, these values represent a heater efficiency of between 76% and 90%. This is in the range of what would be expected due to heat through the borehole as described in the companion paper [1]. Moreover, unlike the injection test, both point and line source models yield good estimation of the distances R, when compared to the tentative values compiled in the same tables of Appendix C. Finally, we note that the estimated parameters given in Table 10 for the third test, which has the longest duration of thermal load td ; show more consistent results than those in the other two tests. Hydrothermal coupling parameter g: The hydrothermal parameter g is estimated to be in the range 425–475 kPa/1C and the effective drainage radius a in the range 0.09–0.34 m (see Tables 14–16 in Appendix C for details using the point and the spherical source models). Although the three heater tests give coherent values for g and a, more consistent results are found in Table 16 for the third heater test, in which the duration of heating td is longer than those of the other two tests. The estimated values of a are reasonable since they are

ARTICLE IN PRESS E. Detournay et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1395–1411

1402

Table 1 Estimated values of thermoporoelastic parameters of Lac du Bonnet granite Types of source solutions

c (m2/s) ð 107 Þ

k (m2/kPa-s) ð 1015 Þ

cn (m2/s) ð 106 Þ

g (kPa/1C) ð 102 Þ

Point Line Spherical

1.66–1.77 0:46  0:49 —

8:20  9:25 6:05  6:34 —

1:32  1:92 1:03  1:22 —

4:36  4:75 — 4:26  4:51

usually more than the radius of the main borehole (0.096 m) and less than the half-length of the test zone (0.397 m). The main results of the parameter estimation are summarized in Table 1.

1

log10 θ

ED PM/P LS/P PM/L LS/L

0

6. Coherence of estimated parameters 6.1. Asymptotic analysis of temperature response

does not depend on distance from the source. In a log/log plot of temperature with respect to time, the temperature response from all thermistors should converge towards ys ðtÞ: For time t much larger than the duration of heating td ; ys ðtÞ will become a straight line with a slope of 1=2: Figs. 10 and 11 shows the temperature response from all six thermistors in the second and the third heater tests, respectively, together with the asymptotic temperature yas ðtÞ: The asymptotic temperature is obtained from the estimated parameters given in Table 9 for the second heater test and Table 10 for the third heater test, respectively. It can be seen from Figs. 10 and 11 that the ‘‘large’’ time has not been reached yet. Unfortunately, the temperature response measured after time t ¼ 106 s was too small and noisy. Nevertheless, it can be concluded that predicted asymptotic behavior of the temperature is consistent with the experimental data. Again, the line source model gives the best approximation. 6.2. Analysis of natural peak of pore pressure response When the duration of the thermal load td in the heater test is large enough, the time of the peak of pore pressure response tp is independent of td : The peak is then called a ‘‘natural’’ peak and the time of the peak in this case according to a point source model is given by tp ¼

½R22  R21  ; ln½R2 =R1 

(26)

-1 5.1

5.4

5.7

6.0

log10 t

Fig. 10. Comparison between all temperature curves from the second heater test and the asymptotic temperature evaluated using estimated parameters from Table 9.

1.0 ED LS/P LS/L

log10 θ

Consider the temperature response given by (A.3) in Appendix A. The asymptotic temperature ys ðtÞ   L0 1 1 ys ¼ (25) pffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi  pffiffi t  td 4pkn pcn t

0.5

0.0 6.65

6.70

6.75 log10 t

6.80

6.85

Fig. 11. Comparison between all temperature curves from the third heater test and the asymptotic temperature evaluated using estimated parameters from Table 10.

where a2 R R 1 þ1 (27) tn ¼ ; R1 ¼ ; R2 ¼ 1 a o 4cn pffiffiffiffiffiffiffiffiffi and o ¼ c=cn ; the square root of the ratio of the hydraulic to thermal diffusivity. It can be seen from the above equation that tp depends on four parameters,

tp ¼

tp ; tn

ARTICLE IN PRESS E. Detournay et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1395–1411 Table 2 Estimated parameters calculated from the natural peak at PZ1, PZ5 and PZ6 in the third heater test

0.30 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38

R1 PZ1

PZ5

PZ6

2.91 3.21 3.58 4.07 4.73 5.70 7.21 9.92 16.23

2.87 3.16 3.52 3.98 4.62 5.52 6.92 9.35 14.71

2.10 2.22 2.35 2.51 2.70 2.93 3.20 3.55 3.99

300 Amount of water (ml)

o

1403

200

100 ED PS SS

0 0

10

20

30

40

50

t (day)

Fig. 12. Accumulated amount of water entering the test zone in the third heater test.

the distance R, the effective radius a and the two diffusivities c and cn : Consider the third heater test in which td is long enough that the natural peak of pore pressure response occurred in PZ1, PZ5 and PZ6. By using the value of cn ¼ 1:32 106 m2 =s obtained from the interpretation of the temperature field and tp observed at those three piezometers, the parameter R1 can be obtained from (20) for different values of ratio of diffusivities o: Table 2 shows the estimated values of R1 using the time of the natural peak tp and (20) for different values of o: From the values of R1 ; the suitable values of o can be obtained. Since the appropriate value of a should be between 0.096 m (the diameter of the borehole) and 0.397 m (half-length of the test zone) and the expected distance R for PZ1, PZ5 and PZ6 computed from the survey are 0.48, 1.82 and 1.07 m, respectively, the square root of the diffusivity ratio can be estimated as o 0:33–0.35. This is consistent with the value of o obtained from the back-analysis, which is about 0.35–0.37 ðc ¼ 1:66  1:77 106 m2 =s from the injection test with a point source model). 6.3. Interpretation of amount of water entering the test zone In all three heater tests, heat is produced at a constant rate over a time td ; while keeping pore pressure in the sealed-off interval (the test zone) at a constant value. The amount of water flowing in the test zone is then measured during testing until the heat power is turned off. The accumulated amount of water can actually be calculated according to the point source model by integrating the fluid source strength whereas, for the spherical source model, it is obtained by integrating the water flux across the boundary. Let V w ðtÞ denote the accumulated amount of water entering the test zone

at time t. Then, rffiffiffi        $ t $2 exp  V w ðtÞ ¼ V n t þ 2$2 Erf pffiffiffi þ 2$ p t t (28) for the point source model, and    rffiffiffi rffiffiffi t t t V w ðtÞ ¼ 4$V n exp Erfc þ 4 4 p

(29)

for the spherical source model. In the above equations, Vn ¼

L0 gktn ; kn ð1  o2 Þ

$¼

o1 ; o

t¼

t ; tn

tn ¼

a2 : 4cn (30)

Fig. 12 presents the comparison of the accumulated amount of water entering the test zone in the third heater test between measured data and theoretical predictions by both point and spherical source models using (28) and (29), respectively. The parameters employed in the predictions by both equations are obtained from Table 10 (cn ¼ 1:32 106 m2 =s and L0 =kn ¼ 105 C m) and Table 16 (least square curve matching). It can be seen once again that a good agreement between experimental and theoretical responses is obtained.

7. Conclusions This paper has focused on the estimation of the hydro-thermal coupling parameter g as well as the permeability k; the thermal diffusivity cn and the hydraulic diffusivity c of the Lac du Bonnet granite using data from the in situ THE conducted at the Underground Research Laboratory (URL) in Manitoba, Canada, which was described in the companion paper [1]. A two-step scheme has been developed in this

ARTICLE IN PRESS 1404

E. Detournay et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1395–1411

paper to minimize the effect of uncertainty on the distances between the center of the test zone and those sensors. In a first step, parameters lumping distance and material constants are estimated by matching a theoretical prediction based on source solutions to experimental data. In a second step, the material constants are extracted from the lumped coefficients using a least square technique involving the tentative distances deduced from the borehole survey. The in situ values of thermoporoelastic constants obtained from the back-analysis of all tests are summarized in Table 1. A good agreement of the estimated values between different types of models except for the value of the hydraulic diffusivity c obtained from the line source model which is lower than the other estimates. In fact the hydraulic diffusivity estimated from the THE is about one order of magnitude smaller than the tentative value used for designing the experiment (see Table 1 of Part I [1]). Interpretation of the data indicates that the ratio of the thermal to hydraulic diffusivity cn =c ’ 10 for the LdB granite. This large ratio implies that significant thermohydraulic coupling is to be expected. Moreover, the coupling parameter g; which controls the maximum thermally induced pore pressure, is estimated to be between 0.4 and 0.5 MPa/1C. In spite of the uncertainty on the distances from the center of the test zone to all sensors, the estimated values obtained from the methodology presented in this paper indicate that very consistent results are obtained from the interpretation of the field data for both peak matching and least square curve matching techniques. In addition, the results are also verified from the redundancy and coherence of the test results, which include the asymptotic analysis, the interpretation of the natural peak of pore pressure response in the third heater test and the interpretation of the data of the amount of water entering the test zone in the third heater test.

Acknowledgements Support for the research reported in this paper has been provided by a research contract funded by Ontario Power Generation as part of the Deep Geologic Repository Technology Program. This support is gratefully acknowledged. The authors would like to express their deep thanks to Dr. Neil Chandler for his steadfast support to this research over several years, his guidance, and his careful review of the manuscript.

thermoporoelasticity. An instantaneous fluid source is used to model the injection test, whereas, in the heater test, a combination of a step heat source and a timevarying fluid sink is employed to enforce the hydraulic boundary conditions. A.1. Injection test A.1.1. Point source model Consider the injection test where a given volume of water O0 is pumped quasi-instantaneously in a sealedoff interval (test zone) of the borehole. Such a test can be modeled with an instantaneous point source of strength O0 located at the center of the test zone. The change in pore pressure in this case can be written as   O0 R2 pffiffiffiffiffiffiffiffiffi exp  pðR; tÞ ¼ ; (A.1) 4ct 8pk pct3 where R is the distance from the center of the test zone and t is time. A.1.2. Line source model To obtain a better representation of the geometry of the test zone, the injection test can be modeled with a line source. Let ps ðR; tÞ denote pore pressure induced by an instantaneous point source of a constant strength O0 ; i.e., ps ðR; tÞ is given by Eq. (A.1). Then, the pore pressure corresponding to the line source can be constructed by integrating a source singularity along the line segment of the test zone, i.e., jzjpL; and r ¼ 0: By assuming that the source density is uniformly distributed along the segment, the change in pore pressure can be written as      O0 r2 zþL pðr; z; tÞ ¼ exp  Erf pffiffiffiffiffiffiffi 16pkLt 4ct 4ct    zL  Erf pffiffiffiffiffiffiffi : ðA:2Þ 4ct A.2. Heater test A.2.1. Temperature field Point source model. The heater test can be modelled using a continuous point heat source of constant strength L0 located at the center of the test zone, acting over the time interval td : The change in temperature is thus given by y¼

L0 ½Fs ðR; tÞ  Fs ðR; t  td Þ; 4pkn

(A.3)

where Appendix A. Source solutions

  HðtÞ R Erfc pffiffiffiffiffiffiffiffiffi : Fs ðR; tÞ ¼ R 4cn t

Interpretation of the injection test and the heater test is based on different source solutions from the theory of

Line source model. Consider the line source model for the heater test. The change in temperature

ARTICLE IN PRESS E. Detournay et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1395–1411

is also given by (A.3) with Fs being replaced by F defined as Z L qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  1 Fðr; z; tÞ ¼ ðz  xÞ2 þ r2 ; t dx (A.4) Fs 2L L Note that the integral in (A.4) has to be evaluated numerically. A.2.2. Pore pressure field Point source model. An approximate solution of the pore pressure response in the heater test can be obtained by superposing a time-varying fluid sink onto the continuous heat source. A fluid sink has to be introduced to satisfy the drained boundary condition on the borehole walls. The strength of this sink, which represents the total fluid discharge entering the borehole, can be calculated from the condition that the pore pressure is imposed at some radial distance a from the source. If the test zone were a perfect sphere of radius R0 ; a would be equal to R0 : For any other geometry, this distance a represents an ‘‘effective’’ radius which is unknown. Therefore, a has to be determined together with material properties in the back-analysis of the experimental data. By using the thermoporoelastic source solution, the pore pressure field induced by the heater produced at a constant rate L0 during a certain period of time td is given by p¼

L0 g ½C1 ðR; tÞ  C1 ðR; t  td Þ; 4pkn ð1  o2 Þ

(A.5)

C1 ðR; tÞ ¼

where   hr c ti r n Cðr; tÞ ¼ Erfc pffiffiffiffiffiffiffiffiffi  exp þ 2 Erfc a a 4cn t  pffiffiffiffiffiffi cn t r

pffiffiffiffiffiffiffiffiffi þ : a 4cn t

The curve matching technique involves a least square method to find the unknown parameters giving the least error in the difference between predicted and measured quantities [3,2]. This technique has successfully been used in the past for the back-analysis of an in situ heater test [8]. The general principle of this approach is described below. Let a measured quantity y be related to a vector of independent variables x and an M-dimensional vector of unknown parameters t by means of y ¼ f ðx; tÞ:

with R r1 ¼ pffiffiffiffiffiffiffiffiffi 4cn t

and

Ra a r2 ¼ pffiffiffiffiffiffiffi þ pffiffiffiffiffiffiffiffiffi : 4cn t 4ct

Spherical source model. The pore pressure field from the heater test cannot be modeled using the line source solution since it is not possible to satisfy the hydraulic boundary condition ðp ¼ 0Þ along the whole segment of the test zone. This condition can only be satisfied by introducing a space- and time-dependent fluid source. The alternative to the combined heat and fluid source is a spherical heat source of an unknown ‘‘effective’’ radius a and a constant heat flux Q0 ¼ L0 =4pa2 ; with drained boundaries ðp ¼ 0Þ: Then, the change in pore pressure corresponding to this source is also given by (A.5) with C1 being replaced by C2 defined as    1 Ra CðR  a; tÞ  C ;t ; (A.6) C2 ðR; tÞ ¼ R o

(B.1)

Consider a series of N measurements. Then i ¼ yi  f ðxi ; tÞ;

(B.2)

where a subscript ‘‘i’’ denote the values associated with ith experiment and i is an experimental error. By minimizing a selected objective function CðtÞ (usually the least squares objective function) the unknown parameters t can be found t

HðtÞ ½Erfcðr1 Þ  Erfcðr2 Þ R

ðA:7Þ

Appendix B. Curve matching

t0 ¼ min T ;

where

1405

(B.3)

where  is N-dimensional vector of errors. The minimizing value of t; t0 ; is usually found by means of an iterative procedure tnew ¼ told þ Dtold :

(B.4)

Various gradient-based methods can be used to compute Dtold : They can all be written in the common form tnew ¼ told  rRg;

(B.5)

where r is a strictly positive scalar (known as a step length), R is an M-dimensional positive definite matrix used to modify the descent direction, and g is the M-dimensional gradient vector of CðtÞ; i.e.,  @CðtÞ g¼ (B.6) @t u¼uold There are a variety of gradient-based methods to be employed in the parameter estimation as discussed below. They are different in the definition of the matrix R:



Steepest descent. In this method, R ¼ I (the identity matrix). Although this method will always converge, it is often very inefficient in obtaining a solution.

ARTICLE IN PRESS 1406





E. Detournay et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1395–1411

Thus, it is not recommended for practical applications. Newton’s method. In this method r ¼ 1 and R ¼ H1 where H is the Hessian matrix defined as  @2 CðtÞ H¼ : (B.7) @t2 t¼told

source), LS/P (least square curve matching with point source), PM/L (peak matching with line source), LS/L (least square curve matching with line source), LS/P (least square curve matching with point source), PM/SP (peak matching with spherical source) and LS/SP (least square curve matching with spherical source).

Although Newton’s method is faster than the steepest descent method, there are two deficiencies in this method. First of all, it requires the evaluation of the second derivatives which may reduce its efficiency if the objective function is very complicated. In addition, there is no guarantee that the algorithm will always converge since the Hessian matrix H is not necessary positive definite except near the minimum. Gauss (or Gauss–Newton) method. By omitting the terms involving the second derivatives in the above equation, an approximation of the Hessian matrix N is used in place of H: Therefore, in this case, i.e. R ¼ N1 where N ¼ 2JT J and J is the Jacobian matrix defined as  @f ðx; tÞ J¼ : (B.8) @t 

C.1. Injection Test

t¼told



Marquardt (or Levenberg–Marquardt) method. In this method, R ¼ ðN þ lIÞ1 ;

(B.9)

where l is a scalar which controls the algorithm. When l ! 0; it approaches the Gauss method whereas, for l ! 1; the matrix in the bracket in the above equation is forced into being diagonally dominant, so the algorithm will proceed in the same direction as in the steepest descent method.

Two methods of parameter estimation are employed in the determination of the lumped coefficients U1 and U2 given by (7). Pore pressure data from two working piezometers, namely, PZ1 and PZ6 are used. For the peak response matching, the parameters are estimated using (9) for the point source model, and from the root of (12) together with (10) for the line source model. For the nonlinear least square method, the number of observations used in the parameter estimation are 690 from each piezometer. Note that all noticeable spikes, which can be seen in Fig. 8 of Part 1 [1], are removed from the experimental data before it is used in the estimation. Estimated values of the parameters U1 and U2 from the back-analysis of the injection test are shown in Table 3. It can be seen that the values of U1 and U2 obtained in the same model from both methods are very close in each piezometer. When compared with different models, they are still within an order of magnitude. Note the very consistent results obtained from PZ6. Once U1 and U2 have been estimated, the hydraulic diffusivity c and the mobility coefficient k are computed by performing a linear regression on the sets flog U1 ; log Rg and flog U2 ; log Rg: Once the values of c and k are determined, the values of the distance R can be computed from (7) in each case. The estimated values of c and k as well as the computed values of the distances R1 and R6 together with the tentative distances estimated from the URL Borehole Collar Survey Information Report are listed in Table 4.

Appendix C. Pore pressure and temperature matching C.2. Heater Test In this appendix, determination of the lumped parameters U1 ; U2 ; X1 ; X2 ; G1 and G2 are presented for the injection test and the three heater tests, together with the determination of c, k; cn ; g: Two types of parameter estimation, namely, peak matching and least square curve matching, are employed to match pore pressure and temperature data with predictions from various types of source solutions from theory of thermoporoelasticity. In the least square method, the minimization of the objective function is performed by using the function nlinfit provided in MATLAB. This function makes use of the Gauss–Newton algorithm with Levenberg–Marquardt modifications for global convergence. The following notation is used in the tables compiled in this appendix: PM/P (peak matching with point

C.2.1. Temperature Response The data from six thermistors (PZ1-T, PZ4-T, PZ6-T, PZ7-T, PZ8-T and PZ9-T) are used to analyze the temperature data to estimate the lumped parameters X1 and X2 given by (13). Both point source and line source models are considered for the three heater tests. Two methods of parameter estimation are employed in the back-analysis. The technique of matching the peak temperature response is not applicable for the third test since tp  td : Therefore, only the temperature data from the first and second heated tests were analysed using peak matching. For the nonlinear least square method, the number of observations used in the parameter estimation from each thermistor are 900, 1000 and 2800, respectively, for the first, second and third heater test. In

ARTICLE IN PRESS E. Detournay et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1395–1411

1407

Table 3 Injection test: estimated lumped parameters U1 and U2 Estimation method

PZ1

PM/P LS/P PM/L LS/L

PZ6

U1 pffiffiffiffiffiffiffi ð 1=sÞ

U2 (m3/kPa-s)

U1 pffiffiffiffiffiffiffi ð 1=sÞ

U2 (m3/kPa-s)

6:76 104 6:33 104 2:28 104 2:26 104

5:90 1015 4:63 1015 3:29 1015 3:19 1015

3:86 104 3:86 104 3:18 104 3:02 104

9:83 1015 9:87 1015 8:29 1015 7:81 1015

Table 4 Injection test: estimated parameters c and k and distances Estimation method

Parameter

Estimated values

R1 (m)

R6 (m)

PM/P

c k c k c k c k

1:77 107 m2 =s 9:25 1015 m2 =kPa s 1:66 107 m2 =s 8:20 1015 m2 =kPa s 4:92 108 m2 =s 6:34 1015 m2 =kPa s 4:63 108 m2 =s 6:05 1015 m2 =kPa s

0.62 0.64 0.64 0.56 0.97 0.52 0.95 0.53 0.58

1.09 1.06 1.05 1.20 0.70 1.31 0.71 1.29 1.17

LS/P PM/L LS/L Tentative distances

addition, the duration of a thermal load td employed in the estimation is 65820 s, 102120 s and 4311720 s, respectively, for the first, second and third heater tests. Tables 5–7 present the estimated values of X1 and X2 from the first, second and third heater test, respectively. The values of X1 and X2 for each thermistor in Tables 5 and 6 from the first two tests indicate that they are in a good agreement when compared within the same source model. When the results from the two models are compared, they are still within one order of magnitude. Note that the line source model always yields smaller values of X1 but higher values of X2 when compared to those from the point source model. More consistent results from different models are found in Table 7 for the third heater test which has the longest duration of a thermal load. It can be seen that the estimated values of X1 and X2 presented in Tables 5–7 show an expected trend. The largest values are found from PZ1 which is the closest piezometer to the heater whereas the values from PZ8, the farthest piezometer, are the smallest ones. The values of the thermal diffusivity cn and the lumped parameter L0 =kn can be obtained by performing a linear least square fit of the log/log plot between X1 and X2 ; respectively, and the tentative distances Ri ði ¼ 1; 4; 6; 7; 8; 9Þ: Once the values of cn and L0 =kn are determined, the values of the distance R can be computed from (13).

The estimated values of cn and L0 =kn as well as the computed value of the distances between the center of the test zone and the six thermistors (R1 ; R4 ; R6 ; R7 ; R8 and R9 ) together with the tentative distances estimated from the URL Borehole Collar Survey Information Report are presented in Tables 8–10 for the first, second and third heater tests, respectively. C.2.2. Pore pressure response The data from four working piezometers (PZ1, PZ3, PZ5 and PZ6) for all three heater tests are used in modeling the pore pressure field to determine the lumped parameters G1 and G2 given by (18). Two methods of parameter estimation are once again employed in the back-analysis for two types of models, point and spherical source models. For the nonlinear least square method, the number of observations used in the parameter estimation from each piezometers are 600, 900 and 2800, respectively, for the first, second and third heater test. In the estimation of G1 and G2 ; the lumped parameters U1 and X1 as well as L0 =kn are taken from the interpretation of the injection and heater tests (temperature field) in the case of a point source model. Note that for PZ1 and PZ6, the modification has been made on the distances for U1 and X1 ; due to the difference in expected distances for each test. For PZ3 and PZ5, the values of

ARTICLE IN PRESS 1408

E. Detournay et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1395–1411

Table 5 First heater test, temperature field: estimated lumped parameters X1 and X2 Estimation method

Lumped parameters

PM/P

X1 X2 X1 X2 X1 X2 X1 X2

LS/P PM/L LS/L

Unit

PZ1-T

PZ4-T

PZ6-T

PZ7-T

PZ8-T

PZ9-T

pffiffiffiffiffiffiffi

103 1=s 1C pffiffiffiffiffiffiffi

103 1=s 1C pffiffiffiffiffiffiffi

103 1=s 1C pffiffiffiffiffiffiffi

103 1=s 1C

2.76 138.51 2.58 157.89 1.75 183.93 1.83 171.84

2.02 117.41 1.97 130.70 1.45 152.49 1.50 139.95

1.25 86.63 1.22 97.41 1.14 93.41 1.10 101.24

1.37 89.07 1.34 100.33 1.19 101.33 1.15 105.30

1.02 72.81 0.99 81.84 0.98 75.73 0.95 83.32

2.36 127.11 2.20 148.11 1.77 156.28 1.76 157.64

Table 6 Second heater test, temperature field: estimated lumped parameters X1 and X2 Estimation method

Lumped parameters

PM/P

X1 X2 X1 X2 X1 X2 X1 X2

LS/P PM/L LS/L

Unit

PZ1-T

PZ4-T

PZ6-T

PZ7-T

PZ8-T

PZ9-T

pffiffiffiffiffiffiffi

103 1=s 1C pffiffiffiffiffiffiffi

103 1=s 1C pffiffiffiffiffiffiffi

103 1=s 1C pffiffiffiffiffiffiffi

103 1=s 1C

2.49 125.35 2.31 140.58 1.44 170.25 1.52 156.87

1.83 105.99 1.74 121.25 1.22 142.10 1.30 126.89

1.16 81.08 1.14 89.46 1.04 88.78 1.01 92.64

1.25 91.23 1.23 101.82 1.03 107.23 1.02 106.99

0.88 87.54 0.94 85.89 0.83 91.58 0.87 89.19

2.25 112.12 1.96 137.39 1.57 139.70 1.53 144.66

Table 7 Third heater test, temperature field: estimated lumped parameters X1 and X2 Estimation method

Lumped parameters

LS/P

X1 X2 X1 X2

LS/L

Unit

PZ1-T

PZ4-T

PZ6-T

PZ7-T

PZ8-T

PZ9-T

pffiffiffiffiffiffiffi

103 1=s 1C pffiffiffiffiffiffiffi

103 1=s 1C

1.92 171.11 1.54 152.06

1.54 138.20 1.30 126.13

1.05 96.05 0.99 93.38

1.12 104.72 1.02 100.17

0.87 82.34 0.85 81.40

1.72 157.0 1.49 145.20

Table 8 First heater test, temperature field: estimated parameters cn and L0 =kn and distances Estimation method

Parameters

Estimated values

R1 (m)

R4 (m)

R6 (m)

R7 (m)

R8 (m)

R9 (m)

PM/P

cn L0 =kn cn L0 =kn cn L0 =kn cn L0 =kn

1:92 106 m2 =s 84.16 1C m 1:77 106 m2 =s 95.26 1C m 1:22 106 m2 =s 99.33 1C m 1:21 106 m2 =s 100.48 1C m

0.50 0.61 0.52 0.60 0.63 0.54 0.60 0.58 0.52

0.69 0.72 0.68 0.73 0.76 0.65 0.73 0.72 0.67

1.11 0.97 1.09 0.98 0.97 1.06 1.00 0.99 1.12

1.01 0.94 0.99 0.95 0.93 0.98 0.96 0.95 0.86

1.36 1.16 1.34 1.16 1.13 1.31 1.16 1.21 1.34

0.59 0.66 0.60 0.64 0.62 0.64 0.62 0.64 0.68

LS/P PM/L LS/L Tentative distances

ARTICLE IN PRESS E. Detournay et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1395–1411

1409

Table 9 Second heater test, temperature field: estimated parameters cn and L0 =kn and distances Estimation method PM/P LS/P PM/L LS/L

cn L0 =kn cn L0 =kn cn L0 =kn cn L0 =kn

Estimated values

R1 (m)

R4 (m)

R6 (m)

R7 (m)

R8 (m)

R9 (m)

1:80 106 m2 =s 86.82 1C m 1:66 106 m2 =s 96.60 1C m 1:03 106 m2 =s 104.57 1C m 1:06 106 m2 =s 102.07 1C m

0.54 0.69 0.56 0.69 0.70 0.61 0.68 0.65 0.57

0.73 0.82 0.74 0.80 0.83 0.74 0.79 0.80 0.72

1.16 1.07 1.13 1.08 0.98 1.18 1.02 1.10 1.18

1.07 0.95 1.05 0.95 0.99 0.98 1.01 0.95 0.91

1.52 0.99 1.37 1.12 1.22 1.14 1.18 1.14 1.38

0.60 0.77 0.66 0.70 0.65 0.75 0.67 0.71 0.73

Tentative distances

Table 10 Third heater test, temperature field: estimated parameters cn and L0 =kn and distances Estimation method

Parameters

Estimated values

R1 (m)

R4 (m)

R6 (m)

R7 (m)

R8 (m)

R9 (m)

LS/P

cn L0 =kn cn L0 =kn

1:32 106 m2 =s 105.32 1C m 1:04 106 m2 =s 98.92 1C m

0.60 0.62 0.66 0.65 0.57

0.75 0.76 0.78 0.78 0.72

1.10 1.10 1.03 1.06 1.18

1.03 1.01 1.00 0.99 0.91

1.32 1.28 1.21 1.22 1.38

0.67 0.67 0.68 0.68 0.73

LS/L Tentative distances

Table 11 First heater test, pore pressure field: estimated lumped parameters G1 and G2 Estimation method

Lumped parameters

Unit

PZ1

PZ3

PZ5

PZ6

PM/P

G1 G2 G1 G2 G1 G2 G1 G2

kPa/1C m — kPa/1C m — kPa/1C m — kPa/1C m —

913.84 0.40 797.23 0.06 847.90 0.15 792.48 0.14

— — 211.52 0.54 — — 187.86 0.19

— — 271.65 0.03 — — 272.78 0.11

504.07 0.48 522.26 0.09 547.92 0.12 507.69 0.17

LS/P PM/SP LS/SP

Table 12 Second heater test, pore pressure field: estimated lumped parameters G1 and G2 Estimation method

Lumped parameters

Unit

PZ1

PZ3

PZ5

PZ6

PM/P

G1 G2 G1 G2 G1 G2 G1 G2

kPa/1C m — kPa/1C m — kPa/1C m — kPa/1C m —

833.19 0.47 695.68 0.23 716.88 0.18 681.83 0.18

— — 248.20 0.36 — — 237.59 0.14

— — 236.13 0.06 — — 240.58 0.18

514.47 0.48 506.95 0.37 492.96 0.14 447.39 0.21

LS/P PM/SP LS/SP

ARTICLE IN PRESS 1410

E. Detournay et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1395–1411

Table 13 Third heater test, pore pressure field: estimated lumped parameters G1 and G2 Estimation method

Lumped parameters

Unit

PZ1

PZ3

PZ5

PZ6

PM/P

G1 G2 G1 G2 G1 G2 G1 G2

kPa/1C m — kPa/1C m — kPa/1C m — kPa/1C m —

758.33 0.25 712.14 0.15 743.42 0.23 776.83 0.25

252.55 0.39 230.17 0.32 219.36 0.24 217.13 0.25

251.42 0.12 246.03 0.12 250.39 0.12 247.03 0.13

427.59 0.14 417.46 0.09 425.44 0.13 426.26 0.14

LS/P PM/SP LS/SP

Table 14 First heater test, pore pressure field: estimated parameters g and a and distances Estimation method

Parameters

Estimated values

R1 (m)

R3 (m)

R5 (m)

R6 (m)

PM/P

g a g a g a g a

475.44 kPa/1C 0.29 m 452.17 kPa/1C 0.11 m 451.44 kPa/1C 0.09 m 435.61 kPa/1C 0.17 m

0.52 0.73 0.57 1.93 0.53 0.59 0.55 1.23 0.43

— — 2.14 0.20 — — 2.32 0.89 2.25

— — 1.66 3.33 — — 1.60 1.49 1.77

0.84 0.61 0.87 1.20 0.82 0.74 0.86 1.02 1.02

LS/P PM/SP LS/SP Tentative distances

Table 15 Second heater test, pore pressure field: estimated parameters g and a and distances Estimation method

Parameters

Estimated values

R1 (m)

R3 (m)

R5 (m)

R6 (m)

PM/P

g a g a g a g a

469.25 kPa/1C 0.34 m 471.19 kPa/1C 0.30 m 425.99 kPa/1C 0.11 m 439.54 kPa/1C 0.21 m

0.56 0.72 0.68 1.33 0.59 0.63 0.64 1.16 0.48

— — 1.90 0.84 — — 1.85 1.53 2.29

— — 1.99 — — — 1.83 1.18 1.82

0.91 0.71 0.93 0.83 0.86 0.81 0.98 1.02 1.07

LS/P PM/SP LS/SP Tentative distances

U1 and X1 ; which are given by (7) and (13), are obtained by using the estimated values of c (from Table 4) and cn (from Tables 8–10), respectively, and the expected distances from the sensors and the test zone. Tables 11–13 present the estimated values of G1 and G2 from pore pressure data of the first, second and third heater test, respectively. Note that the application of peak matching technique at PZ3 and PZ5 in the first and second tests is very difficult due to the uncertainty of selecting a single peak value during the extended period of time over which the peak occurred (see Figs. 10 and 12 in Ref. [1]). Remember that PZ3 and PZ5 are the two furthest piezometers in this case,

and as a result, reasonable results were not obtained from peak matching (e.g., in some cases G2 41; which is not possible). Therefore, the estimated values from PZ3 and PZ5 from the peak matching method in the first two tests were not used in the estimation of g and a. The estimated values of G1 (Tables 11–13) show a more reasonable trend with respect to distance than those of G2 : In most cases, the largest values of G1 are found from PZ1 followed by PZ6, PZ5 and PZ3; that is, from the nearest to furthest piezometer. Once again, the most consistent results are obtained from the longest test, which is the third heater test, as shown in Table 13. In addition, the spherical source model

ARTICLE IN PRESS E. Detournay et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 1395–1411

1411

Table 16 Third heater test, pore pressure field: estimated parameters g and a and distances Estimation method

Parameters

Estimated values

R1 (m)

R3 (m)

R5 (m)

R6 (m)

PM/P

g a g a g a g a

458.14 kPa/1C 0.24 m 435.71 kPa/1C 0.18 m 439.14 kPa/1C 0.20 m 441.57 kPa/1C 0.22 m

0.60 0.98 0.61 1.19 0.59 0.90 0.57 0.89 0.48

1.81 0.62 1.89 0.56 2.00 0.87 2.03 0.88 2.29

1.82 1.97 1.77 1.51 1.75 1.73 1.79 1.71 1.82

1.07 1.78 1.04 2.13 1.03 1.57 1.04 1.60 1.07

LS/P PM/SP LS/SP Tentative distances

yields more consistent results than the point source model. The values of the hydro–thermal coupling coefficient g and the effective drainage radius a can be obtained by performing a least square fit between G1 and G2 ; respectively, and the tentative distances Ri ði ¼ 1; 3; 5; 6Þ: Once the values of g and a are determined, the distance R can be computed from (18). Tables 14–16, list the estimated values of the parameters g and a from the first, second and third heater tests for point and spherical source models. References [1] Berchenko I, Detournay E, Chandler N, Martino J. An in-situ thermo–hydraulic experiment in a saturated granite I: design and results. Int J Rock Mech Min Sci 2004, 41(8).

[2] Beck J, Arnold K. Parameter estimation in engineering and science. New York: Wiley; 1977. [3] Bard Y. Nonlinear parameter estimation. New York: Academic Press; 1974. [4] McTigue DF. Thermoelastic response of fluid-saturated porous rock. J Geophys Res 1986;91(B9):9533–42. [5] Coussy O. Poromechanics. New York: Wiley; 2004. [6] Biot MA. General theory of three-dimensional consolidation. J Appl Phys 1941;12:155–64. [7] Detournay E, Cheng AH-D. Comprehensive rock engineering, vol. 2, Fundamentals of poroelasticity, New York NY: Pergamon, 1993, p. 113–71, [chapter 5]. [8] Chan T, Jeffrey JA. Scale and water-saturation effects for thermal properties of low-porosity rock. In: Association of Engineering Geologists Proceedings of 24th US Symposium on Rock Mechanics; 1983. p. 287–301.