A mathematical model and numerical investigation for determining the hydraulic conductivity of rocks

Int. J. Rock hfech. Min. Sci. Vol. 34, No. 5, pp. 741-159,1997 0 1997 Elscvier ScienceLtd.All righta resewed Printed in GreatBritain 0148~9062197 $17.00+0.00

PII: 50148-9062(~9

A Mathematical Model and Numerical Investigation for Determining the Hydraulic Conductivity of Rocks D. L. D. B. R.

LESNICt ELLIOTT? B. INGHAMtS CLENNELLg J. KNIPEg Transient flow analysis has numerous geophysical and geotechnical applications, and is necessary for the understanding and development of laboratory measurement techniques in rock and soil mechanics, particularly permeability testing. In this paper, a hydraulic pump flow permeability test is mathematically modelled reducing the problem to the time-deparordentdi#usion equation with unknown spacewise dependent material property coe@cients. This equation is solved subject to certain boundary conditions and, with limited additional available boundary measurements in time, of the pressure and/or hydraulic flux. Both experimental and numerically simulated noisy data are inputted into the formulation of the inverse problem. The numerical solution of the inverse problem obtained from a finite-dlrerence method (FDM), combined with an ordinary least-squares approach, is compared with the exact solution for homogeneous and heterogeneous materials with constant, linear and quadratic spatially dependent hydraulic conductivity, subjected to one-dimensional, pump-flow hydraulic boundary conditions. The numerically obtained results, based on an optimal selection of the limited number of time measurements, are consistent with the sensitivity analysis performed prior to inversion. The ident@ability and stability of the inverse problem are also numerically investigated. 0 1997 Elsevier Science Ltd

S” SS X

NOMENCLATURE A

B, C I K

K* Ko, L, L LS M, iv

2” Ri Sens(+)

K2

cross-sectional area of the sample non-dimensional coefficients total number of time measurements hydraulic conductivity of the sample mean average of the hydraulic conductivity hydraulic conductivity coefficients length of the sample least-squares functional number of finite-difference space and time elements space of quadratic polynomials inflow rate ratios of the normal&d sensitivity coefficients normal&d sensitivity coefficients

V” ai, bi, ci,

rii

c

i) m P Pr Pm 40, 41 t* tjt

Xi

6,o V

TDepartment of Applied Mathematical Studies, University of Leeds, Leeds LS2 9JT, U.K. $Author to whom correspondence should be addressed. $Department of Earth Sciences, University of Leeds, Leeds LS2 9JT, U.K.

ox

,

L\t

Greek

symbols

a, P”, Y. 741

compressive storage of the upstream reservoir specific storage of the sample mean average of the specific storage volume of the upstream reservoir finite-difference method coefficients fluid compressibility optimal time measurements fluid mass fluid pressure reference pressure values of the pressure at (0, fm) hydraulic fluxes time at which the pump is switched off nodes of the finite-difference grid time measurements space and time coordinates volumetric rate of flow per unit cross-sectional area space and time finite-difference meshes.

dimensionless parameters

LESNIC et al.:

HYDRAULIC

specific weight of the fluid constants porosity of the sample permeability of the sample zero or unity coefficients dynamic viscosity of the fluid fluid density reference density weight factor in the finite-difference method.

direct input pump pump

numerically obtained data measured data is on is off

1. INTRODUCTION

Hydraulic conductivity is associated with petroleum, civil and mining engineering. However, the petroleum and other areas are not mutually exclusive, and civil and mining engineering also require knowledge of the hydraulic conductivity in order to understand a variety of problems, e.g. possible locations of underground storage sites for contaminated fluids. This paper is written in the context of petroleum engineering, but the results also apply to modelling in civil and mining engineering. The flow of fluid within petroleum reservoirs is strongly dependent upon both the heterogeneity and/or the anisotropy of the rock structure. However, any attempt to accommodate both in an approach which is trying to retrieve internal information from boundary measurements alone would be highly complex and, at the moment, virtually impossible. In addition to the above properties, one would like to be able to handle fractured rocks, as well as intact rocks, and it is believed that an extension of the present work will be capable of this. However, any such inclusion will still be maintained within the mathematical continuum model presented in this paper. Thus, the thrust of this paper is modelling rock as a continuum. From the fluid-flow view, anisotropy and heterogeneity are closely related properties, and can be presented by a single approach if one views the hydraulic conductivity to be a tensor, whose principal directions, magnitudes and variation are functions of the space coordinates. Clearly, what is required in the field is a detailed knowledge of this tensor and the question is whether, and how well, information from laboratory small scale sample experiments and numerical calculations from mathematical models, can help to achieve this. One of the main problems is that experiments, simply by their laboratory modified constraints, e.g. the sides of the sample being confined within an impermeable material, make the scaling up of any results difficult. With any mathematical model needing to be initially validated by experiments, it clearly must be subject to the same restrictions as the experiments. It is well-known that the flow pattern in a rectangular sample is significantly changed by the no-flow condition

CONDUCTIVITY

OF ROCKS

imposed on the side walls, according to both the orientation of the sample axis with respect to the principal directions of the hydraulic conductivity tensor and the ratio of the principal values of the hydraulic conductivity. Hence, it is essential that in discussing any one-dimensional model its limitations in any scale up prediction be understood. In general, laboratory permeability tests are performed based on constant-head [l], transient (pulsedecay) [2] and constant flow-rate (pump flow) [3] methods. In this study, since low permeability (tight) rocks are tested, the pump-flow method is employed. Moreover, in real time the pressure is easy to measure very accurately using electron transducers. Also, a very low flow rate is much easier to enforce using a pump than to measure using some outflow volume change measuring device that must take some average, rather than instantaneous value. The plan of the paper is as follows. The hydraulic equipment, experimental method, and pump-flow system performances and limitations are described in Section 2. In Section 3 a mathematical model which describes the hydraulic pump flow experiment is developed. It is shown that for cases where the fluid flow is creep, the propagation of changes in pressure through a porous rock can be represented by the heat (diffusion) equation recast in terms of the pressure. Although for creep flow it might be possible that a rapid acceleration caused by a step load can yield a linear inertial contribution [4], in most practical cases this linear inertia is negligible. The heat flow paradigm is to be expected and other studies in the heat conduction context by Garnier et al. [5] and Dowding el al. [6] considered similar situations in which the inflow rate of mass (fluid or heat) could be supplied, or not, by a hydraulic pump or heater (on and of?+). However, in the hydraulic experiment tested in this study the boundary conditions are different and also the properties of the material possess a spatial dependence rather than a functional solution (pressure or temperature) dependence, as assumed in the previous studies. Further, the boundary conditions derived in Section 3 are different from those encountered in the pump-flow tests of Morin and Olsen [3], and Olsen et al. [7], where a constant hydraulic flux at the upstream face of the sample (rock) is prescribed. However, the constant hydraulic flux boundary condition is not physically achieveable in a real experiment as the upstream storage of an upstream reservoir which is maintained in hydraulic contact with the sample is never zero. The treatment of Morin and Olsen [3] is strictly valid only for the case where the specific storage of the sample dominates the compressibility of the system, e.g. for large soil samples, and not for the small samples or rigid rock specimens [8] investigated in this study. If the hydraulic properties of the sample are known then the mathematical model described in Section 3 formulates a direct well-posed problem which can be solved numerically using an FDM (Section 4). In Section 5, the sensitivity coefficients with respect to the hydraulic

LESNIC et al.:

HYDRAULIC

properties of the rock are briefly introduced as their ratios may characterise the degree of correlation, i.e. of identifiability, of the unknowns in an inverse identification problem. In Section 6, the inverse formulated problem requires the recovery of heterogeneous properties of materials, e.g. the spatially dependent hydraulic conductivity of the rock, only from boundary measurements of pressure and/or hydraulic flux. Numerous engineering and mathematical researchers have considered problems equivalent to estimating spatially dependent hydraulic conductivity, e.g. Chen and Seinfeld [9], and Kitamura and Nakagiri [lo]. More recently, Huang and Ozisik [l l] have been concerned with the simultaneous estimation of the linearly varying thermal conductivity and heat capacity which are the analogous heat paradigm for the hydraulic conductivity and specific storage. However, these studies used multiple internal spatial and temporal measurements made during the transient diffusion, whilst in this study only boundary measurements are allowed. This restriction on the inverse formulation was previously taken into account by Ciampi et aE. [12] who proposed an approximate direct solution applicable to slightly heterogeneous samples of linear and exponential spatially varying type for the use with flash-type experiments, utilising temporal measurements made only at an insulated back surface of the slab sample. However, our hydraulic pump-flow experiment has different boundary conditions on the sample, i.e. it is not of the flash type and, therefore, this method cannot be utilised. In this study, the limitation of the boundary measurements only is imposed since internal measurements in practice require a special core holder with pressure ports leading to certain points on the outside of the sample and, in fact, internal measurements in real petroleum reservoirs require new wells to be drilled, each costing several million pounds. In addition, measurements in the well are likely to be significantly noisy, being affected by sources of errors such as: (i) uncertainty that the reservoir interval under test has been sealed off perfectly; (ii) unknown wellbore storage effects made worst by gas exsolving in the wellbore; (iii) skin effects: anomalously low permeability in the zone surrounding the well due to the damage and mud infiltration during the drilling process. The inverse problem described in Section 6 is further parameterised by assuming that the hydraulic conductivity is a polynomial function in the space coordinate with unknown coefficients to be determined and a steadystate mathematical analysis is performed. These coefficients are then calculated numerically by employing an ordinary nonlinear, least-squares method which naturally minimises the gap between the measured and iteratively calculated values of the boundary pressure and hydraulic flux at sampling times which are prescribed based on an optimal criterion. Finally, in Section 7, for a wide variety of rocks and technical operating conditions, the numerical results are initially compared with the experimental values for the pressure increase and decay curves at the upstream face

CONDUCTIVITY

OF ROCKS

143

of the sample, showing a reasonable validation of the forward mathematical model introduced in Section 3. Next, the inverse analysis as described in Section 6 is performed for the retrieval of constant, linear and quadratic spacewise hydraulic conductivities. In the numerical investigation the number of time measurements is imposed to be equal to the number of unknowns, in order to analyse the most general case of parameter identifiability and operational modes of the pump. Exact, experimental and noisy numerically simulated measurements are inputted into the inverse formulated problem. The interpretation of the numerical results in comparison with their exact values is discussed and argued on the basis of the sensitivity analysis performed prior to inversion. 2. EQUIPMENT

AND EXPERIMENTAL

METHOD

In the hydraulic experiments performed in this study, the mean average permeability is measured using the constant-rate-of-flow (pump-flow) technique [ 131, which has the advantage of an easy set-up and computer control, a high level of accuracy at low permeabilities and, in general, short testing times, i.e. a short transient period before the steady-state is reached. The syringe infusion pump-flow method also enables the permeant to be easily changed, with only a small volume of per-meant required. The same simple apparatus shown in Fig. 1 can also be used for pulse-decay measurements [2] and other types of transient analysis tests, e.g. oscillatory tests [14]. Samples to be tested are first saturated fully in the chosen permeant solution (pure water or formation brine) by immersion in de-aired water under vacuum for several hours. Both standard core plugs and more irregularly shaped samples can be tested for permeability. Cylindrical samples (2.5-5.5 cm cores, 210 cm long) are mounted in a close fitting polypropylene sleeve and placed in a stainless steel sample holder. Prismatic samples sawn from irregular or fragile cores can also be experimented on by first mounting them in an elastic mouldable polymer which is formed into a cylinder before insertion into the sample holder. Dimensions of the samples are measured using a micrometer to a precision of 0.01 cm and a probable accuracy better than 0.025 cm. The confining pressure is set at a constant value of 100 kPa, and is maintained constant at this value by a servo pump. A constant back-pressure of at least 300 kPa, which is maintained in order to ensure that the sample is fully saturated, is provided via a second pressure generator or doubly regulated air line that feeds an air-water exchange burette at the outflow end of the system. Flow rates between 1.2 x lo-’ and 0.3 cm3/sec can be set via the computer control system, and changed at any time during the test. Flow rates have been calibrated by two independent methods (flow into a graticule and flow into a back-pressured volume change device). The pressure drop induced across the sample by the flow of permeant is measured with an accurate

LESNIC ef al.:

144

HYDRAULIC

differential pressure transducer that has a range of 340 kPa. The flow rate of the pump is generally chosen so that the equilibrium pressure is between 10 and 120 kPa. The back-pressure and differential pressure are monitored continuously during the test and logged to a computer via an A/D converter system with 12-bit precision. An independently calibrated A/D converter gives an additional readout that is used to verify the accuracy of the logged data. The differential pressure precision is 0.1 kPa with an accuracy better than 0.5%, except at pressures lower than 3 kPa, which are not used in the analysis. The back pressure remains constant to within the measurement precision throughout the tests. Temperature is found to affect the pressure readings as water expands and contracts in the pipes, so the laboratory is maintained close to 20°C in winter and 25°C in summer. All the experimental results reported in this paper were acquired for samples with a permeability greater than lo-l3 cm’, which are essentially unaffected by temperature changes. 2.1. System performances and limitations Using the equipment and procedures previously described, it is possible to measure the permeability of rocks over a wide range between lo-l6 and lo-* cm’. The upper limit of this range is set by the syringe volume and the tendency of the flow to be non-darcyan, i.e. at high rates the flow enters the Forscheimer regime where inertial effects become important, whilst the lower limit is imposed by the temperature variation effects in the laboratory, which however could be improved in a sealed, air-conditioned environment. The time taken to reach equilibrium in the constant-rate pump-flow test and the time scales of all the transient phenomena that we describe, are largely dictated by the elastic compliance of the system to fluid

B

Permeant reservoir

CONDUCTIVITY

OF ROCKS

pressure changes. Short testing times are obtained by minimising this compliance and this can be ensured by the design having short, rigid lengths of thick-walled stainless steel tubing and a small upstream dead volume in the sample holder. The upstream compliance is a system constant and can be measured by replacing the sample with a steel blank and pumping permeant into the upstream tubing and dead volume at a constant rate, whilst the upstream pressure is monitored. The downstream compliance of the system described here is essentially infinite as the back-pressure is buffered at a constant value with an air line. Other configurations are possible with the same equipment, but they are not described here. While the compliance of the existing system is quite small for its type, the samples are generally of small size and have low porosity, so that the measured transient response and the overall elastic storage of the system are dominated by the upstream compliance of the apparatus [15,16]. This means that whilst the permeability can be measured with great accuracy, the exact value of the sample elastic storage coefficient and so the specific storage of the rock material, cannot be measured accurately with the existing apparatus. Indeed, for purposes of analysis the sample storage S, will be assumed small and known, and all the meaningful remaining storage will be ascribed to a lumped upstream storage coefficient S.. This is somewhat in contrast with the results of Neuzil et al. [15] for sudden pressure increase experiments at the upstream of the sample, who showed that the measured pressure response is sensitive to both or none of these two storages. However, as will be more clearly numerically shown in Sections 4 and 7 for the pump-flow mathematical model which is developed in the next section, the pressure measurements at the upstream face of the sample are found to be

Confining pressure generator Gauge pressure transducer Sample holder

Back pressure

Syringe infusion pump

’Rock core sample Differential pressure

Fig. 1. Diagram showing the hydraulic experiment and equipment.

LESNIC et al.:

HYDRAULIC

insensitive to the values of S,, whilst being sensitive to the value of S, because of the rate of flow which is supplied at the entrance of the upstream reservoir. This conclusion on the storages will be justified numerically in Section 7, but this apparent limitation of the experiment will be tested experimentally in the future and possibly overcome by designing a new system in which the samples are much larger, and the dead and tube volumes are very much smaller. 3. MATHEMATICAL

CONDUCTIVITY

OF ROCKS

small constant and we 4 of the sample is time equation (4) becomes which is often used in

145

further assume that the porosity independent, then the governing the diffusion (heat) equation, petroleum engineering, i.e.

s, t$=V*(Kvp) where S, is the specific storage of the medium and is defined by s, = yw&.

MODEL

(7)

Equation (6) represents a slightly restricted class of model, being formed from the more general situation, where 4 = 4(x, y, z)exp(X@ - p,)) and K = K(x, y, z)exp(q(p - p,)), by assuming both parameters x 3.1. Governing equation and r~ to be zero. The governing (heat, diffusion) equation which is From the experimental viewpoint, the greater the ratio often used in petroleum engineering, for the laminar, of the sample length to sample diameter the less is the isothermal, and single-phase flow of most fluids of influence of any cross-sectional heterogeneity, but since constant and small compressibility, in an isotropic and it is only possible to determine an overall constant value saturated, heterogeneous, porous medium with the for the hydraulic conductivity the fact that some of the influence of gravity neglected can be obtained by experimental sample diameters are approaching the combining the mass conservation equation, Darcy’s law sample lengths is not relevant. In addition, experimental and an equation of state for the fluid [17], as follows. details are utilised only in order to validate the The conservation of mass equation for a fluid in a mathematical model for the case of constant hydraulic porous medium, i.e. sample (rock), can be expressed as conductivity and to investigate any effect of the specific storage of the sample on that situation. As no a vov)+jj$#v)=o (1) experimental details are available for situations involving axial variation of the hydraulic conductivity, for where p is the fluid density, 4 is the porosity of the which it would be essential to ensure that the length to medium and g is the volumetric rate of flow per unit diameter ratio of the sample was large, a one-dimencross-sectional area. sional mathematical model was assumed, independent of For laminar flow at low Reynolds numbers, Darcy’s the sample diameter, with conditions uniform across the law can be expressed as cross section of the sample. Then, for the experimental set-up described in Section 2 we can consider the flow to g= -Kvp (2) be one-dimensional along the length L of the sample and where p is the fluid pressure, K is the hydraulic equation (6) recasts as the one-dimensional heat conductivity of the medium which is assumed to be equation with spacewise dependent coefficients, i.e. isotropic and, for simplicity, the influence of gravity has been neglected. In equation (2) the hydraulic conduca , S(x) at (x, t) = ax tivity is given by In this section the pump-flow experiment described in the previous section is mathematically modelled.

ap

K = W/P

(3)

where ICis the permeability of the rock, p is the dynamic viscosity of the fluid and yWis the specific weight of the fluid. Combining equations (1) and (2) results in V.(PKvP) =z

In order equation required. expressed

a

(h).

(4)

to complete the description of the governing an equation of state for the pore fluid is For single-phase flow a suitable relation can be as [ 171, P = pr exp(c@ - PJ)

(5)

where p, is the value of the fluid density at some reference pressurep, and c is the fluid compressibility. Equation (5) is a good description for isothermal flow of most fluids of constant and small compressibility c. Then, if c is a

(x7t) E (0, L) x (0, a).

(8)

Equation (8) has to be solved subject to initial and boundary conditions which are prescribed in the next subsection. 3.2. Boundary conditions

The boundary conditions are derived for the one-dimensional transient pump-flow hydraulic test described in Section 2, which is drawn schematically in Fig. 2. In this pump-flow experiment an inflow of mass is applied using a pump at the entrance of an upstream reservoir which is in hydraulic contact with the rock sample, and the induced pressure is measured with a differential transducer which is connected between the upstream and downstream reservoirs. After the transients have past, or even before that, the pump may be switched off in order to record both pressure increase

746

LESNIC er al.:

HYDRAULIC

Differential Pressure Transducer

CONDUCTIVITY

OF ROCKS

where S, is the compressive storage of the upstream reservoir which is defined by S, = Y&V”.

Sample

Upstream storage, Su

At this stage, it is worth noting that the boundary condition (13) was previously used by Lin [ 181 for modelling the transient pulse test designed by Brace et al. [2], but its derivation was argued on the basis of a private communication. In addition, the actual experiment is more general than the pulse experiment modelled by Hsieh et al. [19] for which Qu = 0.

Downstream reservoir

k? Ss?

(14)

G[

3.3. Non-dimensional equations Time Fig. 2. Schematic diagram showing the hydraulic experiment arrangement of the pump-flow test described in Section 2 and illustrated in Fig. 1.

Before performing the numerical calculations, the more general mathematical model to that described in Sections 3.1 and 3.2 is non-dimensionalised, as follows:

t=t=i, and decay measurements which are necessary as additional information in the inverse analysis of retrieving the hydraulic properties of the rock which will be finally undertaken in this study. At the start of the experiment, the pressure distribution in the sample (rock) is constant and is defined to be zero. i.e. p(x, 0) = 0,

x

E

[O, L].

p(L, t) = 0

K = K*I?,

l E (0, Lc]

(10)

where t, represents the time after which the pump is switched off, which is usually taken at the steady-state equilibrium. At the upstream face of the sample, x = 0 the conservation of mass at the sample-reservoir interface is applied. The fluid mass in the upstream reservoir is given by m = Y&V”

(11)

where Vu is the volume of the upstream reservoir. The conservation of mass at the sample-reservoir interface, i.e. that the rate of decrease of fluid mass in the reservoir is equal to the rate of mass flow into the sample, is given

the dimensionless

at x = 0,

t E (0, tm]

S: = ;

t E (0,

S(x) dx,

sII

s0

L4(x)

gb* = ;

dx.

(16)

scl

The governing non-dimensional parameters governing partial differential equation are B

=

in the

K*t=

(yw(x + c)L*($*)’ c = (c + vr)(P=-Pr)

(17)

where l/B represents the coefficient of the partial derivative term with respect to time and C the coefficient of the square of pressure gradient term. Justification for neglecting the latter in order to formulate equation (8) was based upon Cc< 1, which appears acceptable. For 1 = 0 the parameter B can be expressed as B = K*t,l (S? L2) and is the only non-dimensional quantity arising within the governing partial differential equation, the other non-dimensional parameters are being confined to the boundary condition at x = 0. Then, in non-dimensional form equations (8)--(10) and (13) become (x, t) E (0, 1) x (0, l]

p(x,O)=O,

XEP, 11

P(l, t) = 0,

t E (0,

(18)

tzl (13)

11

(19)

(20)

z (0, t) = /&K(O) g (0, t) + yu, t E (0, 11 (21) where, for simplicity, the bars have been dropped and the dimensionless coefficients are given by a = B + K*t,IL’SS*,

ap(0, t) = K(O)A zaP (0, t) + Qu,

quantities,

(12)

where A is the cross-sectional area of the sample, v is the one-dimensional component of the volumetric rate of flow per unit cross-sectional area and QU is the inflow rate at the entrance of the upstream reservoir. Introducing equations (2) and (5) into the conservation of mass given by equation (12), and using equation (1 l), results in Sat

(15)

L

K(x) dx,

,

= pAv - pQu

4 = qY’$

I. K* +

by -$

p=p%p,

S, = SF&

where the bars denote pm = ~(0, tx) and

(9)

The downstream face of’the sample, x = L is maintained at its initial pressure value, i.e. zero, thus

x=LZ,

p,, = K*At,ILS,,

yu = tmQu/Sup,.

(22)

LESNIC et al.:

HYDRAULIC

The nature of the solutions is governed by the magnitude of the non-dimensional parameters tl, flu and yU. As the parameter a is of 0(104) in all the present work, which corresponds to either a relatively large value of K* (Figs 3 and 4) or a relatively long dimensional pressure rise time (Fig. 5), then the transient diffusion

CONDUCTIVITY

747

OF ROCKS

(a) 18 16

(a) 18 16

0

10

20

30

40

50

60

70

80

90

(b)

0

10

20

30

40

50

60

(b)

0

40

20

60

80

100

120

140

time (seconds) Fig. 4. Numerical (-) and experimental (xxx) values for the pressure L = 5.6 cm, decay curves, increase and Pa 113 when and (a) Q. = 1.6 x 10-j cm)/sec; (b) K* = 2.1 x lO-6 cm/set Q. = 8.0 x IO-’ cm’/sec.

0

10

30

20

40

50

6b

(cl 50

process is quasi-static and the pressure at the steady-state is expressible in a closed form for any given K(x) from the steady-state form of equations (18)-(21), i.e.

=

45

0,

40 35 G? 2 30

E

(0, 1)

0, /!?.K(O)$O) + yu= 0.

p(l)=

e 25 P &8 20

x

(23)

(24)

Such steady-state solutions are presented in Section 6.1 for various hydraulic conductivities K(x) and they approach (Figs 3-6) the quasi-static pressure values p(x, t) as f + co. 0

10

20

30

40

50

60

70

time (seconds) Fig. 3. Numerical (-) and experimental (xxx) values for the pressure curves, when L=3cm, increase and decay Pa 0, K* = 1.35 x 10-6cm/sec and (a) Q. = 1.6 x 10-3cm)/sec; (b) Q. = 3.2 x lo-‘cm3/sec; (c) Q. = 4.8 x IO-‘cm’/sec.

4.

DIRECT SOLUTION

The first step in the inverse analysis of identifying the hydraulic properties of the rock is the development of the corresponding direct solution for the problem given by equations (18)-(2 1).

748

LESNIC et al.:

HYDRAULIC CONDUCTIVITY OF ROCKS

(a)

In this study, the numerical method adopted for solving the initial boundary value, direct problem given by equations (18)-(21) is a weighted average finite-difference method (FDM), using a weight factor t? E [0, 11, and uniform space and time grids, similar to that used by Huang and Ozisik [l l] for solving a Neumann inverse problem in heat conduction. If Xi = i Ax for i = 0, M, and tj = j At for j = 0, N, are uniform, linear space and time grids of the interval [0, 11, respectively, then the resulting finite-difference equations (18)-(21), in matrix form, become bo al 0

co b,

0 cl 0

[: 0

.

0

‘0 0 CM-Z

aM-2 h-2 0

...

aM-1

(b)

*IO4

0

b-1 :I -I

3

a 23

-2

3

x

~~~~~~~

=

[jj

t2%

with the initial and boundary conditions given by to) = 0, i= 0,A4

p(Xi,

(26)

/I(x~, tj) = 0, j = 0,N

(27)

(c)

*lo4

where a ! =ae[K(xi)+K(xi-l)l 2Ax ’

i= 1 (M_ 1) ’

(28a)

c,I= aWWJ + K(xi+111 i = 0 (M _ 2) 2Ax ’ ’

(28b)

b. = _a+&Q)

+ K&o)] _ %Xo) Ax

2Ax

(28~)

2 At

I

I*104

0’100

time b 3

70

‘;;

60

t2

50

g

40

Lk-_

30 20

su=1.7x10-4

cmxcm

10

su=1.4x10-4

cmxcm

08100

l’iO4

2*;04

3*io4

Sw2.Oxlo-4

cmxcm

Su=2.3x10-4 Su=2.6xlO-4

cmxcm ~rnx~rn

4*io4

2.6} x IU-‘cm‘.

I

I

3*104

4.104

(seconds)

Fig. 6. The numerical values for the (a) pressure ~(0, I), (b) hydraulic flux go(r) and (c) hydraulic flux ql(r), for constant K(x) = L = 1.48 x lo-‘cm/set (---), linear K(x) = & + Kix = (1.11 + 0.21.x) x IO-scm/sec (-) and quadratic K(x) = & + Rx + KZX*= (0.21 + 0.37x + 0.16~‘) x lo-* cm/set (. . .), hydraulic conductivities.

bi

=

_ae[K(xi+i) + 2K(xI) +

K(xi-

I)]

2Ax

time (seconds) Fig. 5. Numerical (-) and experimental (000) values for the pressure increase and decay curves, I), when __. ~(0, _^1 . S. E { 1.4, 1.7, 2.0, 2.3,

I

21104

_

S(Xi)

AJC

At

’

i= l,(M-

1)

(284

LESNIC et cd.: HYDRAULIC CONDUCTIVITY OF ROCKS

-

e)[md + wxo)l 2Ax

6. INVERSE PROBLEM FORMULATION

P(XO, _ tj) BuAt

‘“Qp

41 -

ey; +w1,(,,,

tj)

(2*e)

”

d/ =

_‘(l - Q[K(Xi)+

K(Xi-

2Ax

-

l)jp(xi_, 9

a(1- e)[K(xi)+ WI+ 111

p(xi+

,,

tj)

tj)

2A.x +

a(1 - e)[K(xi+

1) + X(X~)

-l- K(Xi_ I)]

p(Xi, 2Ax tj),

- Ss(xi) Ax At

i = 1, (M - 1)

In the first instance, the inverse analysis would require the determination of both the spatially varying hydraulic conductivity and specific storage. However, since the main purpose is to determine the permeability of the sample, only the hydraulic conductivity needs to be identified. In addition, for the actual experimental design, the value of the specific storage St cannot be measured exactly, nor can it be identified from the a priori estimates inverse analysis. However, &in < S?l < SmpXcan be experimentally provided. For typical data of the experiment performed in this study, the value of the dimensionless parameter

1

(28f)

where Ax=;,

749

At=$

and 0 is taken to be 2/3, which is the Galerkin weighting. In fact, it can be proved [20] that the above FDM is unconditionally valid, i.e. stable and convergent, for 8 E [l/2, 11.At each time level ti+ I forj = 0, (N - l), the system of linear equations (28) is solved using a Gaussian elimination method. Logarithmic finite-difference grids can be more illustrative, but in this study for simplicity the linear grids were considered sufficiently scaled for the graphical presentations, discussions and conclusions of the numerical results.

SZAL

B=T
results in being of O(10V4) which is too small for identifying the specific storage, see Neuzil et al. [15], Timmer [16], and Wang and Hart [21], who recommended 0.01 < fl < 10. Therefore, in the present investigation the inverse analysis requires, whilst assuming available prior estimates for St, only the determination of the pressure and the hydraulic conductivity, i.e. of the pair (p(x, t), K(x)), which satisfy the problem given by equations (18)--(21). Although uniqueness conditions for this inverse problem are difficult and not yet stated in the literature, it is to be expected that strong additional information should be provided in order to render a feasible unique solution. In the first instance, physical constraints can be imposed, an obvious one being K(x) > 0,

5. SENSlTIVITV COEFFICIENTS

Prior to performing the inverse analysis of identifying the hydraulic properties of the rock, it is useful to calculate the sensitivity coefficients as a function of time. Sensitivity coefficients are the first derivatives of the measured quantities, i.e. pressure and hydraulic flux, with respect to the unknowns, i.e. the hydraulic properties S, and K. They provide indicators of how well-designed is the experiment. In general, the sensitivity coefficients are desired to be uncorrelated, i.e. linearly independent. A sense of the magnitude of the sensitivity coefficients is gained through normalising them with respect to their corresponding unknown differentiation variable, resulting in units of pressure or hydraulic flux for the normalised sensitivity coefficients. The degree of uncorrelation of these coefficients can then be illustrated by the departure of their ratios from a constant function. Based on this criterion we can determine the optimal time instant measurements to be imposed or recorded in order to reduce the ill-posedness of the inverse formulated problem. Therefore, a study of the sensitivity coefficients, prior to performing experiments, can lead to better experimental designs. Sensitivity coefficients for the hydraulic experiment under investigation, modelled mathematically in Section 3, will be further presented in Section 7.

S,,,axAL

x E [O, 11.

(31)

Further, a first attempt is undertaken to analyse the particular case when unknown hydraulic conductivity K(x) belongs to a finite-dimensional space of known dimension, e.g. the space of quadratic polynomials Pz, i.e. P* = {K(x) = Ko + K,x + W}.

(32)

Further, this assumption may generalise for identifying spacewise dependent hydraulic conductivities K(x) which can be approximated by piecewise quadratic interpolation functions. Cubic splines may also be used, but then this increases the dimension of the finite-dimensional space in which the solution is sought and for our inverse formulated problem this approximation was found inappropriate, i.e. too large to retrieve all the unknowns from the very limited boundary measurement information. Alternatively, if the piecewise interpolation approach is considered a too strong constraint for the solution space, then one may use these preliminary identification results as initial guesses in a more general least-squares minimisation problem [l I]. Moreover, it is noted that if the dimension of the solution space is taken as a further unknown in the problem or is infinite, then the inverse identification problem should be reformulated to allow for the pressure p(x, r) to be measurable and known, probably almost everywhere in both space

750

LESNIC et al.:

HYDRAULIC

and time, which is not practically possible or realistic at present. Further, we agree that the number of unknowns should be kept as small as possible and the whole of the available measurement information be used. However, this approach may affect the existence or stability of the solution in the sense that too much or possibly redundant information is added to the inverse problem formulation. Therefore, our approach is simple but rational, by starting with one to three unknowns to be retrieved from the same number of time measurements available, which is the minimal necessary condition for identifiability. Further some preliminary numerical investigations have shown that an increase in the number of unknowns affected the uniqueness of the solution, whilst an increase in the time measurement information did not significantly improve the results and in fact, they may affect the stability of the solution. 6.1. Steady-state mathematical analysis

A simple steady-state analysis of the problem given by equations (23) and (24) shows that the pressure distribution inside the sample at the steady-state time, t + co, is given by (33)

&ln p(x? 60)=AK,

(34)

P(X, co) =

(x - ~‘)(L - ‘*) if* > (x - L*)(L - 2,)>

&In

Qu(L - xl

if A = ’

[tan-i(y)-tan-‘(F)]

if A<0 (35)

K(x) = K,,, i.e. the rock is homogeneous, K(x) = KO+ kx, i.e. the rock is linearly heterogeneous and when K(x) = Ko + kx + Kzx2, i.e. the rock is

when

quadratically heterogeneous, respectively. In equations (33x35) the quantities involved are defined by A

=

I

-K,

- ALI2 )

(38)

Pa 00) = &ln(t:[iI

1:;)

AK,G(“L+ S) 44’2 ,”

ifA > 0

if A = ’

L [ tan-’

(E)-tan-l(i)] c

ifA < 0.

(39) Moreover equation (39) is difficult to use as additional information in the inverse problem, since, in general we do not known a priori the sign of A. However, the pressure measurement given by equation (37) may be used in practice for determining the homogeneous hydraulic conductivity K,, if the steady-state time is not very long. Therefore, the practical use of equations (37x39) is limited if: (i) long steady-state times are encountered; (ii) internal measurements are not available; (iii) heterogeneous samples are tested. Whilst in the present investigation case (i) is avoided by the design of the experiment, cases (ii) and (iii) are very significant due to their practical importance. Case (ii) can alternatively be replaced by recording, instead, transient boundary data, whilst case (iii) can be considered by solving the inverse problem given by equations (18x2 1) numerically, as described in the next section.

c=

In contrast with a similar formulation by Huang and Ozisik [ 1l] of an inverse problem in heat conduction, where reasonably accurate (temperature) readings are in general available to be taken with sensors installed at various locations within the sample, in the present hydraulic experiment, as argued in Section 1, only boundary measurements are available without damaging the sample. Therefore, in practice, possible reliable measurements will involve the temporal values of the pressure at the upstream face of the sample x = 0, i.e. ~(0, t) and the values of the hydraulic fluxes at both ends of the sample x E (0, 1}, i.e. qO(t) = q(0, t) = K(O)(ap/ 8x)(0, t) and q,(t) = qU, t) = K(lWpldx)(l,

2K2

K, 6 =%,

OF ROCKS

6.2. Inverse numerical solution

AKz(L + 6)(x + 6)

” &-A)“*

0

CONDUCTIVITY

(-A)“* 2K2 ,

A=K:-4K,,Kz.

(36)

In this study, only boundary measurements are considered and, therefore, equations (33x35) can provide only a single equation at x = 0 for the unknowns, respectively, pa

00) =

z

(37)

t>.

In a first numerical attempt, based on an ordinary least-squares approach, we consider the least-squares functional, i.e.

WK; 1’) = i 1; @(O, t;)‘4 - p(o, t;)(e))* i= I

+ i i=l

ap Gl K(O);i;; (0, t;)

K(O)zap(0,ti’)

LESNIC

+i

i-l

G

CL

K(1)!2 ax (13tl!)

WI

-

1

[

et al.:

HYDRAULIC

K(1)zap(1, tl)

I) w

RESULTS AND DISCUSSION

We first examine the suitability of the mathematical formulation described in Section 3 to model the practical experiment investigated in this study. In this hydraulic experiment the pressure increase at the upstream face of the sample, x = 0, is measured as a function of time until usually the steady-state is reached. After the equilibrium is achieved the tap which supplies the inflow mass of fluid is closed and, therefore, a pressure decay curve is recorded at the upstream face of the sample, x = 0. The numerical results for the pressure increase curves at x = 0 are reproduced for various constant inflow rates Q”, whilst the pressure decay curves at x = 0, when QU= 0, are obtained under the “initial-decay” condition P(% L) = POW,

x E LO,11

OF ROCKS

751

Cd 2

where A;, A&,and ni are either zero or unity depending on whether we measure the multiplying quantities in equation (40) or not, the symbol f denotes whether the pump is on, off, or on and off, Z is the total number of pressure or hydraulic flux measurements taken at times t’ = (t:) for i = m, the symbol (e) simulates the input measured data and the symbol (d) simulates the FDM numerically obtained values from the direct solution of the problem, as described in Section 4, by using iteratively the estimated values of the unknown hydraulic conductivity vector K = (Ko, K,, K2). Then, in order to identify the vector K, the nonlinear least-squares functional LS given by equation (40) is minimised. An iterative NAG routine, i.e. E04UCF [22] based on quadratic programming, a line-search with an augmented Lagrangian merit function and a quasiNewton update of the approximate Hessian of the Lagrangian function, is utilised for this purpose. In the present situation, this routine minimises the functional LS on the space P2 defined by equation (32). If a sufficient condition for satisfying equation (31) is adopted, bounds are specified for the parameters Ko, K, and K2 between zero and a large positive number. Further, if the sample is assumed homogeneous then & > 0 and K, = K2 = 0, if it is linearly heterogeneous then K,,, K, > 0 and Kz = 0 and finally, if it is quadratically heterogeneous then KO, Kl, Kz > 0. The gradient of the functional LS, the boundary hydraulic flux and the sensitivity coefficients are calculated using forward or backward finite differences. 7. NUMERICAL

CONDUCTIVITY

(41)

where pa(x) is the pressure distribution inside the sample at t = t, obtained before the pump is switched off suddenly. Since an analytical solution, which is similar to that derived by Hsieh et al. [ 191based on Laplace transforms, does not seem available at the present, the discretisation of the inverse problem under investigation has been performed numerically using M = 35 and N = 100,

which corresponds to a dimensional space step of AX = 0.097 cm and a time step of At = 200 set, respectively. A finer mesh discretisation did not significantly affect the accuracy of the numerical results, but only increased the computational time. 7.1. Direct solution In the first set of experiments which are used to validate the mathematical model described in Section 3, small, medium and large yellow bunters (samples) have been experimented upon subject to various inflow rates. Figure 3 shows the numerical and experimental results for the pressure increase and decay, as a function of time, at the upstream face of the sample, x = 0. The sample has length L = 3 cm, cross-sectional area A = 22.6 cm2 and measured mean value of hydraulic conductivity given by K* = 1.35 x 10e6 cm/set. The specific storage of the sample could not be measured exactly, but it should vary between 5 x lo-” cm-’ = S,,,,, < SF < SmaX = low9cm-‘. However, the numerical results have been found insensitive to both any variation of the specific storage between the aforementioned limits and any spacewise dependency, as also discussed in Sections 1 and 6. Based on this numerical investigation, in the remainder of this paper the value S? = 10m9cm-’ was chosen. In Fig. 3 the constant inflow rates which reservoir are are supplied to the upstream QUE {1.6, 3.2,4.8) x 10m3cm3/sec, which give rise to the steady-state times and pressures, i.e. t, E {40,43,48) set and pm E { 16,32,48) kPa, respectively. If the rock is homogeneous, i.e. K(x) = K*, then these pressure equilibrium values can be found from equation (37). The value of the compressive storage of the upstream reservoir was difficult to measure, but it was found that an “adjusted” value of S, = 0.75 x 10m4cm2 produced good agreement between the numerical and experimental results in all the results presented in Figs 3 and 4, although more investigations need to be undertaken in the future. In Fig. 4 the length of the sample is increased to L = 5.6 cm, whilst maintaining the same cross-sectional area and then the measured hydraulic conductivity is given by K* = 2.1 x 10e6 cm/set. In Fig. 4 the inflow rates are QUE {1.6,&O} x 10e3 cm’/sec, which gives rise to the steady-state times and pressures t, E (71,92} set and pm E { 19.5,98} kPa, respectively. Again, the same good agreement between the numerical and experimental results is obtained. Although not further illustrated, the same conclusion was also obtained when the length of the sample was further increased to L = 9.25 cm. In the second experimental test the sample has a low hydraulic conductivity of mean average K* = 1.48 x 10d8cm/set and a very low constant inflow rate of QU= 3.33 x IO-’ cm3/sec is supplied at the upstream reservoir. The length of the sample is L = 3.42 cm, with a cross-sectional area A = 7.66 cm2. The pump was switched off after t, x 20,000 set, with a corresponding pressure value of pm 2 96 kPa. However, if the material was homogeneous, then equation (37) gave the equilibrium pressure as 100.5 kPa. This shows that the

152

LESNIC et al.:

HYDRAULIC

pump might have been switched off slightly earlier before all the transients have completely passed, or possibly that the material could be slightly heterogeneous. At this stage it is worth noting that for this experiment, the compressive storage of the upstream reservoir was measured to be S, z 2 x 10m4cm*, but due to the difficulties encountered in performing this measurement we tend to believe that a more realistic range is S, E (1.4,2.6) x low4cm’. In addition, it has been found that the pressure increase and decay curves are quite sensitive to the value of S,. This can be seen in Fig. 5 where the numerical results for ~(0, t) are compared with the experimental data for various values of SE 11.4, 1.7,2.0,2.3,2.6} x 10e4cm2. The pressure curves are insensitive to the values of the specific storage X E (Smin,S,,,,,) because fi = 0(10-4), see equation (30), as argued in Section 6, but they are sensitive to S, because of the presence of the non-zero dimensionless coefficient yu =

$f” CC= O(1)

see equation (22), in the boundary condition (21). This conclusion is somewhat in contrast to that obtained by Neuzil et al. [15] from the theory of Hsieh et af. [19] in which QU= 0 and, hence yu = 0, but this is because the boundary condition (21) allows the presence of the coefficient yu of O(1) which is important. In fact, the consistency of our results with the previous pulse theory in which QU= 0 can be observed from the pressure decay graphs in Fig. 5, where at the beginning of the graphs when QU= 0 and hence yu = 0, the curves look insensitive to the value of S,. However, this conclusion cannot be extended for the whole of the pressure-decay curve due to the non-zero initial-decay condition (41). Unfortunately, the presence of the non-zero coefficient y,, in equation (21) prevented a similarly based Laplace transform analytical solution to that obtained by Hsieh et al. [19] from being derived and, therefore, at present our arguments are only numerically based. In addition, this study also considers heterogeneous samples for which any analytical solution for (1 I?),regardless of what type of boundary condition is imposed, is unlikely to be obtained using traditional methods. From Fig. 5 it can be seen that the best agreement between the numerical results and the experimental data for the pressure increase and decay are obtained for the values S, = 2.6 x 10V4 and 2.0 x 10m4cm2, respectively. Furvalue measured only the experimentally ther, S, = 2 x 10e4 cm2 will be considered in the inverse analysis, but it might be that a correction to the term S, which is present in the boundary condition (21) should be made and taken into consideration in future work. This correction may involve assuming that the compressive storage S, is time dependent, or has different values when the pump is on and off. Nevertheless, the overall agreement between the numerical and experimental results obtained in Figs 3-5 which have been obtained over a wide range of narameter selection and technical data. confirms the

CONDUCTIVITY

OF ROCKS

reasonable validity of the mathematical formulation presented in Section 3 for modelling the pump-flow experiments under investigation. 7.2. Sensitivity analysis Once the mathematical model has been validated, prior to performing the inverse analysis it is useful to undertake a sensitivity analysis of the input boundary data with respect to the type of heterogeneity of the rock. Figure 6 shows the numerically obtained boundary pressure and hydraulic fluxes curves for constant, linear and quadratic hydraulic conductivities which have the same data and mean average conductivity of the second hydraulic experiment of K* = 1.48 x lo-* cm/ sec. From Fig. 6 it can be seen that the curves are consistent with the steady-state analysis given by equations (37)-(39) as t + 00, and that they are quite sensitive to various sample heterogeneities, thus confirming that the inverse analysis which is undertaken is meaningful. A more rigorous sensitivity analysis can be performed by calculating, as described next, the ratios of the sensitivity coefficients (briefly discussed in Section 5). The (normalised) sensitivity coefficients, as function of time t, are calculated using forward finite differences, i.e. Sens(T; Ki) = Ki g.

I

T(K + h) - T(K) , i = h

>

0

1

2

5 >

(43)

where h = lo-l3 cm/set is a small perturbation of conductivity coefficients the hydraulic and T E (~(0, 0,40(t), 41(t)}. A much enhanced characterisation of the degree of correlation of the sensitivity coefficients can be viewed by calculating the ratios between the sensitivity coefficients, i.e. R)(T) = Sens( T; Ki)/Sens( T; Ki + I),

i = 0, 1,2

(44)

with the obvious circulating convention. The ratios of the sensitivity coefficients will be presented for linear and quadratic spatially dependent hydraulic conductivities, together with their implications, in Sections 7.5 and 7.6. 7.3. Inverse analysis The inverse analysis is applied as described in Section 6.2. With the technical data of the second experimental test, the FDM, combined with the least-squares minimisation technique, is performed for the identification of constant, linear and quadratic spatial dependences of the hydraulic conductivity K(x). The initial guess for all the cases considered was taken to be K = lo-‘* cm/set. In order to investigate the most general case of identifiabilitv of the narameters K. the number of time

LESNIC et al.:

measurement I is taken to be one to constant, linear or quadratic hydraulic respectively. Further, the selection of the t; for i = 17, in the least-squares functional given by equation (40) is following optimal criterion min,,= I.((N-~)/~(ll((‘-) KQ’)I; and

HYDRAULIC

three for the conductivities, sampling time minimisation based on the

LS(K(t’), t’) = min&S(K, r’)}

RMMS M,5--D

Constant

Linear

93 39 5 93

19 19 19 43 22 14

9:

753

Table 2. The retrieval of a homogeneous hydraulic conductivity obtained using the least-squares FDM from only boundary experimental measurements, when the pump is on (+), off (-), or on and off (+1 FDM

Kox6E8

LS(Kb)

8.880 9.702 8.890 9.436 9.136 8.803 8.982 9.141 8.803 8.982 9.310 8.814 9.106 9.310 8.819 9.106 9.139 8.799 8.982 9.241 8.810 9.054

0.000

(45)

Table 1. The optimal valuts of iT+(pump on) and if (pump off) for constant K(x) = KO= I .48 x lo-” cm/set, linear K(x) = Ki + K,x = (1.11 + 0.21x) x lo-* cm/~ and quadratic K(x) = Kb+ Klx + hydraulic conducKG’ = (0.21 + 0.37x + 0.16x2) x lo-* cm/see tivities

._ ‘P G .‘,, G ._ ‘#I

OF ROCKS

where I’ = (tjXrT)j-n

where R’) is the exact or an u priori estimation of the hydraulic conductivity parameterised vector which is sought to be recovered and i;t or i; denotes the minimum point in equation (45) when the pump is on or off, respectively, and T E (~(0, t), qo(t), q,(t)). It is also worth noting that in the computational program the value of tj X(7 is understood as (tj XiT + tm). Although in this study use has been made of the available exact values of the hydraulic conductivity coefficients for estimating P) in equation (49, a priori estimates can be provided in the homogeneous case by using any value greater than the experimental measured equilibrium pressure ~(0, 00) in equation (37) [23], whilst in the heterogeneous cases these prior estimates are assumed to be experimentally available or to be numerically simulated from a hypothetical hydraulic conductivity within a confidence range. The optimal criterion given by equation (45) was chosen taking into account that sufficiently large sampling time, at equidistant finite-difference time node locations, and with no measurement being imposed at the instants when the pump is switched on or off and its control becomes difficult, are preferable to the experimentalist. In other circumstances, the sensitivity analysis (Section 7.2) would offer optimal numbers and locations of the time measurements to be recorded and imposed in equation (40). The optimal values of if and i; , given from equation (45), when exact estimates for K(e)are assumed available, are summarised in Table 1 for constant, linear and quadratic K(x), and for the on (+) and off (-) operations of the pump. These optimal values correspond to individual measurements of either pressure or hydraulic fluxes, and they are also considered when combined measurements of both pressure and hydraulic fluxes are imposed in equation (40). The corresponding values of the optimal sampling times t’ will then be used in the following sections.

ip+

CONDUCTIVITY

Quadratic 15 15 11 19 :z

lE7ElE+3 2E2E9ElEZE9ElE2E4EIEZE4E8E7ElElE3E5E-

I5 14 16 16 13 18 16 13 12 14 12 12 14 12 17 19 12 12 14 12

Iter

8* 6* 6*

I 3 4 3 3 4 3 4 4 3 5 4

1 4 4 3 5 4

With the above notations the least-squares functional given by equation (40) becomes

+ i

i=,

[Agf,(qO(ti i&)‘4- qO(ti i&)(c))2 X

x

7.4. Constant K(x) = Ko For homogeneous rocks, i.e. K(x) = Ko = constant, the experimental data from Fig. 5 for S, = 2 x lo-’ cm* is input for the boundary pressure ~(0, t’)@) taken at the optimal times t;+ = rg = 18,600 set and/or t;- = f,; + t, = 27,800 set (Table l), in equation (40). As yet there is no experimental curve available for the hydraulic fluxes so these were numerically simulated by adding to the numerically calculated values from the FDM solution of the direct problem (Section 4), 5% random Gaussian normal noisy variables. It should be noted that this amount of noise is comparable to the errors observed in Fig. 5 between the experimental and numerical data for the pressure increase and decay curves obtained for S, = 2.6 x lo-‘cm2 and 2.0 x lo-‘cm2, respectively. Table 2 shows the recovery of the homogeneous hydraulic conductivity, &, obtained using the leastsquares FDM from only boundary pressure and/or hydraulic flux measurements, in comparison with the exact value & = K* = 1.48 x lo-* cm/set. When the least-squares functional given by equation (46) contains both boundary pressure and flux measurements, then a factor of 10” has been included for scaling. From Fig. 5 and Table 2 it can be seen that the experimental pressure

154

LESNIC

et al.:

HYDRAULIC

measurements at x = 0 seem to possess too much experimental and numerical noise, and the NAG routine E04UFC used indicates (in Tables 2-6) that the current point cannot be improved because during the final line search a sufficient decrease in the merit function is unattainable. However, this failure can be overcome by recording additional hydraulic flux measurements and then an optimal solution is found in less than four-five iterations. The experimental pressure measurements at x = 0, when the pump is switched off, i.e. 1; = 1, are more accurate in comparison with those encountered when the pump is on (Fig. 5) and, therefore, it is possible to reasonably estimate the exact value of &. This conclusion is also valid for the hydraulic flux, and data provided by the off pump measurements possess more information and, therefore, a better retrieval of & than the on pump measurements. When on and off pump measurements are gathered, the conclusion traced from Table 2 is that the single on pump measurements are better interpreted, whilst the single off measurements are worst interpreted as more noisy data is added to the inverse procedure. Finally, from Table 2 it can be concluded that the single off pump measurement of either boundary pressure or hydraulic fluxes, when a single optimal time measurement, i.e. I= 1 is allowed, contains the best information such that when interpreted, i.e. imposed in equation (40), it provides a good recovery of the exact value of &. In the absence of such a thorough numerical investigation which allows minimal available input data, from the practical point of view, Table 2 also shows (see 10”~~ = 1; = A,‘,= 1) that gathering all the available measurement information still provides a reasonable estimation of the homogeneous hydraulic conductivity of the rock without the peril that significantly, unwanted, redundant data are added to the ill-conditioned, inverse problem. Overall, from Table 2 it can be concluded that the FDM least-squares approach in general retrieved successfully the homogeneous hydraulic conductivity of the rock subject to the hydraulic experiment under investigation, with an accuracy comparable to the errors encountered by the experimental and noisy simulated input boundary data.

CONDUCTIVITY

OF ROCKS

Table 3. The retrieval of a linear heterogeneous hydraulic conductivity obtained using the least-squares FDM from only boundary errorless measurements, when the pump is on (+), off (-), or on and off ( j-) FDM

1015~* P = /y. = A*#I = 1

6.660 3.176 3.176 3.776 6.645 6.596 6.657 4.844 6.794 6.546 6.632 6.608 6.660 4.783 6.79 1 6.599 6.659 6.662 6.659 6.661 6.660 6.659

1.298 3.776 3.776 3.776 1.309 1.341 1.300 2.681 1.200 1.377 1.319 1.333 1.298 2.735 1.203 1.340 1.298 1.296 1.298 1.296 1.297 1.298

LS(Ko, K,) 0.000 3E-6 3E- 18 3E-5 IE17 3E- 17 lE17 3E14 7E16 6E15 6E17 7E- 17 lE18 3E- 14 YE- 16 IE15 3E-21 2E18 2E18 IE- 17 5E- 19 8E-20

Iter 0

7* 4 49* 20 21 11 7 11 9 20 19 12 7 12 12 12 14 8 11 14 7

and K, obtained using the least-squares FDM from boundary pressure and/or hydraulic flux exact or 1% noisy measurements, respectively. From Tables 3 and 4 it can be seen that only boundary pressure measurements do not contain sufficient information for the identification of both K. and K,. Also, although from Table 3 exact flux measurements at x = 1 only appear to improve on the identifiability of K,, and K, , from Table 4 the same measurements when subjected to noise produce an unstable solution. However, Tables 3 and 4 show that the flux measurements at x = 0 only produce an identifiable and stable estimation of both K. and K,. These conclusions may also be argued from Fig. 7 in which the ratios of the sensitivity coefficients, as a function of time, are plotted. Clearly, the sensitivity

Table 4. The retrieval of a linear heterogeneous hydraulic conductivity obtained using the least-squares FDM from only boundary 1% noisy measurements, when the pump is on (+), off (-), or on and off (+) FDM

7.5. Linear K(x) = K0 + K,x For a heterogeneous linear hydraulic conductivity, the values of K. and K, to be retrieved are chosen so as to satisfy equation (16) with the mean average hydraulic conductivity measurement supposition, i.e. K* = K; = 1.48 x lo-* cm/set, e.g.

The boundary values are inputted as exact or noisy measurements by adding 1% random variables to the FDM numerically calculated values at the optimal times 2; , i = 1, I, Z = 2 (Table 1) from the solution of the direct problem (Section 4). Tables 3 and 4 show the recovery of the linear heterogeneous hydraulic conductivity coefficients, K.

& x 6E8 K, x 6E8

,0’51*P = A*a = A*q, = 1

KO x 6E8 K, x 6E8 6.660 3.731 3.763 3.759 6.879 6.684 7.631 4.798 5.541 3.078 7.924 5.535 7.646 4.724 4.297 3.801 6.690 5.900 6.742 6.694 5.905 6.748

1.298 3.731 3.763 3.759 1.188 1.259 0.647 2.789 2.102 4.591 0.449 2.082 0.636 2.854 3.227 3.746 1.334 1.809 1.262 1.331 1.807 1.256

LS(KO, K,) 0.000 lE+O 2E- 1 lE+ 1 2E- 14 5E- 15 2E- I8 5E15 5E- 14 7E- 14 lE14 8E- 14 3E- 14 3E- 14 6E15 3E- 14 lE13 IE13 3E- 12 IE- 13 lE13 3E-12

Iter 0 7* 491 50* 17 17 14 14 43 25 9 I2 8 23 16 11 12 8 11 12 8

LESNIC el al.:

HYDRAULIC

CONDUCTIVITY

OF ROCKS

155

Table 5. The retrieval of a quadratic heterogeneous hydraulic conductivity obtained using the least-squares FDM from only boundary errorless measurements, when the pump is on (+), off (-), or on and off (k) FDM

I

I

1*104

I

2*104

3’104

I 4*104

time (seconds) Fig. 7. The ratios of the sensitivity coefficients, i.e. (-) h(p), (---) &&JO) and (. .) R&I), for a linear heterogeneous hydraulic conductivity.

coefficients for the boundary pressure are strongly correlated giving Z&,(p)= l/Z&(p) x 3.6 and, therefore, we cannot identify both K, and k from this single measurement. Also, the sensitivity coefficients for the flux q,(t) appear slightly correlated, whilst for qo(o(t)they are uncorrelated and therefore the identifiability of both Ko and K, is expected. Finally, from Tables 3 and 4 it can be seen that the stable estimation of & and ZG can also be achieved when the hydraulic flux measurements at x = 0 are gathered with the hydraulic flux measurements at x = 1, or with flux measurements at x = 1 and pressure measurements at x = 0. The overall conclusion from Tables 3 and 4, and argued by Fig. 7 is that the hydraulic flux measurements at the hydraulic contact boundary x = 0 are necessary for a stable identification of the linear heterogeneous hydraulic conductivity of the rock.

&x6E8

Exact ;r,’= 1 A; = 1 if = 1 &= 1 I., = 1 n,t, = 1 n,: = 1 & = I n,: = 1 lO”l+ = A+ = 1 lOlJi[ = 19 = 1 10151;= A? = 1 lO15iC= $! = , = = lO1SA$ *0+ = 19 AP = 1 n,i, &;, =“I, n, = 1, = 1 1’ = 1’ = 1 &+Z~+=~+=l lOl$=$$=l lO$=&$=l D 0” 0,

K, x6E8

Kzx6E8

1.268 1.421 1.421 1.421 1.291 1.272 1.080 1.416 0.971 1.355

2.225 1.421 1.421 1.421 2.248 2.234 1.801 1.495 2.806 2.077

0.976 1.421 1.421 1.421 0.911 0.961 1.803 1.346 1.059 0.955

I.166 I.297 1.140 1.305

2.287 1.998 1.939 2.145

0.873 1.389 1.506 0.980

1.137 1.128 1.344 1.268 1.268 1.344 1.269 1.268

2.457 2.481 2.459 2.224 2.225 2.458 2.221 2.225

1.017 1.013 0.670 0.917 0.976 0.671 0.973 0.976

LS(ki,G,Kz) 0.000 5E- 10 SE- 10 5E-9 9E- 17 2E- 18 2E- 18 4E- 15 2E- 17 2E-20 3E- 19 lE19 9E- 19 8E- 18 lE19 6E- 17 4E15 3E- 18 lE19 6E- I5

IE-

17

lE-

18

Iter 0 4* 4* 4* 9 9 11 2 7 9

10 9

11 15 8 6 10 15 18 9 12 21

and flux ql(t) are correlated as their ratios are constant, or, almost constant functions of time, i.e. &l(p) x Ro(q1) = 1.15,

R,(p) = R,(q,) = 2,

R*(p) x &(q,)

x 0.43

(49)

whilst the sensitivity coefficients for qO(t) are uncorrelated. Table 6 also shows that the hydraulic flux noisy measurements at x = 0 only provide a stable solution, whilst, unlike in the linear case in Table 4, gathering all the available boundary noisy measurement information results in too much unwanted, noisy, redundant data being contained in either ~(0, t) or q,(t), so destroying the stability of the inverse solution. Therefore, from the

7.6. Quadraric K(x) = K, + K,x + Kzx* For a quadratic heterogeneous hydraulic conductivity, the values of &, K, and Kz to be retrieved are chosen so the mean average of K(x) is K* = ZG = 1.48 x 10m8cm/ set, e.g.

&q-,

j&g

&=S.

(48)

Tables 5 and 6 show the recovery of the quadratic heterogeneous hydraulic conductivity coefficients, K,, K, and K2 obtained using the least-squares FDM from boundary pressure and/or hy&aulic flux measurements at the optimal times tl, i = l,Z, Z = 3, respectively. The conclusion for the linear case from Table 3 is also validated for exact data for the quadratic case in Table 5, i.e. that only the hydraulic flux measurements at x = 0, i.e. q&‘), contain the necessary information for the identifiability of all the coefficients I&,, K, and &. This conclusion is again motivated on the basis of the ratios of the sensitivity coefficients which are plotted, as a function of time, in Fig. 8. This figure shows that the sensitivity coefficients for the boundary pressure ~(0, t)

Table 6. The retrieval of a quadratic heterogeneous hydraulic conductivity obtained using the least-squares FDM from only boundary 1% noisy measurements, when the pump is on (+), off (-), or on and off (+) FDM Exact n; = 1 n; = I &+ = 1 n; = 1 A; = l n& = 1 1; = 1 1, = 1 2,: = 1 lO”1+ = A+ = 1 lO’$- = 12 = 1 ,0’5$ = A? = 1 ,O$+ = $ = l = = lOoA: lOl$- = AP A’-’ = 1 A&CPA; f’l A,-,= 1, = 1 n& = A; = 1 10’52+=L+=I+=l 10rs~~=kn’ P 10-“=n~=l “I 10’51’ D =Afa.,=,q I = 1

Kox 6E8 Kt x 6E8 KZ x 6E8 LS(Ko,Ka,K2) 1.268 1.398 1.415 1.414 1.302 1.295 1.184 1.425 1.319 1.641 1.333 1.125 1.251 1.678

2.225 1.398 1.415 1.414 2.128 2.115 1.884 1.508 0.571 0.082 1.912 1.773 2.020 0.077

0.976 1.398 1.415 1.414 1.000 1.025 1.460 1.362 2.511 2.295 1.076 1.698 1.186 2.228

1.447 1.695 1.158 0.865 0.865 0.774 0.922 0.866

0.556 0.052 1.936 1.145 1.144 0.884 1.298 1.147

2.202 2.324 1.501 3.439 3.442 4.376 2.943 3.428

0.000 2E+O 2E+O 2E+ 1 lE18 4E- 18 8E-20

lElE-

14 15

6E-20 lE3E-

14 14

IE5E3EIEIE3E7E4E4ElE-

13 14 14 13 14 18 18 14 14 13

Iter 0

4* 4* 4* 8 9 10 2 7 13 8 9 12 16 9 13 28* 34 40 Ill 31 35

LESNIC et al.:

HYDRAULIC

(a) 2.0

______..____.._____.______.__

1.6

(1) The boundary pressure is easy to measure very accurately using electron transducers. (2) A very low flow rate is much easier to enforce using a pump than some outflow volume change measuring device that must take some average, rather than instantaneous, value. (3) The facility of switching the pump on and off during the experiment enables both pressure increase and decay data to be recorded, which adds more information into the formulation of the inverse problem.

1.2

1

0.8

i

0.6 i

(b) 30

In Section 3, the pump-flow experiment is modelled and the validity of the mathematical model is highlighted in Section 7.1. In comparison with other models for transient pulse-decay or pump-flow tests, the model given by equations (18)-(21) possesses some realistic feature generalisations, e.g.

1

-30

] 0’100

I 1*104

I

I

I 2.104

3*104

4*104

(cl 2.0

. . . . . ..__..................._........

_::,_______

_______

_______

___.

_

,I

I

i

1.5

OF ROCKS

modelled. The advantages of this pump-flow experiment over other existing methods when testing low permeable rocks are highlighted in Section 1. In summary, these advantages are:

_ . . . . . . . . . . . . . . . . . . . . . . . . . . -----...

1

I.0

CONDUCTIVITY

1

(1) It is mathematically derived, whilst the boundary condition given by Lin [ 181 is argued on the basis of a private communication. In addition, it incorporates the enforcing of an inflow rate at the entrance of a reservoir. (2) The compressive storage of the upstream reservoir, which is maintained in hydraulic contact with the rock, is not neglected, as was unrealistically assumed by Morin and Olsen [3]. (3) The rock may possess heterogeneous hydraulic properties, which is always true in reality, whilst Hsieh et al. [19] considered only the homogeneous case for which an analytical solution was readily available. The generality of the mathematical model formulated as a direct problem prevented any analytical solution and, then, a FDM was numerically employed in Section 4. The FDM results compared well with the experimental results (Figs 3-5) for homogeneous hydraulic conductivities. In Section 5, the importance of displaying the non-constant ratios of the sensitivity coefficients as a characterisation of the degree of uncorrelation of the hydraulic property coefficients that can lead to better experimental design, has been highlighted. The sensitivity analysis has been performed prior to the inversion in Section 7.2, and shown in Figs 7 and 8, for the argumentation on the consistency of the numerical results of the inverse identification problem. At this stage it should be noted that if, in a future investigation, too many hydraulic property coefficients are to be retrieved, then the criterion for the ratios of the sensitivity coefficients not to be a constant function of time should be replaced by the maximisation of the determinant of the covariance matrix of the sensitivity coefficients or its condition number [24]. In Section 6, the pump-flow model problem is inversely formulated as the problem of identifying the spacewise dependent hydraulic conductivity of the rock from only boundary measurements. Unfortunately, as

:)1....... ‘,

0.100

1*104

,

2*104

(

3*104

(

4.104

time (seconds)

Fig. 8. The ratios of the sensitivity coefficients, i.e. (-) Ro(T), (- - -) R,(T) and (. .) RI(T), for: (a) the boundary pressure T=p(O, I); (b) the hydraulic flux T= q&); (c) the hydraulic flux T= 91(t); for a quadratic heterogeneous hydraulic conductivity.

investigation performed in Tables 5 and 6 we can conclude that the hydraulic flux measurements at the contact boundary x = 0 are necessary and sufficient for a stable identification of the quadratic heterogeneous hydraulic conductivity of the rock. 8. CONCLUSIONS AND FUTURE WORK

In this paper, a hydraulic pump-flow permeability test (Section 2) is both directly and inversely mathematically

LESNIC et al.: HYDRAULIC CONDUCTIVITY OF ROCKS

previously illustrated by Neuzil et al. [15] and Trimer [ 161, for small samples or rigid rock specimens investigated in this study, from boundary measurements only, the specific storage of the sample cannot be numerically identified, but experimentally can be measured within some confidence range and therefore it was assumed to be known. Unlike in the pulse tests, in the pump-flow experiment of this study the pressure curves were found to be quite sensitive to the values of the compressive storage of the upstream reservoir (Fig. 5), as argued in Section 7.1. Further, in contrast with other inverse formulations which often lack physical motivation and use internal measurements, e.g. Huang and Ozisik [I 11, it should be noted that only boundary pressure and/or hydraulic flux are used as available measurements in the inverse model since we want to maintain the generality of the problem, and to take into account both laboratory and real life difficulties of undertaking internal measurements, e.g. (1) Internal measurements in our laboratory experiments would require a special core holder with pressure ports leading to certain points on the outside of the sample to be available. (2) In real petroleum reservoirs, internal measurements would require new wells to be drilled, each costing several million pounds. Once the direct mathematical model has been validated by the reasonable comparisons shown in Figs 3-5, the inverse analysis of identifying the hydraulic conductivity from only boundary pressure and/or hydraulic flux measurements is performed in Sections 7.4-7.6. In the first instance, the hydraulic conductivity has been assumed to be a constant, linear or quadratic polynomial function in the space variable x, but future work will be concerned with more general spatial dependences of K(x). Since the steady-state mathematical analysis given by equations (33x35) has practical limitations as argued in Section 6.1, the coefficients of these polynomial functions were numerically recovered by using an ordinary least-squares minimisation as given by equation (4O), subject to positivity constraints on the unknowns. The results obtained are shown in Tables 26. A sophisticated iterative NAG routine, E04UCF, is used for solving this constrained minimisation problem, whilst allowing arbitrary initial guesses and good control for finding an optimal solution. As argued in Section 6, in the inverse formulated problem the number of time measurements equalled the number of unknowns in order to numerically investigate, for the first time, the most general case of parameter identifiability under various modes of on, off, and on and off operations of the flow pump. In those circumstances a new criterion for selecting the optimal time boundary measurements of the pressure or hydraulic fluxes, equation (49, was introduced. At this stage it is fair to say that this optimal criterion assumes the availability of exact or prior estimates on the unknowns, but this is also the situation investigated by Carrera and Neuman [25] who showed that using prior estimates improves on the globality of

151

the minimum of a modified least-squares functional. However, in this study the prior estimates are used in equation (45) only for determining optimal time instants at which boundary measurements should be taken and, in fact, this is the procedure which everyone employs when undertaking the sensitivity analysis prior to inversion. Further, another positive reason for choosing the optimal criterion (45) is that sufficiently large sampling times, at equidistant finite-difference time nodes, are preferably allowed to the experimentalist. Finally, it will be interesting to see in a future study how other numerical methods perform when accompanied by the optimal selection of the time measurements given by the minimisation (45). The optimal individual sampling times of the boundary pressures or fluxes presented in Table 1, were then also used in the minimisation of the least-squares functional given by equation (40) when collected measurement information of pressure and fluxes is inverted. For the constant homogeneous hydraulic conductivity, which has been retrieved in Section 7.4 from experimental measurements of the boundary pressure and 5% noisy numerically simulated measurements of the hydraulic flux, it was concluded that: (1) The measurements recorded when the pump was switched off contain more information and, therefore, a better retrieval of the unknown homogeneous coefficient than for the on pump measurements. (2) Collecting all the available data still provides a reasonable estimation of the homogeneous hydraulic conductivity of the rock without the peril that significantly, unwanted, redundant data are added to the inverse problem. (3) Any boundary measurement, individual or gathered, of the pressure and/or hydraulic flux, formulates an inverse problem whose FDM least-squares numerical solution is stable and has an accuracy comparable to the errors encountered by experimental and noisy simulated input data. For the linear heterogeneous hydraulic conductivity, which has been retrieved in Section 7.5 from exact or 1% noisy numerically simulated boundary measurements, it was concluded that: (1) The inverse formulations based on individual measurements of the boundary pressure at x = 0, hydraulic flux at x = 1 and hydraulic flux at x = 0, generate an unidentifiable problem, an unstable and stable FDM least-squares numerical solution, respectively, since the sensitivity coefficients are strongly correlated, slightly correlated and uncorrelated, respectively. (2) Gathering measurements which include hydraulic flux measurements at the active contact boundary x = 0, guarantees the identifiability of the inverse problem and improves the stability of the numerical solution. For the quadratic heterogeneous hydraulic conductivity which has been retrieved in Section 7.6 from exact or 1% noisy numerically simulated boundary measure-

758

LESNIC et al.:

HYDRAULIC

menus, the same conclusion (1) and partially the identifiability of the problem from (2), as for the linear case were maintained. However, in the quadratic case Table 6 has shown that a stable solution can be achieved only when flux measurements at the hydraulic contact boundary x = 0 are inverted and no other redundant information, which increases the ill-posedness of the problem and destroys the stability of the numerical solution, supplied by other boundary measurement, is inputted in the least-squares minimisation procedure. The paper has attempted to present a mathematical technique, which checks against experimental data for the case of constant hydraulic conductivity. By making the appropriate boundary measurements from this, it is possible to establish the internal structure of the hydraulic conductivity, when it is varying quadratically in the axial direction. Obviously, scaling up of the information is difficult and limited, but no worse than any other traditional one-dimensional experiment/analysis. Despite the restrictions that the information on the macroscale is only applicable to situations where one of the principal values of the hydraulic conductivity greatly exceeds the other values and the principal direction associated with that values remains in a fixed direction, further development of the technique has great potential. The paper is in no way claiming that the higher the order of the polynomial taken, the better it will be at inferring the hydraulic conductivity. It is just setting out to establish whether, in a sample within which the hydraulic conductivity is expressible as an unknown quadratic function, its value can be retrieved from boundary measurements alone, and what those measurements are and where they should be taken. As it is now possible to predict the one-dimensional variation of the hydraulic conductivity, at least on a quadratic level, it is envisaged that the analysis may be capable in a future study of handling piecewise constant, linear or quadratic polynomials for whose identification it might be necessary to measure internal values of the pressure inside the rock specimen. In particular, the piecewise constant interpolation approximation will be very useful in the situation of several one-dimensional samples of different, but homogeneous, hydraulic conductivities which are butted together, and so form a representative model of faults within one-dimensional fractured homogeneous rocks. In this situation, predicting the position of the faults and the constant, but different, hydraulic conductivities on either side of the faults, for which experimental data is already available, will be an important development of the present investigation. Finally, the pump-flow experiment described in Section 2 will be further improved on in the future, as follows. 8.1. A future improved transient pump-flow test system An advanced transient permeameter will be custom designed to enable the recovery of both specific storage and hydraulic conductivity coefficients as a function of position within the rock sample. The design is based

CONDUCTIVITY

OF ROCKS

around a 5.5 cm diameter core holder that accepts cores up to 10 cm long. The sample holder has three pressure ports that enable measurement of the pore pressure in time at three interior wells along the length of the sample, in addition to the standard boundary measurements. High precision and low compliance differential pressure transducers are used for the upstream and pressure port measurements, and the downstream pressure/back-pressure is monitored by a high precision gauge pressure transducer. Both upstream and downstream faces can be sealed off to minimise the upstream and downstream compliance for pulse-decay type tests. The new syringe pump is similar to the existing device, but accepts a low compliance stainless steel syringe that plumbs directly into the tubing system. A high confining pressure (3.5 MPa) and high back-pressure (1.3 MPa) can be used to ensure that the pore fluid is linearly compressible with respect to relatively small pressure changes. The back-pressure system can be replaced by a dynamic pore-pressure controller that enables downstream flux to be measured. The dynamic controller can also be used to establish constant-head or oscillating-head (sinusoidal or arbitrary waveform) boundary conditions across the sample. A common computer system controls the syringe infusion pump and the dynamic pressure controller, and acquires all the transducer and instrument test data. Acknowledgements-The authors would like to acknowledge critical, but very constructive comments of the referees.

Accepted for publication 24 December

the

1996.

REFERENCES Al-Dhahir, Z. A. and Tan, D. C., A note on one-dimensional constant-head permeability tests. Geotechn., 1968, 18, 499-505. Brace, W. F., Walsh, J. B. and Frangos, W. T., Permeability of granite under high pressure. J. Geophys. Res., 1968,73,2225-2236. Morin, R. H. and Olsen, H. W., Theoretical analysis of the transient pressure response from a constant flow rate hydraulic conductivity test. Water Resources Res., 1987, 23, 1461-1470. 4. Johnson, D. L., Koplik, J. and Dashen, R., Theory of dynamic permeability and tortuosity in fluid-saturated porous media. J. Fluid Mech., 1987, 176, 379-402. 5. Garnier, B., Delaunay, D. and Beck, J. V., Estimation of thermal properties of composite material without instrumentation inside the samples. Inr. J. Thermophysics, 1992, 13, 1097-1111. 6. Dowding, K., Beck, J., Ulbrich, A., Blackwell, B. and Hayes, J., Estimation of thermal properties and surface heat flux in carbon-carbon composite. J. Thermophys. Heat Transfer, 1995,9, 345351. 7. Olsen, H. W., Morin, R. H. and Nichols, R. W., Flow pump applications in triaxial testing. In Advanced Triaxial Testing of Soil and Rock, eds R. T. Donaghe, R. C. Chaney and M. L. Silver.

ASTM STP 977, Philadelphia, 1988, pp. 68-81. 8. Kataja, M., Rybin, A. and Timonen, I., Permeability of highly compressible porous medium. J. Appl. Phys., 1992,72,1271-1274. 9. Chen, W. H. and Seinfeld, J. H., Estimation of spatially varying parameters in partial differential equations. Int. J. Control, 1972,

15, 487-495. 10. Kitamura, S. and Nakagiri, S., Identifiability of spatially-varying and constant parameters in distributed systems of parabolic type. J. Control Optimiz., 1977, 15, 785-802.

11. Huang, C. H. and Ozisik, M. N., A direct integration approach for simultaneously estimating spatially varying thermal conductivity and heat capacity. Inl. J. Hear Fluid Flow, 1990, 11,262-268.

LESNIC et al.:

HYDRAULIC

12. Ciampi, M., Grassi, W. and Tuoni, G., The flash method and the measure of the thermal diffusivity of nonhomogeneous samples. Termotecnica, 1983, 37, 43-48. 13. Olson, R. E. and Daniel, D. E., Measurement of the hydraulic conductivity of fine grained soils. In Permeability and Groundwater Contaminant Transport, eds T. F. Zimmie and C. 0. Riggs. ASTM STP 746, Philadelphia, 1981, pp. 18-64. 14. Kranz, R. L., Saltzman, J. S. and Blacic, J. D., Hydraulic diffusivity measurements on laboratory rock samples using an oscillating pore pressure method. Inr. J. Rock Mech. Min. Sci. & Geomech. Abstr., 1990, 27, 345-352. 15. Neuzil, C. E., Cooley, C., Silliman, S. E., Ijredehoeft, J. D. and Hsieh, P. A., A transient laboratory method for determining the hydraulic properties of “tight” rocks--I. Application. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr., 1981, 18, 253-258. 16. Trimmer, D., Laboratory measurements of ultraflow permeability geological materials, Re;. Sci. Instrument, 1982, 53,-1246-1254. 17. Atkinson. C.. Hammond. P. S.. Sheooard. M. and Sobev. I. J.. Some mathematical problems fro& _the oil service industry. Proceedings of the London Mathematical Symposium, eds. R. J. Knops and A. A. Lacey, Lecture Note Series 122. Cambridge University Press, Cambridge, 1986. 18. Lin, W., Compressible fluid flow through rocks of variable

CONDUCTIVITY

19.

20.

21.

22.

23.

24. 25.

OF ROCKS

759

permeability, Lawrence Livermore Laboratory. Livermore, California, Report UCRL-52304, 1977. Hsieh, P. A., Tracy, J. V., Neuzil, C. E., Bredehoeft, J. D. and Silliman, S. E., A transient laboratory method for determining the hydraulic properties of “tight” rocks--I. Theory. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr., 1981, 18, 242-252. Smith, G. D., Numerical Solution o/Partial Dtrerential Equations: Finite Difference Methods. 2nd edn. Clarendon Press, Oxford. 1978. .Wang, H. F. and Hart, D. J., Experimental error for permeability and specific storage from pulse decay measurements. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr., 1993, 30, 1173-1176. Gill, P. E., Hammarling, S. J., Murray, W., Saunders, M. A. and Wright, M. H., User’s Guide for LSSOL, Version 1.0. Report SOL 86-1, Stanford University, 1986. Zeynaly-Andabily, E. M. and Rahman, S. S., Measurement of permeability of tight rocks. Meas. Sci. Technol., 1995, 6, 1519-3527. Beck, J. V. and Arnold, K. J., Parameter Estimation in Engineering and Science. J. Wiley, New York, 1977. Carrera, J. and Neuman, S. P., Estimation of aquifer parameters under transient and steady state conditions: 2. Uniqueness, stability and solution algorithms. Water Resources Res., 1986, 22, 21 l-227.