03/02/2010

Modélisation des
réservoirs souterrains
A. Dassargues
[fichier ModRes2 2010]

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

Rappel et extension des

équations d’écoulement en
„ hypothèses et concept d’EVR milieu saturé
„ hydrostatique
„ porosités
„ conductivité hydraulique et loi de Darcy
„ transmissivité
„ écoulement permanent en milieu poreux
saturé
„ coefficients d’emmagasinements
„ écoulement transitoire en milieu poreux
saturé
„ Conditions aux frontières

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

1
 03/02/2010

Hypothèses et concept d’EVR

l’ EVR (REV) est le volume considéré de milieu
souterrain pour lequel les propriétés vont être
quantifiées (par des valeurs moyennes, équivalentes)

„ assez grand pour être au delà de l’échelle microscopique

„ assez petit pour éviter des lissages nuisibles à la
représentation du processus

… ce concept sous-entend implicitement que le

milieu est continu (et poreux)

„ l’ EVR dépend du problème étudié et des objectifs de

l’étude
„ l’ EVR est utilisé pour les écoulements souterrains
mais aussi pour le transport de contaminants
… ou tout autre ‘quantification’ dans un milieu
souterrain
Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

Hypothèses, concept d’EVR et échelle du problème

Echelle

Micro Dynamique des fluides et étude des pores

Homogénéisation

Macro Essais de laboratoire (du dm à quelques m)

Homogénéisation
Modèles numériques détaillés – Simulation
Méga ‘physiquement significative’ du réservoir

Homogénéisation
Modèles de bilans, modèles linéaires,
Giga fonctions de transfert, etc.

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

2
 03/02/2010

Hydrostatique

charge hydraulique, hauteur piézométrique

Φ tot p
h= =z +
g ρ. g
Profondeur jusqu’à la nappe

Charge de pression Hauteur piézométrique

p p
hp = h= z+
ρ .g ρ .g

Elévation du point P
Point considéré
hg = z
Plan de référence

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

Hydrostatique
… lien direct entre hauteur piézométrique h et
pression interstitielle p: p = ( h − z ) ρ .g

∂p ∂h
= ρ .g
∂t ∂t
pour les problèmes d’écoulements souterrains, la pression
et/ou la hauteur piézométrique
variable principale

Les hauteurs piézométriques sont comparables si et

seulement si: même température
même contenu en sel

si ce n’est pas le cas, … la densité de l’eau va varier et les

mesures de hauteurs piézométriques ne sont plus
comparables : il faut une mesure de la salinité et
faire une correction Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

3
 03/02/2010

Porosité

… la porosité totale est constituée de deux composantes:

nc + S r = n avec nc =
Vm
Sr =
Vim
Vt Vt
porosité de drainage (‘drainage porosity’)
= porosité efficace
ne ≅nc
l’eau qui peut se libérer par drainage gravitaire
(également eau mobile ou eau libre) ….
également = “specific yield” S y
50 %

porosité totale
40 %

30 %

20 %
porosité efficace
= ‘specific yield’
10 %
capacité de rétention
0
diamètre moyen des grains
0.0001 0.001 0.01 0.1 1 10 100 mm Hydrogéologie et Géologie
Terug naar eerste pagina
argile fine argile silt sable fin sable gravier fin gravier grossier de l’Environnement

Conductivité hydraulique et loi de Darcy

… loi expérimentale

quantité d’eau par unité de temps à travers un milieu

poreux: ∆h
Q = K . A.
L

K le coefficient de perméabilité, la conductivité

hydraulique, perméabilité à l’eau (par abus de
language: perméabilité) du milieu poreux (m/s)

Q
… le débit spécifique: q=
A
en m3/(m2.s) donc en m/s Hydrogéologie et Géologie
Terug naar eerste pagina
de l’Environnement

4
 03/02/2010

Conductivité hydraulique et loi de Darcy

Débit spécifique improprement appelé ‘vitesse de Darcy’

… il ne s’agit que d’un débit Q divisé par une surface A

cette surface n’est pas la section

réelle d’écoulement

la section réelle d’écoulement est : [Link]

… pour obtenir une valeur moyenne (sur l’ EVR) de la
vitesse d’écoulement:
q K ∆h
ve = = . m/s ‘ vitesse d’advection’
ne ne L (vitesse effective)
Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

Conductivité hydraulique et perméabilité intrinsèque

K dépend :
„ des propriétés du fluide concerné par
les écoulements: la viscosité
le poid spécifique
„ des propriétés du milieu poreux où
l’écoulement a lieu:
granulométrie, forme des grains,
répartition et forme des pores,
porosité intergranulaire

perméabilité intrinsèque/ perméabilité (m2)

masse spécifique du fluide (kg/m3)
k .ρ.g
K= accélération de la pesanteur (m/s2)
µ
viscosité dynamique (kg/(m.s), N.s/m2 ou Pa/s
Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

5
 03/02/2010

Généralisation de la loi de Darcy

… gradient de la hauteur piézométrique: ⎛ ∂h ∂h ∂h ⎞

grad h = ⎜⎜ , , ⎟⎟
⎝ ∂x ∂y ∂z ⎠
dans un EVR isotrope, la conductivité hydraulique est un scalaire

mais K peut prendre différentes valeurs en fonction de la

direction considérée
la conductivité hydraulique et
la perméabilité intrinsèque
sont donc décrites par des tenseurs: K et k
q = − K grad h cote du point considéré
(par rapport au plan de référence)
k
q = − ( grad p + ρ g grad z )
µ
poids spécifique (constant dans l’EVR)
Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

Transport et application de la loi de Darcy

Exemple: calcul d’un temps de transport par advection pure

q K ∆h
ve = = . m/s
ne ne L

‘ vitesse d’advection’
(vitesse effective)
(Feflow ©, 2002)
Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

6
 03/02/2010

Conductivité hydraulique
valeurs équivalentes, valeurs moyennes

milieu poreux ‘uniformément’ hétérogène

propriété log K distribution normale
„ valeur moyenne (équivalente sur l’EVR) =
n moyenne géométrique des K
K mg = n
∏K
i =1
i

valable également en conditions anisotropes

applications: beaucoup de mesures nécessaires

ne pas oublier les structures
géologiques
exemple: milieu stratifié horizontal
moyenne arithmétique
Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

Transmissivité

… pour un aquifère captif

transmissivité
(m2/s) en un point

e (x , y )
T (x , y ) = ∫ K ( x , y ).dz
0

valeur moyenne de la
épaisseur de l’aquifère conductivité hydraulique
captif en ce point sur la verticale en ce point

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

7
 03/02/2010

Transmissivité

… pour un aquifère libre

z1
y z2
h(x,y)
z3
x

l’épaisseur saturée dans l’aquifère libre

au point de coordonnées horizontales x et y
dépend de la
h (x , y ) hauteur
T (x , y ) = ∫ K ( x , y ).dz piézométrique !
0

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

Ecoulement permanent en milieu poreux saturé

… principe de conservation de la masse: entrée = sortie

− div ( ρ q ) − ρ q = 0
z
le débit
spécifique (m/s) qz
débit de sollicitation extérieure
(‘sink/source flow rate’) par q
unité de volume (s-1),
positif pour pompages, etc.
y
et négatif pour infiltration,
injection, etc. x qy

qx

⎛ ∂q ∂q y ∂q z ⎞
divq = ⎜⎜ x + + ⎟⎟
⎝ ∂ x ∂ y ∂ z ⎠
Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

8
 03/02/2010

Equation d’écoulement en régime permanent

− div ( ρ q ) − ρ q = 0

… en utilisant la loi de Darcy:

[ ( )]
div ρ. K grad h − ρ q = 0
∂ ⎛⎜ ∂h ⎞⎟
K ij − qi = 0
∂xi ⎝ ∂x j ⎟⎠
⎜ h variable principale

⎛k
(
div⎜⎜ grad p + ρ g grad z )⎞⎟⎟ − q = 0
⎝µ ⎠

p variable principale

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

Coefficient d’emmagasinement spécifique

… en écoulement transitoire,
∂ (n ρ )
variation de l’emmagasinement en fonction du temps:
∂t
coefficient d’emmagasinement spécifique
(‘specific storage coefficient’) (m-1)

∂ (n ρ ) ∂h ∂h
= ρ 2 g (α + nβ s + nβ w ) = ρ S s
∂t ∂t ∂t
compressibilité
volumique du compressibilié
milieu poreux de l’eau (Pa-1)
(Pa-1) compressibilité
des grains solides
(Pa-1)
Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

9
 03/02/2010

Coefficient d’emmagasinement spécifique

… souvent, on peut négliger l’influence de la compressibilité

de l’eau et de grains solides par rapport à la compressibilité
volumique du milieu

Ss = ρ g α
… cette relation entre la compressibilité volumique et le
coefficient d’emmagasinement spécifique, démontre
l’existence du couplage direct entre les écoulements non
stationaires et la géomécanique en milieu souterrain
compressible

la compressibilité volumique est dépendante de

la variation de la contrainte effective
la contrainte effective de préconsolidation du
milieu
Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

Ecoulement transitoire en milieu saturé

… principe de conservation
de la masse: entrée = sortie + variation d’emmagasinement

∂ (n ρ ) ∂h
− div ( ρ q ) − ρ q = = ρ. S s
∂t ∂t
∂h
( )
div K grad h − q = S s
∂t
en conditions
strictement
saturées
∂ ⎛⎜ ∂h ⎞⎟ ∂h
K ij . − qi = Ss
⎜
∂xi ⎝ ⎟
∂x j ⎠ ∂t

⎛k
(
div⎜⎜ grad p + ρ.g grad z )⎞⎟⎟ − q = S ∂ h ⎛ S s ⎞ ∂p
=⎜ ⎟.
∂t ⎜⎝ ρ.g ⎟⎠ ∂t
⎝µ
s
⎠

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

10
 03/02/2010

Coefficient d’emmagasinement
… le volume d’eau (m3) libéré ou stocké par unité de
surface de l’aquifère (m2) pour une variation
unitaire de la hauteur piézométrique (m)
… intégration sur la verticale
e( x, y )
S ( x, y ) = ∫ S s ( x, y ).dz nappe captive
0

S = S s .e
h
S = ne + ∫ S s .dz nappe libre
z1

la composante la plus importante de l’emmagasinement

est due au drainage du milieu poreux qui passe
de l’état saturé à non saturé (ou vice-versa)
Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

Coefficient d’emmagasinement
nappe libre nappe captive
S = ne + S s .h
S ≅ ne

plan de référence
= base de la
formation
aquifère

S = S s .e

drainage expulsion
Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

11
 03/02/2010

Equations d’écoulement
∂h
nappe captive ( )
div K grad h − q = S s
∂t
(3D, général)

∂ ⎛⎜ ∂h ⎞⎟ ∂h
K ij . − qi = Ss (3D, forme tensorielle)
⎜
∂xi ⎝ ⎟
∂x j ⎠ ∂t

∂ ⎛ ∂h ⎞ ∂ ⎛ ∂h ⎞ ∂ ⎛ ∂h ⎞ ∂h
⎜ K xx ⎟ + ⎜⎜ K yy ⎟⎟ + ⎜ K zz ⎟ − q = S s
∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂z ⎝ ∂z ⎠ ∂t
(3D, anisotropie)
∂h
div (T grad h ) − q' = S (2D horizontal)
∂t
∂ ⎛⎜ ∂h ⎞⎟ ∂h
Tij − q'i = S (2D horizontal, forme tensorielle)
⎜ ⎟
∂xi ⎝ ∂x j ⎠ ∂t
∂ ⎛ ∂h ⎞ ∂ ⎛ ∂h ⎞ ∂h (2D horizontal,
⎜ Txx ⎟ + ⎜ Tyy ⎟ − q' = S
∂x ⎝ ∂x ⎠ ∂y ⎜⎝ ∂y ⎟⎠ ∂t anisotropie)
Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

Equations d’écoulement

nappe libre … en 3D il est impossible d’écrire les équations

d’écoulement sans tenir compte des processus
ayant lieu dans la zone non saturée

( )
div T (h)grad h − q' = ne.
∂h
∂t
(2D horizontal)

le tenseur de 2ème ordre de la transmissivité

dont les composantes dépendent de la
hauteur piézométrique (équation non lineaire)
∂ ⎛⎜ ∂h ⎞⎟ ∂h
Tij (h) − q'i = ne (2D horizontal, forme tensorielle)
⎜
∂xi ⎝ ⎟
∂x j ⎠ ∂t

∂ ⎛ ∂h ⎞ ∂ ⎛ ∂h ⎞ ∂h
⎜ Txx (h) ⎟ + ⎜⎜ Tyy (h) ⎟⎟ − q' = ne (2D horizontal,
∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂t
anisotropie)

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

12
 03/02/2010

Equations d’écoulement : autres formulations

La notation vectorielle est souvent utilisée, …

avec l’opérateur Nabla ∇, vecteur défini par
(∂/∂x, ∂/∂y,∂/∂z)
Si on fait le produit scalaire de 2 vecteurs, on obtient:
(2D horizontal)

∂ ∂ ∂ ∂J ∂J y ∂J z
∇⋅J = ( , , ) ⋅ (J x , J y , J z ) = x + +
∂x ∂y ∂z ∂x ∂y ∂z
On peut donc écrire l'équation d'écoulement:

∂ ( ρθ ) ∂ (ρθ )
−∇⋅J = − ∇ ⋅ ( ρq) =
∂t ∂t

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

L'équation peut aussi être écrite comme

∂ (ρθ ) ∂ ( ρ qx ) ∂ ( ρ q y ) ∂ ( ρ qz ) ∂ ( ρθ )
− div ( ρq ) = ou − − − =
∂t ∂x ∂y ∂z ∂t
On doit ensuite définir l'équation en terme d'un paramètre
mesurable, qui est
la charge hydraulique h [L]. ∂h
qi = − K ij pour i,j=x,y,z
Utilise la loi de Darcy ∂x j

où Kij est le tenseur de conductivité hydraulique [L T-1]

En 3D, on a
⎛ qx ⎞ ⎡ K xx K xy K xz ⎤⎛ ∂h / ∂x ⎞
⎜ ⎟ ⎢ ⎥⎜ ⎟
⎜ q y ⎟ = − ⎢ K yx K yy K yz ⎥⎜ ∂h / ∂y ⎟
⎜q ⎟ ⎢ K zx K zz ⎥⎦⎜⎝ ∂h / ∂z ⎟⎠
avec ⎝ z⎠ ⎣ K zy

⎛ ∂h ∂h ∂h ⎞
q x = −⎜⎜ K xx + K xy + K xz ⎟⎟ même définition pour
⎝ ∂x ∂y ∂z ⎠ q y , qz Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

13
 03/02/2010

L'emmagasinement est approximé par

∂ ( ρθ ) ∂h
≈ ρS s
∂t ∂t
où Ss est le coefficient d'emmagasinement
spécifique [L-1] défini par
S s = ρg (α + θβ )
et où α et β sont, respectivement, la compressibilité
du milieu poreux et du fluide [L M-1 T2 ]
→ unités inverses de la pression (superficie/force)

− dVmp / Vmp − dVfluide / Vfluide

α= ; β=
dσ e dP
Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

On modifie l'équation

∂ ( ρq x ) ∂ ( ρq y ) ∂ ( ρq y ) ∂ ( ρθ )
− − − =
∂x ∂y ∂z ∂t

On suppose que la densité du fluide est constante

et on obtient ⎛∂ ⎞
⎜ K ij ∂h ⎟ = S s ∂h
∂xi ⎜ ∂x j ⎟⎠ ∂t
⎝
Si les directions principales de K sont alignées
avec les coordonnées x, y, z, on obtient

∂ ⎛ ∂h ⎞ ∂ ⎛ ∂h ⎞ ∂ ⎛ ∂h ⎞ ∂h
⎜ K xx ⎟ + ⎜⎜ K yy ⎟⎟ + ⎜ K zz ⎟ = S s
∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂z ⎝ ∂z ⎠ ∂t

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

14
 03/02/2010

Écoulement 2D horizontal - confiné

Intègre verticalement l'équation 3D sur l'épaisseur de
l'aquifère confiné
∂ ⎛ ∂h ⎞ ∂ ⎛ ∂h ⎞ ∂h
⎜ Txx ⎟ + ⎜⎜ Tyy ⎟⎟ = S −R+L
∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂t

où T : transmissivité [L2 T-1], K × épaisseur

S : coefficient d'emmagasinement [-], Ss × épaisseur
R : taux de recharge [L3 L-2 T-1] = [L T-1]
L : taux de drainance [L T-1]
recharge (R)

écoulement
Épaisseur (b)

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

Exemple de calcul de drainance

Système multiaquifères
Drainance (écoulement vertical) à travers un aquitard
Loi de Darcy
K z' ( hs − h )
L=−
b'
où K'z : conductivité hydraulique verticale de l'aquitard
b' : épaisseur aquitard
hs : charge hydraulique aquifère source
h : charge hydraulique dans aquifère (variable inconnue)
Aquifère
hs
Aquitard b', K'z
h
Aquifère écoulement

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

15
 03/02/2010

Écoulement 2D horizontal - libre

Intègre verticalement l'équation 3D sur épaisseur de
l'aquifère non confiné
Équation non linéaire - Dupuit

∂ ⎛ ∂h ⎞ ∂ ⎛ ∂h ⎞ ∂h
⎜ K xx h ⎟ + ⎜⎜ K yy h ⎟⎟ = S y −R
∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂t

où Sy : porosité de drainage [-]

recharge (R)

écoulement
Épaisseur saturée = h

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

Écoulement 2D vertical (confiné)

Intègre l'équation 3D sur largeur horizontale de

l'aquifère confiné
∂ ⎛ ∂h ⎞ ∂ ⎛ ∂h ⎞ ∂h
⎜ K xx ⎟ + ⎜ K zz ⎟ = S s − R'
∂x ⎝ ∂x ⎠ ∂z ⎝ ∂z ⎠ ∂t

où R' : taux de recharge par largeur unitaire [T-1]

recharge (R')

épaisseur
unitaire

écoulement

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

16
 03/02/2010

Conditions aux frontières

Modèles d’écoulement

„ Conditions de Dirichlet ou de hauteur

piézométrique imposée
„ Conditions de Neumann ou de flux imposé
„ Conditions de Cauchy ou de flux dépendant
d’une hauteur piézométrique

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

Hauteur piézométrique imposée

(condition de Dirichlet)
La hauteur piézométrique est spécifiée sur la frontière:

h(x, y, z, t) = f ' (x, y, z, t)

f ' peut varier dans l’espace et le temps
(une valeur par nœud concerné et par pas de temps)

le programme calcule alors un flux en chaque nœud

concerné

mathématiquement, pour que le problème soit défini

de façon univoque il faut minimum une valeur de
hauteur piézométrique imposée dans le domaine
simulé
Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

17
 03/02/2010

Hauteurs
piézométriques Hauteurs
imposées piézométriques
imposées

Hauteurs
piézométriques
imposées

Limites de la zone où
les résultats sont
significatifs Hydrogéologie et Géologie
Terug naar eerste pagina
de l’Environnement

Flux imposé (condition de Neumann)

La dérivée première de la hauteur piézométrique est
spécifiée sur la frontière concernée:
∂h
( x, y, z, t ) = f ' ' ( x, y, z, t )
∂n
f ' ' le gradient piézométrique normal à la frontière concernée,
sa valeur peut varier dans l’espace et le temps
(une valeur par nœud concerné et par pas de temps)
Par application de la loi de Darcy, c’est une manière
d’imposer un flux à travers la frontière:
∂h
− K. ( x, y , z , t ) = q ' ' ( x, y , z , t )
∂n
q ' ' : flux imposé à travers la frontière (m/s)
cas particulier: f ''= 0
Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

18
 03/02/2010

Flux imposé (2)

si le flux est imposé non nul, la condition de débit à travers
la surface de la frontière peut s’écrire:
∂h
∫ K . ∂n ( x, y, z, t ).dS = Q( x, y, z, t )
S

Q débit d’eau souterraine à travers la frontière (m3/s)

Flux imposé

K << K >>

Flux imposé ?

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

Flux imposé (3)

… autre cas particulier :

flux en provenance des eaux de surface

K'
K'
q ' ' = (hs − b) e'
e'
hs b

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

19
 03/02/2010

Flux dépendant de la hauteur piézométrique

(condition mixte ou de Fourier/Cauchy)

Une combinaison (relation linéaire) de la hauteur piézométrique

et de sa dérivée première est spécifiée sur la frontière:
∂h
a. (x, y, z, t) + b.h(x, y, z, t) = f ' ' ' (x, y, z, t)
∂n
f ' ' 'peut varier dans l’espace et dans le temps
(une valeur par nœud concerné et par pas de temps)

„ interactions entre eaux de surface et eaux

souterraines
„ interactions entre différents aquifères

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

Flux dépendant de la hauteur

K' piézométrique (2)
q' ' =(hs − h)
e' K'
∂h K ' K'
− K . + .h = .hs
e'

∂n e' e'
hs h

K'
K' e'
q' ' = (h1 − h)
e' h h1

Hydrogéologie
Terug et Géologie
naar eerste pagina
de l’Environnement

20