D Thiery
email: d.thiery@brgm.fr

SALTWATER INTRUSION MODELLING WITH AN EFFICIENT MULTIPHASE APPROACH: THEORY AND SEVERAL FIELD APPLICATIONS

Keywords: Sea water intrusion, groundwater model, multiphase flow, sharp interface, validation, field cases

Groundwater resources in Mediterranean coasts are often limited by saltwater intrusion. In order to assess accurately the available resources it is usually necessary to model the consequences of implementing new wells in the aquifer. A common approach for taking into account the presence of saltwater is to use a variable density model. This classical method however is plagued by several drawbacks: first the system must be modelled in 3D in order to account for the density effects, even when the aquifer could be considered as monolayer; second this scheme usually introduces artificial smearing due to numerical dispersion. An alternative very efficient method is to use a special multiphase approach with two separated fluids, freshwater and saltwater coupled with an equivalent retention law, which respects the saltwater and freshwater equilibrium in each model cell. This method, which may be used in steady state as well as in transient state has several decisive advantages: it may be used for monolayer aquifers, or if needed for multilayer systems, the interface between the two fluids is sharp without any smearing, and the calculations are very fast (often on the order of magnitude of 30 to 50 times faster than with the classical variable density method). The method has been verified with a set of available analytical solutions. It has been applied with success to field cases. Several cases are presented: in the island of Malta and in Camargue (South of France).

Introduction

Saltwater intrusion in aquifers has been studied for more than a century by a number of authors and many exact or analytical solutions were developed mainly before the conception of the numerical models (see [START_REF] Herzberg | Die Wasserversorgung einiger Nordsee bader[END_REF][START_REF] Glover | The pattern of fresh water flow in a coastal aquifer[END_REF][START_REF] Keulegan | Ninth progress report on model laws for density currents; an example of density current flow in permeable media[END_REF][START_REF] Henry | Transitory movements of the salt-water front in an extensive artesian aquifer[END_REF][START_REF] Henry | Interfaces between saltwater and freshwater in coastal aquifers[END_REF][START_REF] Hantush | Unsteady movement of fresh water in thick unconfined saline aquifers[END_REF][START_REF] Bear | Dynamics of fluids in Porous media[END_REF][START_REF] Bear | Hydraulics of groundwater[END_REF]Strack, 1976, among others). When numerical models became available around 1970, simplified approaches with two distinct fluids, freshwater and saltwater, were derived having in mind the common sharpness of the fluids interface. More recently, models solving the coupled variable-density groundwater flow and solute-transport equation were developed and used to study saltwater encroachment. For instance MOCDENS3D [START_REF] Oude Essink | MOC3D adapted to simulate density-dependent groundwater flow[END_REF], SEAWAT [START_REF] Guo | User's guide to SEAWAT: A computer program for simulation of threedimensionnal variable density ground-water flow[END_REF], SUTRA [START_REF] Voss | Variable density flow and solute transport simulations of regional aquifers containing a narrow freshwater saltwater transition zone[END_REF], MARTHE [START_REF] Thiéry | MARTHE: Modélisation d'aquifères par maillage rectangulaire en régime transitoire pour le calcul hydrodynamique des écoulements. Version 4.3. (MARTHE: modelling of aquifers by[END_REF][START_REF] Thiéry | Modélisation des aquifères complexes -prise en compte de la zone non saturée et de la salinité; calcul des intervalles de confiance. (Complex aquifers modelling -Unsaturated zone, salinity. Determination of confidence limits)[END_REF] FEFLOW [START_REF] Diersh | Interactive, graphics-based finite element simulation system FEFLOW for modelling groundwater flow contaminant mass and heat transport processes[END_REF], and others. This single-phase approach is appropriate when there is a large vertical transition zone due to mixing by dispersion. However it has some drawbacks:

• It is generally necessary to adopt a 3D or multilayer scheme to model the density field, even when the aquifer is monolayer with quasi horizontal flows, which implies a large number of cells and long execution time • Due to numerical dispersion, and the usually smaller vertical extension of the layers compared to the horizontal direction, the computed mixing very often tends to be too large, especially for shallow aquifers. • The coupled variable-density flow and solute-transport is prone to oscillations and not always adapted to direct steady state calculations, for instance when using Total Variation Diminishing (TVD) methods for transport.

Actual two-phase flow models were also developed for instance: INTRANS, a monolayer model [START_REF] Bonnet | L'exploitation des eaux souterraines en domaine littoral. Exemples des cotes du bassin aquitain. Présentation du modèle INTRANS[END_REF], SHARP, a quasi-three-dimensional model [START_REF] Essaid | A multilayered sharp interface model of coupled freshwater and saltwater flow in coastal systems -model development and application[END_REF], MARTHE, a full 3D model both multiphase flow or density driven model [START_REF] Thiéry | MARTHE: Modélisation d'aquifères par maillage rectangulaire en régime transitoire pour le calcul hydrodynamique des écoulements. Version 4.3. (MARTHE: modelling of aquifers by[END_REF][START_REF] Thiéry | Modélisation des aquifères complexes -prise en compte de la zone non saturée et de la salinité; calcul des intervalles de confiance. (Complex aquifers modelling -Unsaturated zone, salinity. Determination of confidence limits)[END_REF]. This approach, which applies only when the width of the freshwater-saltwater transition zone is small relative to the thickness of the aquifer, which is very often the case, has many very important attractive features:

• Possibility to adopt a monolayer scheme when physically acceptable.

• No numerical dispersion.

• Very short execution time and direct steady state calculation possible.

This paper describes the principles of the calculations and shows how a general multiphase flow model, oil-water for instance, may be transformed into a stratified flow model to simulate fresh-and saltwater flows. Then, the verification of the method with the MARTHE model is presented by comparison with analytical solutions. Finally, applications to real field cases are presented.

Description of the numerical method

The MARTHE model is a versatile, full three-dimensional finite volume model. It is possible to choose between a classical single phase coupled variable-density flow and solute-transport scheme, or a full 3D two phase approach integrating seepage faces when necessary. It is possible to combine a full 3D approach in some layers and a quasi-3D approach in others by introducing semi-pervious aquitard layers between some aquifer layers.

Governing equations

In each model cell there are two governing equations: one for freshwater and one for saltwater. Each equation, which states the conservation of mass of water (fresh or salt), is very similar to the standard equation of groundwater flow in non-saturated medium. These two equations are coupled by the hydrostatic relation which determines the position of the interface and hence the saturation of each phase.

Each fluid being of constant density, the exchanges of mass are proportional to the volumetric flows. The following equations are written for a finite difference scheme:

Freshwater equation: three dimensional conservation of freshwater mass: (subscript f)

Saltwater equation: three dimensional conservation of saltwater mass: (subscript s)

with:

T = exchange coefficient [L 2 T -1 ] H = "true" hydraulic head [L] Q = inflow of fluid [L 3 T -1 ] STO = stored volume of fluid [L 3 ] t = Time [T]
Indices: f = freshwater; s = saltwater c = centre cell; i = index of one of the 6 cells adjacent to the centre cell in 3D p = previous time step (centre cell)

T fi = (K fi .A i / dx i ) averaged between cell i and cell c T si = (K si .A i / dx i ) averaged between cell i and cell c K fi = k . (ρ w .g / µ f ) . = pressure head at the centre of the cell expressed in freshwater height [L] with:

A = exchange area between cells [L 2 ] x = distance [L] z = altitude of the model cell centre [L] dz = vertical size of the cell [L] Vol = volume of the cell [L 3 ] S S = specific storage coefficient [L -1 ] g = gravity acceleration [LT -2 ] k = permeability [L 2 ] ρ = fluid density [ML -3 ] µ = viscosity [ML -1 T -1 ] d s = saltwater relative density = ρ s / ρ f [-]
ρ w = water density in standard conditions of pressure and temperature [ML -3 ] θ f = freshwater volumetric content in a cell [-] θ s = saltwater volumetric content in a cell [-]

Equations coupling

Using a Finite Difference scheme, it is necessary to compute the pressure heads h f and h s in the centre of the cell, at an elevation dz/2 above the cell floor, dz being the height of the cell. In each cell, when freshwater and saltwater both are present, freshwater is situated above saltwater and the pressure distribution is hydrostatic within the cell. This gives the following relations:

• Confined cells: θ f + θ s = θ sat = porosity. See Figure 1. The interface elevation above the floor of the cell is z i = dz . θ s / θ sat . At this location the pressure heads h i are equal in both phases. At mid elevation the hydrostatic relation yields:

• Freshwater: h f = h i -(dz/2 -z i ) . 1 • Saltwater: h s = h i -(dz/2 -z i ) . d s h s -h f = (θ s / θ sat -0.

5) . (d s -1) . dz

Knowing h s and h f , the relation yields There are three phases: air + freshwater + saltwater, but assuming that air is at atmospheric pressure only two phases are considered. The air interface elevation is z a = dz . (θ s +θ f ) / θ sat . At this location freshwater is in contact with air and the pressure head is equal to zero which enables to compute the freshwater pressure head at mid elevation dz/2:

θ s hence θ f = θ sat -θ s . • When (h s -h f ) ≥ 0.5 . (d s -1) . dz ⇒θ s / θ sat = 1.0 ⇒ the freshwater phase disappears • When (h f -h s ) ≥ 0.5 . (d s -1) . dz ⇒θ s / θ sat = 0.0 ⇒ the saltwater phase disappears • Unconfined cells: θ f + θ s = θ liq < θ sat . See Figure 2.
h f = 0 -(dz/2 -z a ) . 1 = [(θ f + θ s ) / θ sat -0.5] . dz [1]
As for confined cells, the hydrostatic relation between freshwater and saltwater yields the same equation:

h s -h f = (θ s / θ sat -0.

5) . (d s -1) . dz [2]

Knowing h s and h f , the first relation yields (θ f + θ s ), and the second one yields θ f . Hence in both phases, the water contents θ f and θ s are known.

• When h f ≤ -0.5 . dz ⇒ (θ f + θ s ) = 0.0 ⇒ no more fluid: both freshwater and saltwater phases disappear. • When h f ≥ 0.5 . dz ⇒ (θ f + θ s ) = θ sat ⇒ the cell becomes confined In both cases, confined or unconfined, letting h_capil = h s -h f , one gets a relation in the form: h_capil = f(θ f ) which is equivalent to a retention law in a classical multiphase flow approach. There are however two important differences from the classical multiphase formulations:

• Instead of being always positive, the capillary pressure (or «suction») is positive or negative depending on the freshwater saturation θ f . This implies that when the freshwater relative saturation is less than 50 %, saltwater has a higher pressure head and plays the role of the Non Aqueous Phase (NAQ). On the other hand, when the freshwater relative saturation is greater than 50 %, it is freshwater which has a higher pressure head and plays the role of the NAQ phase. • The parameter of the "equivalent" retention law does not depend only on the porous medium and on the fluids characteristics. It also depends on the vertical dimension of the model cell.

Solution of the coupled system of equations

The equations are solved with the 3D finite volume code MARTHE described by [START_REF] Thiéry | MARTHE: Modélisation d'aquifères par maillage rectangulaire en régime transitoire pour le calcul hydrodynamique des écoulements. Version 4.3. (MARTHE: modelling of aquifers by[END_REF][START_REF] Thiéry | Modélisation des aquifères complexes -prise en compte de la zone non saturée et de la salinité; calcul des intervalles de confiance. (Complex aquifers modelling -Unsaturated zone, salinity. Determination of confidence limits)[END_REF][START_REF] Thiéry | Multiphase modelling of a gas storage in aquifer with automatic calibration and confidence limits. calibration and reliability in groundwater modelling ModelCARE 99[END_REF]. Three formulations are available for the unknowns: "freshwater head -freshwater content" (H f , θ f ) or "saltwater head -freshwater content" (H s , θ f ) or "freshwater head -saltwater head" (H f , H s ). In steady state however, only the "freshwater head -saltwater head" formulation is possible. Picard iterations are used to take care of the non-linearities. The transfer coefficient between two adjacent model cells uses an upstream weighting of the adjacent water contents (i.e. it is the water content of the upstream cell which is used). This scheme has proved to be very stable and efficient. The default formulation is the (H f , H s ) formulation, which simultaneously solves the two systems of equations which has proved most of the time to be the best.

Phase disappearance

As outlined above, according to the capillary pressure one phase (or even both phases) may disappear. For instance in the upstream region, or upper region, there is often only freshwater. On the other hand, in the downstream, or lower region, there is often only saltwater. This is taken into account classically in multiphase flow models and causes no problems. One method is to consider that there is always a very small minimum relative saturation in each phase, on the order of magnitude of 10 -5 , in order to insure hydraulic continuity and flow in the "disappeared" phase and to permit a future re-saturation. Another method is to consider that there is always a very small minimal equivalent exchange coefficient.

Verification

The presented multiphase approach has been verified on a number of analytical solutions in order to verify and outline its applicability. Two very classical saltwater intrusion problems, namely Henry test [START_REF] Henry | Transitory movements of the salt-water front in an extensive artesian aquifer[END_REF][START_REF] Henry | Interfaces between saltwater and freshwater in coastal aquifers[END_REF] without diffusion and the steady state "Glover wedge" [START_REF] Glover | The pattern of fresh water flow in a coastal aquifer[END_REF]) could be simulated without problems and are not presented here. Some verification examples are presented below in 1D, 2D or 3D, confined or unconfined aquifer, steady state or transient state.

Confined aquifer, transient state: Keulegan test

This test described by [START_REF] Keulegan | Ninth progress report on model laws for density currents; an example of density current flow in permeable media[END_REF], who gives an analytical solution, is a difficult one that has been studied by [START_REF] Shamir | Motion of the seawater interface in coastal aquifers: a numerical solution[END_REF], and by [START_REF] Mercer | Simulation of salt-water interface motion[END_REF]. It is relative to a confined aquifer of uniform thickness and infinite lateral extension. Initially there is a vertical interface, maintained by a gate, separating two compartments with freshwater in the left one and saltwater in the right one. At t = 0 the gate is removed and the interface begins to move owing to the density difference. The analytical solution neglects vertical head losses and the scheme is 1D. The test is a difficult one because there is a shock at initial time and the movement is fast. The following parameters have been selected: aquifer thickness = 10 m, hydraulic conductivity = 39.8 m/day, porosity = 30 %. The test has been modelled in transient state with adaptable time steps from 0.05 to 10 days. The 1D grid which is 140 m in length is composed of 44 cells of variable width, mainly from 2 to 4 m. Figure 3 shows that the simulated position of the interface and freshwater head are quite close to the analytical solution. It has been verified that the Ghyben-Herzberg [START_REF] Herzberg | Die Wasserversorgung einiger Nordsee bader[END_REF] approximation doesn't apply at all, which is obvious from the graph: the interface variation amplitude is 10 m which is double than the freshwater head amplitude multiplied by 40 = 1 / (ds -1).

Unconfined aquifer, steady state: Strack test

This test is described by [START_REF] Strack | A single-potential solution for regional interface problems in coastal aquifers[END_REF] who gives an analytical solution. The test relates to a shallow coastal unconfined aquifer in direct connection to the sea. There is an upstream flow of freshwater towards the sea and a single pumping well is situated in the freshwater body. Figure 4 illustrates this test: the bottom of the aquifer is 20 m below the sea level, the pumping well is situated 600 m from the coast and is screened over the whole thickness. The analytical solution neglects vertical head losses and the scheme is 2D. The following parameters have been selected for this test: hydraulic conductivity = 70 m/day, upstream flow = 1 m 3 /s per unit of width, pumping rate = 400 m 3 /s, porosity = 10%. The test has been modelled in steady state with a plane 2D grid comprising 47 rows of 63 columns of variable size. Due to the symmetry only one half of the domain was modelled with an extension of 2000 m in each direction. Figure 5 shows that the simulation of this case with 3 phases is quite accurate. The effect of the pumping well can be readily observed.

Unconfined aquifer, transient state: Hantush test

This test is described by [START_REF] Hantush | Unsteady movement of fresh water in thick unconfined saline aquifers[END_REF] who gives an approximate analytical solution. The test relates to unsteady intrusion into a thick unconfined coastal aquifer. Initially the aquifer is saturated with saltwater. At time t = 0, a recharge is applied over a strip parallel to the sea. As shown on Figure 6, a freshwater lens begins to develop, reaches the seashore after some time and then extends landwards in the opposite direction before reaching a steady state. The analytical solution neglects vertical head losses in the saturated region and the scheme is 1D. The following parameters have been selected for this test: aquifer thickness = 50 m, lateral width = 100 m, recharge strip 10 m wide situated 45 m from the coast, recharge rate 0.02 m/d, porosity = 10%, hydraulic conductivity = 100 m/d. The evolution of the interface position has been simulated during 5000 days. The test has been modelled in transient state with adaptable time steps. The 1D grid, which is 100 m in length, is composed of 102 cells of 1 m width except for two finer cells near the sea limit. Figure 7 compares the simulated interface after 10, 25, 50, 100, 200 and 5000 days with the approximate solution. This figure shows a correct simulation of this unconfined transient test.

8

Saltwater intrusion modelling with an efficient multiphase approach: theory and several field applications

3D Steady state Henry test with a seepage face and several interfaces

This test, which has no known analytical solution, has been inspired from [START_REF] Henry | Interfaces between saltwater and freshwater in coastal aquifers[END_REF] test. It was designed to illustrate the possibility of having complicated vertical geometry with a seepage face and several interfaces. The following parameters have been selected for this test: confined aquifer of uniform thickness = 100 m and lateral extension = 200 m, hydraulic conductivity = 0.01 m/s, porosity = 35%.

There is an upstream flow = 6.6 10 -3 m 3 /s per unit width at the left limit (Darcy velocity = 6.6 10 -5 m/s), a vertical seepage face at the right limit and an impermeable horizontal barrier extending from x = 105 m to 200 m at a depth between 45 m and 50 m. The modelling, which takes explicitly into account the vertical velocities, is 3D; however only a 2D cross section is modelled. The test has been modelled in steady state. The grid is composed of 20 layers of 5 m thickness and 41 columns of 5 m width, except for one cell 0.1 m wide at the sea limit. Figure 8 describes the test and displays the two computed salt wedges.

Conclusion of the verifications tests

The verification tests proved the efficiency of the described approach. It is possible to model flows in 3D, taking into account the vertical velocity components, in confined or unconfined formations, in steady state or in transient state. Seepage faces may be incorporated. Pumping wells or injection wells may be introduced in the freshwater phase or in the saltwater or in both. The calculations are very fast and it is possible to model simple aquifers as a monolayer system when the flows are quasi-horizontal.

Field case

The Malta island aquifer system This is the Malta island aquifer that was studied and modelled by [START_REF] Gutierrez | Évaluation des ressources en eau souterraine de l'île de Malte[END_REF][START_REF] Gutierrez | Évaluation des ressources en eau souterraine de l'île de Malte[END_REF]. The aquifer system was first modelled with a variable-density scheme in 3D with 15 layers and a total of 36 000 square cells of variable thickness during 24 years from 1966 to 1990. Later, a two-phase approach was used with a monolayer geometry. The aquifer system was modelled in transient state with adaptable time steps from 2 to 25 days. The 2D grid is composed of 48 rows and 50 columns of square cells of 500 m size. Figure 9 shows the island with the pumping system by deep galleries (appearing as stars on the figure) and the position of a vertical cross-section through a pumping zone. The figure displays the interface depth calculated with the two-phase approach. It shows (pale grey colour) that the interface is very high near the pumping galleries. Figure 10 displays the vertical cross-section of the heads and the time variations of the saltwater interface. Obviously the interface rose significantly between 1966 and 1990 due to the increase in the pumping rate. Both freshwater and saltwater were pumped simultaneously in some 10 Saltwater intrusion modelling with an efficient multiphase approach: theory and several field applications galleries, which was simulated as observed. The results obtained with the two-phase modelling are quite similar, with a number of cell 15 times smaller and a computation time divided by 50, to the results obtained with the variable density scheme considering the 50% salinity isoline.

The Crau aquifer (Camargue -France)

The southern border of this aquifer is the Mediterranean Sea. As can be seen in Figure 11, the hydraulic head near the sea is very flat: it is the Crau plain. A strip larger than 20 km has a hydraulic head between 0 and 1 m. As a consequence of this disposition, there is a large saltwater intrusion as shown on Figure 12. Even with this very flat configuration there were no difficulties in the calculations. The simulation was performed in steady state with a grid composed of 88 rows and 120 columns of square cells of 500 m size.

Conclusions

It has been shown that the two-phase flow method described in the paper is efficient and accurate.

Although it can be applied only when the width of the freshwater-saltwater transition zone is small, because it neglects dispersion, it displays a number of advantages. It is possible to adopt a monolayer scheme when it is physically acceptable, which dramatically reduces the number of cells. There is no numerical dispersion. The resolution scheme is stable and the computation times are short, often 30 to 50 times less than with the variable density method. It is possible to compute directly a steady state. The twophase method presented in this paper can be also used to simulate Dense Non Aqueous Phase Liquid (DNAPL) flows in aquifers.

12 Saltwater intrusion modelling with an efficient multiphase approach: theory and several field applications

Figure 1 .

 1 Figure 1. Phases relations in a confined cell.

Figure 2 .

 2 Figure 2. Phases relations in an unconfined cell.

Figure 3 .

 3 Figure 3. Keulegan test: comparison with the analytical solution. Time variations t = 1, 5, 25 and 100 days. Top = interface; Bottom = freshwater head.

Figure 4 .

 4 Figure 4. Strack test description.

Figure 5 .

 5 Figure 5. Comparison with the analytical solution of Strack (1976).

Figure 6 .

 6 Figure 6. Hantush test: top: system description. Bottom: interface position at two dates.

Figure 7 .

 7 Figure 7. Hantush test: comparison of the simulated interface after 10, 25, 50, 100, 200 and 5000 days with the approximate analytical solution.

Figure 8 .

 8 Figure 8. Theoretical test with multiple interfaces and a seepage face.

Figure 9 .

 9 Figure 9. Malta island modelling. Saltwater interface depth (m) in 1966.

Figure 10 .

 10 Figure 10. Malta island modelling. Vertical cross-section displaying: the hydraulic head (top) and the saltwater interface variations between 1966 and 1990.

Figure 11 .

 11 Figure 11. Crau (Camargue -France) aquifer: computed hydraulic head field.

Figure 12 .

 12 Figure 12. Crau (Camargue -France) aquifer: computed percent of saltwater.

SWIM. Cartagena, 2004, Spain D. THIERY

Saltwater intrusion modelling with an efficient multiphase approach: theory and several field applications

SWIM. Cartagena, 2004, Spain D. THIERY 58 Thiery 10/3/05 12:20 P gina 5

18SWIM. Cartagena, 2004, Spain D. THIERY

Thiery 10/3/05 12:20 P gina 9

18SWIM. Cartagena, 2004, Spain D. THIERY

Acknowledgements

The research described in this paper was financially supported by BRGM research project HYDROSYS.


Dominique Thiery
email: d.thiery@brgm.fr

Karine Guedeney

Multiphase modelling of gas storage in aquifer with automatic calibration and confidence limits

Multiphase flow modelling involving gas and water is widely used in gas dissolution in aquifers or in aquifer gas storage. The parameters related to the gas are usually well known but the parameters of the aquifer system are not. In order to obtain reliable forecasts, it is necessary to calibrate the multiphase model on monitored data. This can be done by automatic calibration followed by the determination of the confidence limits of the parameters, and of the confidence limits of the forecasts. This paper presents a 3D-multiphase flow model for the simulation of Gas-Water, or NAPL-Water couples. It integrates automatic calibration procedures based on Marquardt's algorithm for all types of parameters: hydraulic conductivity, which is classical, and also multiphase flow constitutive relations, aquifer compressibility and boundary limits. The model may be calibrated simultaneously on observation of different kind of data e.g. gas pressure, hydraulic heads, water saturation etc. A sensitivity analysis method, which integrates spatial or temporal dependence of the observations, follows the automatic calibration. The model is applied to the data from an actual gas storage situated in a gas cap in an aquifer near Paris (France) monitored by Gaz de France national company during more than 30 years. The automatic calibration enabled to determine the hydraulic conductivity, the rock compressibility and the lateral boundary limits. It resulted in a very accurate simulation of the gas pressure resulting from the periodical injections and withdrawals of gas. The sensitivity analysis enabled to determine the parameters that were the most important.

Notations

INTRODUCTION

Depending on the season, the consumption of natural gas undergoes significant fluctuations. The fact that sources supply tend to be remote complicates the task of adapting deliveries to seasonal and daily variations in demand. Early on, the need thus arose to constitute natural gas reserves. As far back as 1954, Gaz de France began developing techniques for storing gas underground. Underground storage facilities help balance out the sharp contrasts in annual demand between winter and summer. They also serve to meet peaks in demand on the coldest winter days. Moreover, they have to be ready to deal with exceptional situations like the record cold spells that occur statistically twice in every century. Last, they compensate for the possible failure of a supply source. The inventory reflects the evolution of the demand that has risen steadily and especially in winter.

To answer this demand, attention is to be focused both on the quality of the data, in order to have a very precise knowledge of the storage at all times, and on the quality of the reservoir modelling tool that simulates the storage behaviour and forecasts its future performance. For the forecasts to be valid, a numerical simulation should first adequately reproduce the past performance of the reservoir. The task of finding a reservoir description yielding a correct simulation of the past performance of the reservoir is called history matching (HM). For the HM to be correct, we need to have accurate measured pressure, appropriate reservoir modelling tool and suitable HM technique by automatic calibration. The difference between simulated outputs and measured data, i.e. the data residuals, will then be acceptably small, and the performance forecasts reliable.

MULTIPHASE MODELLING OF FLOWS

Equations to be solved

At each time step it is necessary to solve a system of non linear equations resulting from (1) the volume conservation of water in each cell, (2) the mass conservation of gas in each cell and (3) the equations of state for gas and water. These equations, expressed in finite volume, are :

Volume conservation of water (subscript w) :

Mass conservation of gas (subscript g) :

With : p = Index for the beginning of the time step i = Index of an adjacent cell ρ gi = Weighted gas density with adjacent cell i d gi = Av. relat. density with cell i

T wi = k ρ w / µ w K rw A / dx T gi = k (ρ w / µ g) K rg A / dx STO w = Vol [(θ w -θ wp )] + θ w (β w + β por )(H w -H wp )] ( ) ( ) [ ] ( ) dt STOg d Q z z d h h T g i i gi g i g gi i g = + ∑ - + - ρ ( ) dt w STO d Q w i H w i H wi T = + ∑       - STO g = Vol [(θ g ρ g -θ gp ρ gp ) + ρ g θ g β por (h g -h gp )]

Resolution of the system of equations

The equations are solved with the 3D finite volume code MARTHE described by [START_REF] Thiéry | Modélisation d'Aquifères par maillage Rectangulaire en régime Transitoire pour le calcul Hydrodynamique des Ecoulements. Version 4.3[END_REF]Thiéry ( , 1993b)). Three formulations are available for the unknowns : (water head -volumetric content) or (Gas pressure -volumetric content) or (water head -gas pressure). Picard iterations are used to take care of the non-linearities. The transfer coefficient between two adjacent cells uses an upstream weighting of the adjacent relative conductivities (i.e. it is the relative conductivity of the upstream cell, which is used). This scheme has proved to be very stable and efficient. For the present case a "water head -gas content" scheme has been chosen. In each cell the unknowns are the head of the water phase Hw and the gas content θg.

The volume and mass conservation equations are written with these variables converting θw to θg using the relation ω = θw + θg and converting hg to Hw using the suction h c definition.

Then θg is eliminated by adding the 2 nd equation to the 1 st one multiplied by ρg. One obtains a system of n equations, n being the number of cells, with the only variable Hw which is solved classically. From the computed values of Hw the variables θw are immediately obtained explicitly in each cell from equation 1. Then in each cell hg is obtained immediately from suction definition and θg from ω = θw + θg.

Geometry of the model, boundary conditions and initial parameters

The underground gas storage of Saint-Illiers-la-Ville is located 50 km west of Paris (France) in the Upper Oxfordian. It is a fairly radial sandstone reservoir 30 m thick situated 340 m deep below sea level. The aquifer gas storage now presents a central gas cap surrounded by an active aquifer. The maximum inventory is about 1500 millions of normal cubic meters (MMm3 (n)) and the working gas volume is roughly 650 MMm3 (n). Figure 1 shows a top view of the reservoir. The reservoir has been modelled in radial coordinates. The grid is as follow : a circular cell of radius 500 m, where gas injection or withdrawal take place, surrounded by 22 ring cells 50 m wide. Then larger ring cells extend the domain to a radius of 82 km. (See [START_REF] Ory | A Semi-Automatic History Matching Technique Applied to Aquifer Gas Storages[END_REF]. The vertical thickness of each ring cell is 30 m and the top of the formation is decreasing from -340 m (below msl) in the centre to -460 m at a distance of 1.7 km. In each cell the initial value of permeability k has been fixed to 2 Darcys (1.97x10 -12 m 2 ) and the initial value of the porosity ω to 28 %. The compressibility of the pores β por has been fixed to 2x10 -5 m -1 and the compressibility of the water β w to 0.5x10 -5 m -1 . The constitutive laws for the relative permeability, the capillary pressure (suction), the correction of the gas compressibility with the pressure (the natural gas is not a "perfect gas") and the gas viscosity have been provided by Gaz de France. The gas pressure is monitored at the head of the observation well SI9 that corresponds approximately to the 2 nd cell of the model (525 m from the injection). The pressure history is corrected to a reference datum of 360 m below sea.

The gas rate of injection or withdrawal and the gas pressure measurements are available during a long period of 33 years extending from 1965 (when natural gas injections started) to 1997. 1656 pressure measurements are available which corresponds to weekly observations. At the beginning of the storage (1965) there is no gas and the water phase is immobile with a hydraulic head estimated to 141.11 m which corresponds to a pressure of 616.11 m (60.44 bars) at the outer limit where the mean altitude is -475 m. At the outer limit the hydraulic head is prescribed and the gas content is prescribed at 0. Considering the uncertainty of the hydraulic head at the outer limit, and the relative altitude where the gas pressure is monitored, a correction has been applied to this head with an initial value H0 = +30 m (+2.94 bars). The simulation with MARTHE code has been performed with automatic time steps prescribed to range in the interval [0.5 -4 days]. On a classical PC computer (Pentium II) approximately 6 min. of CPU time are needed for a simulation of the total period.

AUTOMATIC CALIBRATION OF THE PARAMETERS

Marquardt algorithm [START_REF] Marquardt | An algorithm for Least Squares Estimation of non linear parameters[END_REF]Press et al., 1992), which has been implemented in MARTHE code, has been used to calibrate four parameters : k, ω, β por and H0. This algorithm, which is more efficient than the former method described by Thiéry (1993a[START_REF] Thiéry | Calibration of groundwater models by optimisation of parameters in homogeneous geological zones[END_REF], necessitates the computation of the gradients of the observations (the monitored pressure) with respect to each parameter in order to get a Hessian. The gradients have been computed by finite difference, which necessitates, for each iteration of the algorithm, one simulation per parameter plus a reference simulation. [START_REF] Anterion | Use of Parameters Gradients for Reservoir History Matching[END_REF] present a faster method for the direct calculation of all the gradients in a single simulation which could be used.

Table 1 displays the evolution of the calibration of the four parameters. The standard deviation of the observations σ obs is 72.05 m (7.07 bars) and the initial Root Mean Square of Error RMSE is 2.53 bars. The calibration criterion is defined as (RMSE / σ obs ) 2 which is the percentage of residual errors (unexplained variance). After two iterations the criterion is reduced from 12.8% to 0.96%. After four iterations the parameters are calibrated and the RMSE is equal to 0.6 bar. Figure 2 displays the comparison of the monitored gas pressures to the simulated ones. As the simulation with MARTHE code is very accurate the two curves can hardly separated on the graph. It should be noted that the low pressures due to a strong depletion in the winter of 1984-1985 and the peak resulting from a significant increases in maximum inventory in the summer of 1995 are well simulated. Figure 3 displays the initial and final series of pressure deviations.

SENSITIVITY ANALYSIS

The sensitivity analysis is based on a local linearisation of the model in the vicinity of the optimal set of parameters. The method, which uses also the Hessian matrix, has been described among others by Thiéry (1993b, c), [START_REF] Leijnse | Evaluation of pumping tests : identification of parameter values and their reliability. Improvements of methods of Long Term Prediction of Variations in Groundwater Resources and Regimes due to Human Activity Proceeding of the Exeter Symposium[END_REF]. Classically the method is applied assuming that the observations, or the simulations residuals, are independent which is not acceptable most of the time. Instead Thiéry (1993b, c) shows that it is necessary to integrate their spatial or temporal dependence. This could be done using the covariance matrix of observations. However this not feasible considering the large number of observations which would yield a huge matrix. An approximate method has been used instead. The analysis of the monitored pressures shows that two pressures separated by 10 time lags (approx. 70 days) have an autocorrelation coefficient not less than 0.5. It can then be roughly considered that a set of 10 consecutive observations forms only one independent observation. The total number of independent observations reduces from 1656 to 166, which is realistic because it corresponds to 5.2 independent observations per year during 31.8 years. Using this method the standard deviation of the parameters has been calculated. Three parameters, k, ω and β por had a logarithmic transformation. For these parameters relative standard deviation is displayed. The relative standard deviation is based on the standard deviation of the logarithm. It defines the range which has a probability of 68.2% to contain the true value. 19671969197119731975197719791981198319851987[START_REF] Anterion | Use of Parameters Gradients for Reservoir History Matching[END_REF]1991199319951997 Deviation (Bars) strongly correlated k and β por (correlation -0.923). The negative sign indicates that a large compressibility has a similar effect as a large permeability : both parameters tend to dampen pressure variations. The sensitivity analysis shows that the rock compressibility is not very well identified. This parameter could be fixed to a reference value. The other parameters would then be still better identified.

CONCLUSIONS

The multiphase flow numerical method used in MARTHE code performs well and the automatic calibration scheme is efficient even in highly non linear systems. It enables to simulate very accurately the variations of gas pressure in a storage during a long monitored period which proves its ability to yield reliable predictions. The presented sensitivity analysis is useful in determining the most significant parameters.

  hydraulic pressure head [L] H w = Water hydraulic head = hw + z [L] STO g = Stored mass of gas [M] STO w = Stored volume of water [Gas density [ML -3 ] ρ w = Water density [ML -3 ] ω = Porosity = θ w + θ g [ L 3 L -3 ]d g = Gas relative density = ρ g / ρ w [-] Q w = Volumetric flow of water [L 3 T -1 ]Q g = Mass flow of gas [MT -1 ] g = Gravity acceleration [LT -2 ]

Fig. 1

 1 Fig. 1 Top view of Saint-Illiers reservoir: isobathes and observation well. The production wells (not shown here) are located at the top of the reservoir.

Fig. 2

 2 Fig. 2 Monitored and simulated gas pressure at the observation well.

Fig. 3

 3 Fig. 3 Initial (large amplitude) and final (small amplitude) series of pressure deviations.

Table 1

 1 Automatic calibration of parameters.

		I t e r a t i o n :	01234
	Permeability	(D)	2	0.519	0.669	0.688	0.729
	Porosity	(%)	28	31.3	29.6	24.5	22.45
	Head correction	(m)	30	21.2	10.4	-2.78	-7.89
	Pore Compressibility (m -1 )	2	7.82	2.68	1.36	0.857
	Criterion	(%)	12.8	3.60	0.963	0.713	0.693
	RMSE	(m)	25.83	13.66	7.07	6.09	5.99
	RMSE	(bars)	2.53	1.34	0.69	0.597	0.588

Table 2

 2 Sensitivity analysis of the parameters.

				Head	Compressi-	Optimal	Standard	Relative
		Permeability Porosity	Correction	bility	Value	Deviation	Standard
								Deviation
	Permeability	1				0.729		4.89 %
	Porosity	0.064	1			22.45		2.74 %
	Head correction	0.619	0.791	1		-7.89	2.55 m
	Compressibility	-0.923*	0.054	-0.462	1	0.857		34.5 %

Acknowledgements

The research described in this paper was financially supported by the Research Division of BRGM (publication n° 99021). The research has been realised in consultation with Gaz de France Research and Development Division.


Wolfram Kloppmann

Dominique Thiéry

Christophe Kervévan

Anne Bourguignon

Philippe Négrel

Joël Casanova

Chemistry-transport coupled modelling of the Äspö groundwater

come

INTRODUCTION

The long term prediction of the geochemical evolution of groundwater bodies within potential host formations is of particular concern for site assessment in the field of nuclear waste disposal (Boulton et al. in press). Evaluating the impact of changes in climatic/hydrological surface conditions on the groundwater quality has to make use of a combined approach of flow and transport modelling taking into account chemical reactions. An adapted methodology of chemistry-transport coupled modelling based on Specific Chemical Simulators (SCS) has been developed in the last few years (Kervevan et al. 1998). It was applied to the site of Äspö (Sweden) for which a detailed site investigation yielded huge quantities of geological, hydrogeological and geochemical data from the crystalline basement up to depths of >1000 m (Laaksoharju et al., 1999a). The paleohydrological/hydrochemical evolution since the last glacial maximum is particularly complex due to the presence of deep seated brines, meltwater injection, several phases of immersion of Äspö island by the subsequent stages of the Baltic Sea and a final emergence of the site allowing mixing with meteoric waters at the surface (Laaksoharju et al. 1999a).

In the following, we describe, as an example for the application of the SCS approach, the attempt to reconstitute the geochemical and flow conditions at Äspö from the last deglaciation to the modern natural conditions prevailing before the construction of the underground laboratory.

FLOW AND TRANSPORT MODELLING

The modelled zone is a square of 8 km length and extends to a depth of 1600 m below the present day sea level. All geometric (elevation) and hydraulic data (3 components of conductivity, hydraulic head, salinity) were provided by U. Svensson (1997). A compromise had to be found between computing time, spatial and temporal resolution of the transient state flow modelling and the complexity of the geochemical model. The final configuration adopted is a horizontal irregular grid of 26 x 26 cells and a vertical grid of 11 layers from 3 to 300 m of thickness.

The hydrodynamic code used is the finite difference model MARTHE (Thiéry, 1993a, b) taking into account an unsaturated zone, density effects, direc Chemistry-transport coupled modelling of the Äspö groundwater system (Sweden) since the last glaciation W. Kloppmann, D. Thiéry, C. Kervévan, A. Bourguignon, P. Négrel, J. Casanova BRGM, BP 6009, F-45060 Orléans Cedex 2, France Water-Rock Interaction 2001, Cidu (ed.), A.A. Balkema (Swets & Zeitlinger), Lisse, 181-184. ABSTRACT: Based on the specific coupled modelling approach developed at BRGM, a reaction-transport model of the Äspö groundwater system (Sweden) has been developed; it simulates a continuous scenario from 12.000 BP to the natural modern situation. Examples of the results (evolution of Na, 18 O, calcite) are presented. Na, Ca, Cl and stable isotopes are mainly related to mixing between the different mixing endmembers subsequently entering the system (glacial meltwater, Baltic Sea stages, meteoric water), calcite contents, pH and silica are determined by water rock interaction.

tional hydraulic conductivity (in fractures) and multi-component transfer integrating advection and mechanical dispersion.

The initial concentrations in chemical species have been considered as constant in the whole domain depending only on whether a cell is considered as crossed by fractures or as composed of matrix only. The 18 O signature is -8.9‰ in the initial dense brine; -20.45‰ in the initial mixed brine and -32.0 ‰ in the initial glacial meltwater. The sea level, the concentrations and the stable isotope composition of the Baltic varied considerably with time. The boundary conditions were accordingly adjusted every 1000 year period. Salinity and hydraulic head calculations were performed 320 times for 1000 years with a frequency decreasing from the beginning to the end of each 1000 year stress period. The hydrochemical calculations were performed every 1000 years. A test with twice as much time steps yielded no significant difference for the concentrations of aquatic species.

GEOCHEMICAL MODELLING

Development of the "PAGEPA" Specific Chemical Simulator (SCS)

The general approach chosen by BRGM is to design a dedicated geochemical code (so called "SCS") developed specifically for each particular application that only takes into account the relevant processes to be considered [START_REF] Kervévan | SCS: Specific Chemical Simulators dedicated to chemistry-transport coupled modelling Part III -Coupling of SCS with the hydrotransport modelling software MARTHE[END_REF]. One of the advantages of this approach is a significant improvement in computer efficiency that in turn allows relatively complex 3D coupled modelling.

This approach thus requires us to develop a new SCS for each new application. The combined use of ALLAN 3.1 and the NEPTUNIX 4 software package automatically generates the coded part of the simulators, as long as the processes involved in the physical system being modelled can be described by a set of algebraic and/or ordinary differential equations.

NEPTUNIX 4 is a general solver which automatically generates a FORTRAN code solving the selected set of equations. ALLAN 3.1 functionalities include graphical interface, description language, translator for NEPTUNIX 4, graphical assembling of models (including the generation of coupling equations), exploitation of the simulators and postprocessing of the generated data.

The SCS developed for the Äspö site includes 24 aqueous reactions of 12 elements (H, C, O, Ca, Al, Si, Mg, Na, Cl, K, Fe, S ox ), the 18 O isotope, 7 precipitation/dissolution reactions (calcite, kaolinite, clinochlore, goethite, quartz, K-feldspar, anorthite) and two cation exchange reactions (K and Na versus Ca).

Chemical boundary conditions

All calculations are based on the reference waters defined by Laaksoharju et al. (1999b) as the groundwater samples closest to the hypothetic pure mixing endmembers. In order to characterise the geochemical system in terms of major species to be taken into account, we performed preliminary EQ 3/6 (Wollery, 1992) calculations of chemical speciation and saturation states of all the reference waters. The species whose relative concentration was less than 10 -2 with respect to the major species were not considered in the SCS.

The different inputs of marine and continental waters into the system during the past 12,000 years were defined and taken into account in the coupled simulations. In order to finalize our selection of aqueous species actually taken into account in the SCS, EQ3/6 mixing calculations were also performed to simulate chemical compositions, speciation and saturation states of the different Baltic Sea stages.

The three predominant lithologies encountered in the four main boreholes and in the tunnel are the Äspö diorite, the Småland granite and the Finegrained granite (Mazurek et al., 1995 ;Stanfors et al. 1999). Extensive studies characterise the mineralogy of the fracture fillings (Mazurek et al., 1995;Tullborg, 1989;Wikman et al., 1988 ;[START_REF] Landström | Results from a geochemical study of zone NE-1, based on samples from the Äspö tunnel and drillcore KAS 16 (395 to 451M)[END_REF]. A selection of the occurring minerals, taking into account the constraints related to computing, was integrated into the SCS.

Chemical evolution

As a general observation we can state that certain parameters, the conservative species like Cl, but also Na and bivalent cations, reflect essentially the mixing processes in the course of the subsequent stages of glacial/meteoric/marine waters injection whereas others (Si, pH, mineral precipitation) are closely related to the geometry of the reservoir and in particular of major fractures.

Figure 2 shows, as example of the modelling results, the distribution of the Na concentrations of the δ 18 O and of calcite for 11,000, 6000 and 1000 years BP. All block diagrams show E-W sections cutting through Äspö Island. Na: The patterns of the Na distribution are, much like those of Ca, closely related to the chlorinity. The initial scenario at 12,000 BP represents the stratification of glacial meltwater over deep brines. The fluctuations between 12,000 and 10,000 BP can be attributed to density effects. The influence of Yoldia sea becomes visible at the 9000 BP stage by increasing salinity of the surface layers. Ancylus lake waters (8000 BP) dilute the superficial strata and the Litorina stage increases again the the concentration in the surface waters. From 5000 BP onward, the progressive invasion of the system by low salinity meteoric waters from the West is driven by the increasing hydraulic gradient due to lowering of the sea level. The Baltic 2 stage introduces a second low Na plume from 1000 BP to present. Low chlorinities of the modern Baltic Sea (since 1000 BP) and of meteoric water determine the present day conditions in the surface layers. We observe a residual local salinity high near the Äspö island where Baltic 1 (or Litorina sea) water is conserved between the intrusion of meteoric water and modern Baltic Sea water. Calcite: According to the defined initial conditions, calcite is limited to major fractures at 12,000 BP. Calcite precipitation in the fissure porosity of the matrix starts from the surface during the Baltic ice lake stage (Fig. 2). During the whole evolution of the site, we observe a rather continuous progress of a horizontal "calcite precipitation front" moving from the surface to depth. Obviously, calcite precipitation is fairly independent from the subsequent marine/meteoric entries into the system and from the major ion concentration distribution. From 3000 BP to present, calcite contents in the surface layers reach values comparable with the initial concentrations defined for the major fractures.

Silica: Si concentrations do not vary significantly during the evolution of the site. This is only partly due to the fact that Si variations in the initial fluids entering the system are small. The crucial factor seems to be a buffering of the concentrations by mineral dissolution/precipitation (according to the chemical model, silica concentrations are controlled by quartz).

pH: The total range of pH variations is relatively restrained (around 1 pH unit). The pH is definitely not an indicator of mixing but it's distribution is closely related to mineral equilibria. The pH distribution is, much like silica concentrations and calcite contents, determined by the geometry of major fractures. In the fractures we observe the highest values around 9, which is in the uppermost part of the range of measured field values.

δ 18 O: Initially, the system is invaded by isotopically extremely depleted glacial meltwater, which mixes at depth with the more enriched brines (Fig. 2). The pattern of isotopic contents till 10,000 BP remains rather constant and is very similar to the sodium/chloride distribution (Fig 2). Both tracers can in fact be regarded as conservative under the conditions prevailing in the system. The impact of the 18 O-enriched Yoldia sea changes profoundly the isotopic contents of the system and the influence of the Yoldia water reaches depths of >900 m. At 8000 BP, the isotopic composition of the system is ho-mogenised up to the bottom of the block where slightly enriched brine remains. The impact of isotopically enriched Litorina sea (δ 18 O = -5.5‰) becomes predominant from 7000 to 5000 BP and the isotopic composition of the Baltic Sea is constant up to present day conditions around -6‰. Meteoric water which is depleted (δ 18 O = -10.5‰) with respect to Litorina seawater progressively invades the system from the continent. Present day conditions are somewhat similar to chloride distribution but contrasts of the isotopic composition are much lower than for chlorinity.

CONCLUSIONS

The coupled hydrodynamic-SCS modelling approach has been successfully applied to a complex setting of a fractured coastal aquifer influenced by glacial melt water, meteoric water and successive stages of brackish to saline seawater intrusions over 12,000 y. The number of chemical reactions taken into account (in both space and time) are at present mainly limited by the constraints of computing capacity. Some simplifications of the natural system that were adopted in the present study, as neglecting redox processes though they seem to play a particular important role in this site [START_REF] Banwart | The redox experiment in block scale. Final reporting of results from the three year project[END_REF], could be overcome in a near future considering the fast evolution of computer capacities.

-11,000 yr -6,000 yr -1,000 yr (25-88-11). Stockholm Wollery T.J.1992. EQ3/6, A Software Package for Geochemical of Aqueous Systems : Package overview and installation Guide (Version 7.0). Report UCRL-MA-110662-PT-I. Livermore, California: Lawrence Livermore National Laboratory.

Figure 1 .

 1 Figure 1. Situation of the study site

  Svensson, U. 1997. A regional analysis of groundwater flow and salinity distribution in the Aspö area. SKB Technical Report 97-09. Stockholm. Thiéry D. 1993a. Modélisation des aquifères complexes -Prise en compte de la zone non saturée et de la salinité. Calcul des intervalles de confiance. Hydrogéologie 4: 325-336. Thiéry, D. 1993b): Tridimensional and multi-layer modelling of transfers in unsaturated porous medium. GEOCONFINE Symposium international Géologie et confinement des dé-chets toxiques. Montpellier, juin 1993: 467-472. Rotterdam: Balkema. Tullborg, E-L. 1989. Fracture fillings in the drillcores KAS05-KAS08 from Äspö, Southeastern Sweden. SKB HRL Prog. Rep. (25-89-16). Stockholm. Wikman, H., Kornfält K.A., Sgu L., Riad L., Munier R. & E.V. Tullborg 1988. Detailed investigation of the drillcores KAS 02, KAS 03 and KAS 04 on Äspö island and KLX 01 at Laxemar. SKB HRL Prog. Rep.

	Calcite
	2.0E-03
	1.6E-03
	1.2E-03
	8.1E-04
	4.1E-04
	1.0E-05
	mol/m 3 of
	soil
	Figure 2. E-W section through the study area (lateral extension: 8 km x 8km, heith of the block: 1800 m): Distribution of solid cal-cite, dissolved Na and δ 18 O at 11,000, 6000and 1000 BP.
	Mazurek M., Bossart P.& T. Eliasson 1995. Classification and
	characterization of waterconducting features at Äspö: re-
	sults of phase I investigations. SKB HRL Prog. Rep. (25-95-
	03). Stockholm.
	Stanfors R., Rhén I., Tullborg E.-L. & P. Wikberg 1999. Over-
	viev of the geological and hydrogeological conditions of
	the Äspö hard rock laboratory site. Appl. Geochem. 14:
	819-834.

Acknowledgements: The work was conducted under the PAGEPA project (F14W-CT96-0038), a European research project financed through the EU 4 th Framework Program.


Dominique Thiéry

Catherine Golaz

Alexis Gutierrez

Wieslaw Fialkiewicz

Coralie Darsy

Christophe Mouvet

Igor G Dubus
email: i.dubus@brgm.fr

Refinements to the MARTHE model to enable the simulation of the fate of agricultural

Keywords: Model, verification, soil, unsaturated zone, saturated zone, pesticide

Introduction

Most pesticide leaching models used to assess the risk of compounds being transferred to groundwater resources within the context of pesticide registration or research on pesticide fate include a detailed description of processes affecting pesticide fate in the soil, but their remit is limited to 1D simulation in the first few metres of soil. Transport and reactive processes occurring in the deeper unsaturated and saturated zones are ignored. Conversely, deterministic groundwater models (e.g. MODFLOW, SUTRA) usually do not allow the simulation of the fate of pesticides at the soilatmosphere boundary and in soil, where attenuation processes key to the overall fate of pesticides occur. The present paper reports on the refinement of the 3D groundwater model MARTHE to allow the simulation of pesticide fate from the soil surface to and in groundwater. A rigorous verification exercise to ensure that the new code was correctly implemented was undertaken.

The MARTHE model

MARTHE (Modelling Aquifers with an irregular Rectangular grid, Transport, Hydrodynamics and Exchanges) is a versatile 3D groundwater model designed to enable the hydrodynamic and hydrodispersive modelling of groundwater flow in hydrosystems [START_REF] Thiéry | Logiciel MARTHE -Modélisation d'Aquifère avec maillage Rectangulaire, Transport et HydrodynamiquE, version 4.3[END_REF][START_REF] Thiéry | Modélisation des aquifères complexes. Prise en compte de la zone non saturée et de la salinité. Calcul des intervalles de confiance[END_REF][START_REF] Thiéry | Modélisation 3D des écoulements en zone non saturée avec le logiciel MARTHE Version 5.4[END_REF][START_REF] Thiéry | Modélisation 3D du transport de masse avec le logiciel MARTHE Version 5.4[END_REF][START_REF] Thiéry | Consideration of vegetation effects in version 6.2 of the MARTHE model[END_REF]. Hydrodynamic calculations are based on a fully implicit scheme with a 3D finite volumes discretization. The system simulated may be represented with a classical confined/unconfined scheme or by a unsaturated zone -saturated zone continuum. In this latter case Richards' equation is solved with retention and relative hydraulic conductivity relations selected from a number of formulations (Van Genuchten, Brusaert, Brooks and Corey, Gardner, homographic, power law etc.). Advective, diffusive and dispersive transport can be simulated using three different techniques (donor cell, total variation diminishing, method of characteristics using particles) depending on the situation considered (predominance of convection or dispersion). Calculations for energy, temperature, mass and water fluxes are simultaneously fully coupled within the model. Interactions with the atmosphere (rainfall, potential evapotranspiration) may be integrated through coupling to hydroclimatic balance models which compute actual evapotranspiration, infiltration and runoff. The aquifer system may be coupled to a river or a drain network which receives runoff and exchanges flow along the reaches. The horizontal discretisation uses an irregular rectangular grid which may integrate nested grids in zones of particular interest or near singularities. The vertical discretisation may be represented as a fully 3D system or a complex multilayer system with possibilities of disappearing layers and short circuits. MARTHE integrates pre-and post-processors to facilitate the preparation, management and graphical processing of input data and modelling results. The system also allows imports/exports to be made with the GIS software MapInfo and integrates its own tool for interpolation of scattered data. MARTHE has been applied to over 200 aquifer systems since its first release in the late 1980's. Fields of application have included i) the management of groundwater resources (e.g. water balance assessments, evaluation of the potential impacts of abstractions); ii) civil engineering and mining work (e.g. impact of underground activities); and, iii) environmental assessments (both point and diffuse pollution). The scale of applications of the model has ranged from detailed investigations into water movement in soil (soil volumes of a few cm with cells of a few mm) to simulation of large multilayer aquifer systems (150,000 km 2 , 1,500-m thickness). The temporal scale of application ranged from seconds to 15,000 years for long-term geochemical predictions.

Addition of new functionalities enabling the simulation of the transfer of water and pesticides in the root zone -unsaturated zone -saturated zone continuum

MARTHE was initially designed as a 3D groundwater model for water transport in groundwater. Although the remits of the model have largely broaden since the early versions of the model (e.g. through simulation of fluxes in the unsaturated zone and of solute transport or the possibility to couple MARTHE to stand-alone (geo)chemical models), the model was found in 2000 to be unsuitable for adequately simulating the fate of pesticides from the soil surface to and in groundwater. Also, MARTHE was implementing a simplified description of the water balance at the upper boundary and descriptions of physical and biological processes known to drive the overall fate of pesticides (sorption, degradation) were not included in the model. MARTHE was therefore adapted as part of the European project PEGASE [START_REF] Mouvet | PEGASE: Pesticides in European Groundwaters: detailed study of representative aquifers and simulation of possible evolution scenarios[END_REF] through the addition of new code for cropping and pesticide sorption and degradation.

New functionalities for cropping

The presence of vegetation significantly influences water and pesticide fluxes in soil through i) plant transpiration; and, ii) plant uptake of water and solutes. The following features were implemented in MARTHE to improve water and pesticide fluxes close to the soil surface: i) breakdown of evapotranspiration into soil evaporation and transpiration by vegetation; ii) water uptake by vegetation over a variable depth that depends on the condition of root development and vertical root distribution; iii) restriction of water uptake according to the state of hydric stress experienced by vegetation; and, iv) pesticide uptake through the root system, depending on the product considered.

New functionalities for pesticide sorption

The schemes selected for description of sorption for pesticides were those of Langmuir and Freundlich. The Langmuir sorption isotherm expressing the relation between the solid phase concentration S M and liquid phase concentration C at equilibrium is given as:

C K 1 C K Q S L L M ⋅ + ⋅ ⋅ =
where: S M = Solid phase mass concentration [kg/kg], C = Liquid phase volumetric concentration [kg/ m 3 ], Q = Total sites for solid phase concentration [kg/kg] and K L = Langmuir constant [m 3 /kg] The Freundlich sorption isotherm writes: S M = K F . C B where: B = Freundlich exponent (generally < 1), K F = Freundlich constant in unit [m 3 /kg] B .

New functionalities for pesticide degradation

Two types of degradation schemes were integrated in MARTHE. The first scheme is sequential (or "linear"), in that compound #1 is assumed to be degraded into compound #2, which in turn is degraded into compound #3, etc. The second scheme enables the simulation of the degradation of a compound into multiple metabolites. The current version of the model does not allow the combination of the two degradation schemes although modifications to the code will enable the simulation of complex degradation pathways using complex kinetics in the near future.

New functionalities for the effect of water content and temperature on degradation

Pesticide degradation is influenced by temperature and soil humidity conditions. A review of the different equations used to describe the effect of temperature and soil moisture variations on degradation was undertaken and supported the selection of five correction schemes for temperature (LEACHP, AGRIFLUX, MACRO, WAVE and PELMO) and six schemes for humidity (LEACHP, AGRIFLUX aerobic, AGRIFLUX anaerobic, MACRO, WAVE and PELMO).

Verification exercise

Cropping New subroutines for crop development were evaluated using data from the La Côte St André (Isère, France) experimental site [START_REF] Normand | Nitrogen-15 tracers combined with tensio-neutronic method to estimate the nitrogen balance of irrigated maize[END_REF] where water fluxes on a bare soil plot and a plot cropped with maize were monitored over three successive cropping seasons. All data supporting the evaluation were provided by B. Normand (LTHE, Grenoble, France). Figure 1 presents measured and simulated water contents for the two plots. ) in a horizontal column which received an input of NH 4 + at its left boundary (Figure 3). , 2000) were carried out to demonstrate the effects of humidity and temperature corrections on model predictions (Figure 4).

Conclusions

The verification exercise undertaken using analytical solutions published in the literature (sorption, degradation) and experimental data (cropping, degradation) demonstrated that the new code added to MARTHE has been correctly implemented in the model. In addition, the close fits obtained between predicted and theoretical/measured data confirmed the sound numerics of the model. These refinements and their testing mean that that model can now be used with confidence to simulate the fate of pesticides from the soil surface to and in the groundwater. The model is expected to be a useful tool in the modelling of pesticide fate through the whole continuum root zone -unsaturated zonesaturated zone, in 1, 2 or 3 dimensions. In particular, the model may be used to i) assess the risk of pesticide leaching to groundwater resources at a range of scales; ii) provide explanations for detection / non-detection of pesticides in monitoring programmes; and, iii) manage inputs of pesticide to a particular area within a context of sustainability of groundwater resources.

Figure 1 .Figure 2 .

 12 Figure 1. Observed (symbols) and simulated (line) water contents at La Côte Saint André in 1992 (left: bare soil; right: maize) Sorption The correct implementation of the Langmuir and Freundlich equations were tested against solutions provided by Zheng C. et al. (1998) and Huang K. et al. (1998), respectively (Figure 2).Langmuir sorption

Figure 3 .

 3 Figure 3. Comparison of MARTHE predictions (lines) with analytical calculations (symbols) obtained using the CHAIN code from Van Genuchten. (squares: NH 4 + ; diamonds: NO 2 -; triangles: NO 3 -) Effect of temperature and humidity variations on pesticide degradation A number of simulations based on the description of the leaching of bentazone in a 2-m deep Dutch soil (Vredepeel dataset; Boesten J.J.T.I., 2000) were carried out to demonstrate the effects of humidity and temperature corrections on model predictions (Figure4).

Figure 4 .

 4 Figure 4. Influence on MARTHE predictions for concentrations of bentazone in the Vredepeel soil profile of the effects of temperature and saturation corrections of degradation

Acknowledgements

This development and evaluation work was undertaken within the scope of the European project PEGASE (PEGASE: Pesticides in European Groundwaters: detailed study of representative Aquifers and Simulation of possible Evolution scenarios; http://www.brgm.fr/pegase). The funding of the research by the European Commission under the 5 th Framework Program (EU contract number: EVK1-CT1999-00028) is gratefully acknowledged.


D Thiery

SEMI-GLOBAL HYDROLOGICAL MODELLING BY COUPLED SIMULATION OF THE VARIATIONS OF WATER STORAGE IN THE UNSATURATED ZONE AND OF THE PIEZOMETRIC LEVEL EVOLUTION

Keywords: Hydrology models, unsaturated zone, water resources data processing parameter estimation

A new type of hydrological model -called "semi-global" has been created to attempt at a better understanding of the physical behaviour of a drainage basin. This model. called "Gardensol" simulates the basin globally. in its spatial extension, but it reproduces. at least approximately, the real phenomena with physical parameters capable of being measured in the field :

-vertical permeability at saturation. -residual water content, corresponding to an almost nil vertical infiltration permeability. As it can be calibrated simultaneously on the water content variations in the unsaturated zone and on the piezometric variations, the model is much more reliable and more robust than a rainfall-aquifer level model where the aquifer levels are reproduced only from a usually badly known storage coefficient. The estimation of the recharge, in particular, must be much more exact and the model can be used by extending data to analyse the resource fluctuations according to the climatic variabili ty. Three application examples are given : coupled simulation of the water storage in the unsaturated zone and of the piezometric level at Kezieres-Ies-Clery (Loiret department, France). -simulation of the water storage during 8 years in a small area of the Orgeval experimental basin of the C~GREF, -coupled simulation of the water storage in the unsaturated zone and of the piezometric level on the Lille-Sainghin (France) experimental site.

INTRODUCTION

The rainfall-runoff and rainfall-aquifer level lumped hydrologic models are used classically to extend the data in a drainage basin or to estimate the natural recharge of an aquifer by the rain. The Gardenia [START_REF] Roche | SimUlation globale de bassins hydrologiques. Introduction a la modelisation et description du rr~dele GARDENIA[END_REF] hydrological tank model, for example. enables the reliable simUlation, from 4 lumped parameters, of the diSCharge of a stream or a spring at the outlet of a drainage basin. It makes a model possible, frOm 6 lumped parameters. of the level variation at one point of an aquifer recharged by rainfall (Thiery, 1988a), After calibration of the parameters of the model on the piezometric level observations, the groundwater recharge value can be obtained. This kind of model is, in fact. doubly lumped: -lumped for the spatial extent ; as it is considered that a whole basin is described by the same parameters ; -lumped in its working, for the model parameters cannot be measured in the field and convey global phenomena.

The model parameters can therefore, only be determined by calibration on series of observations. Thus, it may happen that several sets of parameters produce almost identical simulations, in a normal period. but risk to give different results in reaction to extreme climatic conditions.

Various applications -in particular [START_REF] Thiery | Analysis of long-duration piezometric records from Burkina-Faso used to determine aquifer recharge[END_REF]Thiery ( . 1988c)-also show that, in the case of slowly varying groundwater, whose storage coefficient is badly knoWll, a satisfactory calibration can sometimes be obtained for several sets of parameters if data of surface run-off on the drainage basin are not available, evapotranspiration with the saturation state was supposed :

(5) The upper soil layer of a basin generally gives rise to a "retention" phenomenon (partly due to interception by vegetation and in the depressions at the surface of the soil). This possible phenomenon can be simulated by a (low) interception capacity. When the possible interception capacity is saturated. the run-off has been chosen proportional to the water content of the soil

The model parameters can thus be predetermined approximately and adjusted so as to be able to reproduce simultaneously the recording of variations of piezometric levels in the groundwater, or of discharge at the basin outlet, and those of water storage (or of average water content) in a given width of the unsaturated zone.

Models like that described by [START_REF] Rambal | Evolution de 1 'occupation des terres et ressources en eau en region mediterraneenne karstique[END_REF] or by [START_REF] Johansson | Methodes for estimation of natural groundwater recharge directly from precipitations -comparative studies in sandy ti11-NATO Advances Research Workshop on Estimation of natural recharge of groundwater[END_REF], also enable modelling of the water storage variations, but the originality of the Gardensol model is that it enables the simultaneous modelling of the eVOlution of the water storage in the unsaturated zone and of the piezometric level at one point in the groundwater.

The model is much more reliable and robust than a classic rainfall-aqUifer level one, which reproduces levels only from a generally badly known storage coefficient. The recharge assessment must be much more reliable, in particular, and the model can be used to extend the data for analysing the resource fluctuation according to the climatic variability.

The diagram for basin simulation is To obtain integrable solutions, very simple laws have been chosen.

ILT-1 J lLT-I J [LT-I]
An unsaturated soil may be described by two characteristic relations of the state of saturation : -a permeability/water content relation, -a suction/water content relation.

Percolation b k

= 1 was chosen, i.e. that the permeability decreases linearly from the Ks value at saturation to a value of 0 for a water content Ok. Moreover a unitary gradient was supposed in all cases. The mean annual balance for the 5 years 1979-1983 can, therefore, be written (in mm/year)

PERC = K = Ks (6 -6k)/(6s -6k) if 6 > 6k (
-u j-----r~---j----T-----t---_,_ '" , '" IVIJ\j~t'\i~)\'v\\t'fAJv1 , '" • ~'" 0 / no m~QScrftm!nt$ ~'"
fiQ 2• SITE Of MEZIERES lES ClERY Exemple 2 -Lille-Sainghin site \," j -----------j -----------, - ~1%1 t"v1hk--I\-----I'i\IJ'o',I-\-",,,t'\-------; ; ~USIl +----''-''---'-,,--+-+-----''II-',jI1-:--+-+ 1; :;; U~C""'-'"'''Y<---C:::-------I-------C:::-------I 1977 \97a
FigUre 2 shows that a very good adjustement is obtained :

correlation coefficient -correlation coefficient storage.

A calibration was made by adjusting the model -at a daily time step-simultaneously on the piezometric series and the water storage series of the layer between 0 to 2.4 m depth. This high permeability value must not be taken as an absolute value, considering the drains which are responsible for the decrease in the water storage.

r ----j ---r ---j ---r ----j --T ---j ----r e 11f".l -" ~10.."'J ! .. 9N ~eN b~"="~b'="*OO~~=O~•O=dbc=c==bc=c==4-==__+--=~'~-J.•-'==__+--=,-
Sainghin IUJ t----j-~T--_j_--t----j---t----r-----. S;>J E control wlibrgl Wiod '" .I---+=~c---l--_
The application shows that the model enables a good simulation of a long series with, in particular, very dry years (1976 had a very low storage) and wet years 0972. 1978, 1979). The simulation is distinctly better than which would be obtained with a single model with a retention capacity that dries up in the same way every summer.

The mean annual balance of the 8 years, 1972-1979, is given in table 3.

The water storage data come from tube nO. 10 of The Gardensol model, which is described in this paper, enables a correct representation of the water content variations of ground exposed to rain, evapotranspiration and natural drainage, It thus enables a simultaneous assessment of the water content variations in the unsaturated zone and the piezometric level variations in the underlying groundwater. Three application examples were presented and a study of robustness showed that the calibration is stable since a calibration carried out over 4 years (at Helarchez-Orgeval) enabled the accurate extrapolation of the water contents over the other 4 years, The model can, therefore, be used for extending the data so as to determine -after calibration over a few years of observations-the resource variability from a long series of climatic data, Other more recent applications (e.g, [START_REF] Bonin | Evaluation par modele semi-global de la recharge d'un aquifere par les pluies a travers la zone non-saturee. Application A un site de stockage en surface de dechets radioactifs de faible et moyenne activite. Collogue international lIydrogeologie et sOrete des dep8ts de dechets radioactifs et industriels toxigues[END_REF] confirmed this model's interest. In the cases, though quite rare, where there are also available data of suction measurements in the unsaturated zone and permeability data in the unsaturated zone, it is surely more interesting to use e physical model according to Richards law like that described by Thiery (1988b) even though it is more delicate and time consuming,

FLINFIL

Fig

  Fig. 1. Gardensol model -Floy diagram

  account of the fact that the suction increases when the water content decreases and partly counterbalances the decreases in permeability, a linear decrease of the actual where PERC ek = percolation = water content corresponding to an almost nil permeability Three application examples are given -coupled simulation of the water storage in the unsaturated zone and of the piezometric level at Mezieres-les-CI~ry (Loiret, France), simulation of the water storage during 8 years in a parcel of the CEHAGREF experimental basin of Orgeval (France), coupled simulation of the water storage in the unsaturated zone and of the piezometric level on the Lille-5ainghin (France) experimental site. Example I -Experimental basin of Mezieresles-Clery (Loiret)

I

  the Boissy-Ie-Chatel parcel(C. Louaa&De, 1984). This drained parcel (drains 0.70 m under the I ground) of 615 m2 which has no appreciable surface I run-off, is part of the M~larchez sub-basin (7 km2). included in the CEMAGREF's experimental1--------------t------t------I----1 basin of the Orgeval. The water contents were -SurfAce zone recorded 200 times by neutron probe down to a . ~:~:~~~~nl:~Interception S S ) b 67 depth of 3.05 m during the 8 years of the


D Thiery

ANALYSIS OF PUMPING TESTS PERFORMED IN A HORIZONTAL RECTAr&ULAR FRACTURE

The analytical solution has been derived for a constant rate uniform flux pumping test in a single horizontal rectangular fracture. The solu tion is given for the drawdown at the pumped well and at an observation well inside or outside the fracture. A set of type curves which depend on two dimensionlesr parameters has been drawn. These type curves give the dimensionless drawdown at the pumped well versus the dimensionless time according to the shape factor of the fracture and the dimensioniess aquifer thickness. Another set gives the average drawdown along the vertical at an observation well versus dimensionless time depending on the dimensionless distance and the shape factor.

For both the drawdown at the pumped well and at an observation well the solution has also been studied for early times and long times in order to derive approximative solutions which are compared to the well known solutions in iso tropic aquifers.

INTRODUCTION

For all these schemes (except the vertical fracture) only the pressure at the pumped well has been computed. These scheaes arc interesting but they are not always adapted to real fractures mainly because of anisotropy of permeability in the horizontal directions.

A horizontal fracture is more likely to be elliptical than circular. In order to analyse drawdowns in such a fracture we derived analytical solutions for a pumping (or injection) test in a rectangular fracture which is a closer npproximation to an elliptical fracture.

The drawdown has been calculated at the center of the fracture but also at an observation well.

DESCRIPTION OF THE SYSTEM

The geometry of the system is defined as fol low :

the aquifer has an infinite lateral extension; its characteristics are :

-constant thickness h -constant Btorage coefficient S (specific storage coefficient S -S/h) -anisotropic permeability : K;:, Ky, Kz.

The fracture , which is situated In the middle of the aquifer, has a rectangular shape with :

-a length 2 xf -a width 2 yf -a negligible thickness

BOUNDARY CONDITIO INSIDE THE FRACTURE

Two types of boundary conditions may be set inside the fracture :

-an infinite conductivity fracture -a uniform flux fra'.ture.

a) An infinite conductivity fracture

It is a fracture which has an equivalent hy draulic conductivity (or permeability) much higher than the conductivity of the surrounding formation. As a matter of fact its transmissivity (K.h pre duct) is much higher then the rock transraissivity This condition is usually met in thick fractures in a low permeability formation. The head gradient in an infinite conductivity fracture is equal to aero.

b) A uniform flux fracture

The discharge from the aquifer to Che frac ture is performed at a uniform rate per unit area of fracture. This boundary condition is probably not exactly fulfilled but, as will be shown later, it is a close approximation of the infinite con ductivity boundary condition.

The uniform flux boundary condition is very interesting because it makes it possible to com pute the drawdowns analytically by the "source functions method" . This uethod has been descri bed by (GRINGARTEN A.C. and RAMEY H.J. 1973) 6 .

The source functioi associated to the rec tangular fracture is obtained as the product of 3 elementary source functions as shown below.

'r*7

The product is then simplified considering that the thickness zf of the fracture is negligi ble compared to the thickness h of the aquifer-

c) Qimenslonless notations

The following dimensionless variables are used :

-coordinates : Nota : All these notation? are consistent with Theis's notations.

Xp = x/xf
The shape factor F is the ratio of the length of the fracture to its width (corrected fay the ho rizontal anisotropy factor) For an horizontally isotropic aquifer : F -I applies to a square fracture (it is the shape which is the closest to the circular fractu re) F • 10 applies to a rectangular 'vacture with a length eqcal to 10 times its width.

The dimensionless thickness h_ is the ratio of the aquifer width to the fracture extension (corrected by the vertical anisotropy factor) .The fracture extension is the geometrical mean of its length and its width.

For a fully isotropic aquifer : h_ = Q.l applies to a fracture with an exten sion equal to 10 times the thickness of the aqui fer h_ -10 applies to an aquifer of thickness equal to 10 times the extension of the fracture.

DRAWDOWN AT THE CENTER OF THE FRACTURE

The drawdown at the center of the fracture is given by the following expressions which are easy LO compute numerically : This equation corresponds to a vertical flow from the formation to the fracture. It does not depend on the shape factor F because there is no flow in the horizontal directions.

For late dimensionless times the drawdown is a linear function of the logarithm of th« time according to the following equation : Equations (It and (2J refer to the center of the fracture, but as a matter of fact, it is a clo se approximation of the drawdown at any location Inside the fracture, because as it will be shown later the drawdown inside the fracture is nearly equal everywhere.

DRAWDOWN AT AN OBSERVATION HELL

f difficulty arises because the drawdown field is 3 dimensionnal : the drawdown is diffe rent at each altitude of an observation well. This difficulty has been solved considering a small diameter observation well with perforations along the whole aquifer thickness. The drawdown measured by such in obsei 'ation well is the average of the drawdown along the vertical direction. This avera ge , which has been computed analytically does not depend on the dimensionless thickness h .

The average drawdown at an observation wel1 situated along the ox axis is given by the follo wing expression s _ J_ 1 -i

•erf ---dt r >> . V /T

The drawdown at an observation we!I situated aloni the oy axis is given by the sumt-expression after replacing F by 1/F. Expression ( r >) has been calculated numerically for the following values of the paraiaetei s : xp = 1, 1.1, 1.2, 1.3, 1.4, 1.5, 2, 3, 10 F =0.1, 0.5, 1, 2, 5, 10, r )0 Two examples of type curves are drawn on graph 3 and 4. They show that for F > ! (observation wells situated along the ox axis), the curves are very close to Theis's solution (x = ">) as soon as X> 2 even for early dimensionless times. For 1a-'_'• dimensionless tines the flow is radial and is described by equation ( 4) rallying JAfOB':; .ipproximation .

DRAHDOMN AT AN OBSERVATION HELL IN CONTACT WITH THE FHACTURE

Th* average of the drawdown along the vertical hat bean computed because it doesn't depend on o any more. The calculations have been performed for various dinensionless tines.

For early dinensionless times the drawdown is given by the following equations :

(6) (7)
The initial drawdown averaged along the vertical is linear with time s it is represented by a straight line of slope 1 in log-log coordinates.

Equations ( 6) and ( 7) are valid, with ,i maxi mum error of IQt , before the following value of dlmensionless time is reached :

Equation (•!) is

1

It appears that JACOB'S approxiiuition applies as soon as the dimensionless time is greater than abou*-10 (as in a pumping test without any frac ture) excepted for x « 1, i.e., at the contact of the fracture, where this minimal time is increa sed to 13 to 17. exnression (8) shows that the linear drawdown is observed during a longer time when the shape factor is smaller i.e. when the shaoe of the frac ture is closer to a square.

For late dimensionless times the variation of the drawdown is given by equation 14) . This equa tion , which describes JACOB'S approximation, is valid as soon as the following dimensionless time is reached -107-It. appears Ui.it JACOB'S approximations is v.ili'l only after •* long (Hnensionless time (grea ter thm> I'i')) for very rectangular fractures i.e. iructutcH wits u large F shape factor. The draw down inside the fracture* has been drawn for diffe rent dimension less:; times on graph c >. This graph stiiiws th.it th*-'tcvmiovm is approximative^ cons tant ni'side the f rarture. Tt moans that the erjuiv.tl Tit tr<ii>K'niusi vity of tl e fracture is very h i gh. Drawdown inside a horizontal rectan gular fracture (average along tlio aquifer width) Trie drawdown corresponding to a pumnin^ test porjorrnivi in a ir.^inite transmissivity horizontal rectangular fracture is then closely approximated by the formula that we have derived for a uniform flux horizontal rectangular fracture.

^vX D' 0 1 R.' ' 0.6 ! 1 ' 1 ' 5.2 ' 5

?. CONCLUSION

A new scheme has been derived for pumping tests in fractured aquifers. This scheme is cha racterized by a single horizontal flat rectangular fracture.

The drawdown has been computed at the center of the fracture and also at an observation well situated along the principal directions. Types cur ves have been drawn which makes it possible to perform pumping test analysis or to compute the drawdown corresponding to a known rectangular

  /Q with T -h VfcxKy " transmis-' geometry factors F -(xf.'yf) VTTK (shape factor) hD«(h/ Vx^). V^TKZ (dimensionle*J thicknear. of the aquifer) -time tap -nTt/x-yF S (relative fo the fracture area)t_x -4 K ht/x 2 S (relative to the distance to the center of the frac Cure).

  o ...r lt ,t «p ( .afg2,] « T Expression (1) is easier to compute for lai values of dimensionless time and expression (2) ir early values. The calculations have been performed for ,e following values of the geometry factors : = 0.1, 0.5, 1, 1.5, 2, 5, 10, 50. 1, 2, 5, 10, 50 which corresponds also to I. ,0.2 0.1, (1 02 . Two kinds of type curves have been drawn : -type curves for constant shape factor F (see graph 1) -typo curvas for constant dlmensionless thic)cneii h <sse graph 2) 10* W 10* Graph 2 -Drawdown at the center of a horizontal rectangular fracture For early dimensionless times the curves show a typical "one half slope" of equation

Hquation

  is the drawdown at time t_, This equation shows that, after a time t the variation of the drawdown does net depend on the geometrical factors F and h any more.Graph I -Drawdown at the center of a horizontal "~ rectangular fractureIt is a radial flov/ which agrees with JACOB'S approximation. The drawdown variation is plotted as a. straight Line on Keni-'ogarithmic paper but Th-.'is curve can not bo recognized on log-log pap3r because of the constant s

  Graph 3 ~ Drawdown at an observation well (hori zontal rectangular fracture)

  CL.), FEUGA (B.l, NOYER (M.L.), TI1IERY (D.

  GRINGARTEN (A.C.), RAMEY (H.J.I. 1974 -Unstea dy state pressure distribution created by a well with a single infinite conductivity ver tical fracture . -Soc. Pet. Eng. J. (Aug. 74) 3. GRINGARTEN {A.C.J, RAMEY (H.J.). -1974 -Uns teady state pressure distribution created by a well with a single lion^ontal fra.T.i-e, par tly ;•.-•• -*r*i* .o-' . r msLiicted entry. -Socie-c.-of pr-ToiPUi' r:ig?nec:s journal -(Augubt 74) 4. CJNCO, RAMEY (H.J.), MILLER . -1975 -Unsteady state pressure distribution created by a well with an inclined fracture. Paper SPE 5591 pre sented at SPE AIHE 50th. Annual fall technical conference wc" exhibition -(Dallas, Sept-Oct. 1975). 5. RAGHAVAN (R.), URAIET (A.), THOMAS (G.W.). 197£ Verfical fracture height effect on transient flow behavior. Paper SPE 6016 presented at the SPE AIHE 51th. Annual fall technical conference and exhibition (New-0rl6ans, Oct. 1976). 6. GRINGARTEN (A.C.), RAMEY (H.J.). -1973 -The use of source and Green's functions in solving unsteady flow problems in reservoirs -Soc. Pet. Eng. J. (0?t. 1973).

	). -1980 -Roches chaudes haute tem
	perature ("Hot Dry Rocks") Contribution a la
	r.i^thodologie de la determination des caracte'ris-
	tiques hydrauliques des milieux rocheux fractu
	res naturellement ou artificiellement -Commis
	sion des Communaut£s EuropSenneo -Contrat
	n° 563-78 EGF {rapport BRGM 80 SGN 029 GEG).
	2.