Cours Modelisation [PDF]

Zitiervorschau

Modélisation hydrogéologique Master 2 HSE

1 NOTIONS DE BASE Voici quelques notions de base qui interviennent dans la description des écoulements en milieux poreux.

1.1 Volume Elémentaire Représentatif (VER) Un volume élémentaire représentatif est un volume pour lequel les propriétés caractéristiques moyennes (comme la porosité, la perméabilité dans le cas d’un milieu poreux) peuvent être déduites. En réalité un milieu poreux est constitué de graines solides et vides pour lesquelles il n’est pas possible d’attribuer des notions comme la porosité et la perméabilité qu’à partir d’une échelle supérieure de plusieurs ordres de grandeur à l’échelle des pores. Le choix du VER doit donc répondre aux critères suivants (de Marsily 1994) : • Le VER doit contenir un grand nombre de pores afin d’avoir une moyenne globale significative ; • Le VER doit être suffisamment petit pour que les variations des propriétés d’un domaine au domaine voisin puissent être approchées par des fonctions continues pour pouvoir introduire l’analyse infinitésimale, sans introduire d’erreur décelable par les instruments de mesure à l’échelle macroscopique. D’après les critères ci-dessus, un VER dépend non seulement de la structure du milieu poreux, mais aussi des phénomènes physiques étudiés. Un VER doit être assez grand pour représenter la structure du milieu poreux, mais aussi petit pour que les variations des propriétés, parfois nonlinéaires (teneur en eau), soient continues. Une telle définition appliquée à l’hydrogéologie, est sûrement subjective car l’hétérogénéité existe à toutes les échelles d’un milieu poreux naturel, et plusieurs hypothèses de modélisation existent pour chaque problème.

1.2 Définitions • Porosité φ : rapport du volume des vides au volume total Vt du sol

• Teneur en eau volumique θ (m3 /m3) : rapport du volume de l'eau Veau contenue dans les pores (ou vides) du sol, au volume total Vt du sol considéré :

Saturation S : rapport du volume d'eau au volume des vides dans le sol :

• La conductivité hydraulique et la perméabilité

La loi de Darcy n'est strictement applicable que pour des milieux homogènes où l'écoulement de l'eau est laminaire. Elle ne peut être utilisée en particulier pour les réseaux karstiques. 1

Modélisation hydrogéologique Master 2 HSE Le coefficient de perméabilité est propre à chaque réservoir; il dépend notamment de la porosité efficace et de la viscosité du fluide; il augmente avec la profondeur (l'augmentation de température diminue la viscosité). 

Extension aux différents liquides et aux grandes profondeurs (pression et température)

Le coefficient de perméabilité doit tenir compte de toutes les caractéristiques du milieu poreux (structure du réservoir) et de l’eau qui le traverse (viscosité dynamique et poids volumique). Ce coefficient de perméabilité au sens large ne correspond plus à la conception initiale d'Henri Darcy valable pour les eaux souterraines normales de faibles profondeurs. K. Hubbert a développé en 1969 une expression généralisée applicable dans tous les cas en différenciant le réservoir du fluide.

Figure 8: Principe d'écoulement d'un fluide dans un réservoir

V  N.d102 .

g .i k.i 

2 avec kint perméabilité intrinsèque : k int  N .d10 en m² ou en Darcy

d’où

k  k int .       

g 

v : vitesse d’écoulement en m/s ;  : viscosité dynamique du liquide en centipoise (=  x viscosité cinématique) ;  : masse volumique ; N : facteur de forme sans dimension (100 en moyenne) ; i : gradient hydraulique ; d10 : diamètre efficace des grains en cm ; g : force de gravité.

2

Modélisation hydrogéologique Master 2 HSE 

La perméabilité (Kint) on la dénome aussi perméabilité intrinsèque ou spécifique car elle ne dépend que de la structure et de la connectivité des pores et non du fluide se déplacant dans les pores. k est exprimée en m2 ou cm2 est très faible d'où l'emploi du Darcy; 1 Darcy vaut 10-8 cm2 et 10-12 m2. On peut estimer k à partir de la granulométrie d'un milieu poreux (équation de Hazen):

k = 6.54 10-4 d102 Il est nécessaire de rechercher les paramètres spécifiques du fluide lorsque que l’on sort de l’étude des nappes superficielles classiques. Ces derniers évoluent en fonction de sa viscosité (hydrocarbures par exemple), de sa pression et de sa température (grandes profondeurs). Pour les réservoirs, les propriétés propres seront ainsi exprimées par sa perméabilité intrinsèque. 

Transmissivité (T)

Pour les écoulements en milieux poreux, le paramètre d’importance est la conductivité hydraulique. Dans le cas des aquifères, il est souvent préférable de faire usage de la transmissivité (T). Ce paramètre est la mesure de la quantité d'eau qui peut être transmise horizontalement sur toute l'épaisseur d'un aquifère saturé (b), sous un gradient unitaire. Ainsi : T = kb [L2T-1] (habituellement m2/j) Dans le cas des aquifères composés de plusieurs couches, la transmissivité totale est la somme de la transmissivité de chacune des n couches perméables, soit : T=

∑Ti

Ce paramètre est souvent obtenu lors d’essais de pompage.

 Coefficients de stockage La variation en fonction du temps des charges hydrauliques des eaux souterraines conduit à la libération de l’eau emmagasinée dans l’aquifère ou conduit à l’accumulation de l’eau dans cet aquifère. Deux mécanismes peuvent les induire : la libération de l’eau de la nappe libre entraîne le changement du niveau de cette nappe et/ou la libération de l’eau d’une nappe captive qui va entraîner un ajustement de la pression dans le système d’aquifère. (se référer à Bear (1979) et Fetter (2001) pour plus de détails et discussions sur les concepts physiques de base des coefficients de stockage). Les coefficients de stockage principaux sont la porosité de drainage et le coefficient d’emmagasinement. La porosité de drainage Sy est le volume d’eau drainé par unité de surface plan d’un aquifère pour que sa charge hydraulique baisse d’une unité. Elle est sans dimension. Le coefficient d’emmagasinement spécifique Ss est le volume d’eau libéré par unité de volume d’aquifère pour que sa charge hydraulique baisse d’une unité. Sa dimension est L-1. Les gammes typiques de variation de ces deux coefficients de stockage sont : 10-6 ≤ Ss ≤ 10-3 en nappe captive 10-3 ≤ Sy ≤ 0,25 en nappe libre  Coefficient d'emmagasinement (S) 3

Modélisation hydrogéologique Master 2 HSE Suite à des changements du potentiel hydraulique dans un aquifère,l'eau peut être libérée ou emmagasinée. Dans le but de modéliser ce phénomène, l’on définit un paramètre approprié, soit le coefficient d'emmagasinement (S, pour storavity). S est défini comme la quantité d'eau qu'un aquifère relâche ou emmagasine par unité de surface d'aquifère d'épaisseur b (volume b m3), lorsque la charge hydraulique varie d'une unité. C'est un coefficient adimensionnel.  Dans le cas d'un aquifère à nappe captive, l'eau est mobilisée de la façon suivante: une variation de charge engendre une variation de pression d'eau interstitielle, donc une variation de contrainte effective, et donc un tassement (ou expansion) des solides de l'aquifère. La valeur de S est, par conséquent, petite dans un aquifère à nappe captive, typiquement de l'ordre de 10-5 à 10-3.  Dans un aquifère à nappe libre, l'eau est mobilisée par drainage gravitaire à l'intérieur d'un volume unité (surface unité x hauteur unité correspondant à la variation d'une unité de charge). Le coefficient d'emmagasinement S dans un aquifère à nappe libre correspond donc à la notion de porosité efficace (drainage gravitaire). La valeur de S est alors typiquement de l'ordre de 1% à 30%. En conclusion, les deux schéma ci-dessous résume la définition du Coefficient d’emmagasinement et de la porosité

4

Modélisation hydrogéologique Master 2 HSE

 Charge hydraulique, sens de l’écoulement Une particule fluide animée d’une vitesse u , à une pression p et située à l’altitude z, est caractérisée par sa charge h qui est égale à son énergie par unité de masse.

Le premier terme correspond à l’énergie cinétique, les deux autres à l’énergie potentielle. Pour les écoulements souterrains, l’énergie cinétique peut être négligée car les vitesses de circulation sont toujours très faibles. L’équation (ci-dessus) se ramène alors à :

Alors, la charge hydraulique d’un élément de volume d’aquifère est la cote de la surface libre de l'eau dans un tube ouvert dont l'orifice inférieur pénètre cet élément de volume d’aquifère (piézomètre). La charge hydraulique s’écrit :

ou P est la pression,

ρ

la masse volumique du liquide, g l’accélération de la pesanteur, et

z la hauteur au-dessus du niveau de référence (Figure ci-dessous). Toutes les pressions sont relatives à la pression atmosphérique. Lorsque la masse volumique ρ peut être considérée comme constante, la charge hydraulique s'identifie au niveau piézométrique.



Charge hydraulique et cote piézométrique : le cas incompressible 5

Modélisation hydrogéologique Master 2 HSE En un point donné de la conduite, l’énergie d’un fluide incompressible se compose : Ehydraulique=Epotentielle +Epression+Ecinétique 

de son énergie potentielle



de son énergie de pression

, ,

 de son énergie cinétique . L’énergie par unité de volume E est donnée par la relation suivante : Cette énergie par unité de volume est aussi appelée la charge. Or il faut savoir qu’entre deux points A et B de la conduite, il existe une dissipation de l’énergie hydraulique. Cette perte d’énergie (ou perte de charge) est due aux frottements sur la conduite, aux phénomènes turbulents ou encore à des phénomènes locaux tels que les contractions ou les élargissements brusques de la veine liquide.

Ce phénomène se traduit par la relation de Bernoulli : EB = EA-JAB avec :  EA : charge en A,  EB : charge en B,  JAB : charge entre A et B. Dans les métiers de l’eau, il est d’usage d’exprimer la charge en mètres de colonne d’eau [mCE]. On parle alors de charge hydraulique :

Soit :

où les 2 premiers termes de la somme correspondent à la charge statique ou charge piézométrique. Remarque, la cote z d'un point M est défini par rapport à une cote de référence. On définit la charge hydraulique d’un fluide incompressible et soumis à la seule gravité par : 2

h=

V P + +z 2 g ρg

où V est la vitesse au point considéré, et z sa cote. En hydraulique souterraine les vitesses sont V2 lentes, et on peut généralement négliger le facteur de charge dynamique 2 g

6

Modélisation hydrogéologique Master 2 HSE et l’on a donc, pour le fluide, une charge hydraulique équivalente à la charge statique : P h= + z ρg Si l’on creuse un trou, jusqu’au point en question à la cote z, s’élève au-dessus de ce point une colonne d’eau, dont le sommet est en équilibre avec la pression atmosphérique. Le sommet de cette colonne d’eau est à la cote h. Ceci est appelé cote piézométrique du point considéré. Il faut pour faire ce raisonnement avoir supposé que la compressibilité de l’eau est négligeable (sans quoi on ne peut pas définir la charge hydraulique). 

Coefficient de compressibilité

Définition : le coefficient de compressibilité de l'eau β=

Ou

β (beta) est défini par :

−∆ V 1 . V ∆p ∆V V

représente la variation relative de volume provoquée par la variation de pression

∆p Dimension : le coefficient de compressibilité est homogène à l'Inverse d'une pression soit β = (ML-1 T-2)-1 = M-1 L T2 Unités : En Europe C.G.S. : (barye)-1 mais on utilise également des unités non homogènes ; - le (kg/cm2)-l - le (m d'eau)-l = (1/10 )cg/cm2)-l - l'unité de volume/atmosphère ou (atm,)-1  Pressions Définition : En ce qui concerne les écoulements souterrains, on désigne sous le terme pression (p) la pression qui règne dans le fluide contenu dans les pores du milieu poreux. En général on mesure ces pressions par référence à la pression atmosphérique prise égale à 0. Dans ces conditions : - la pression en un point P du milieu poreux s'exprime en fonction de la hauteur d'eau dans un piézomètre par : p=ρ . g .h p

7

Modélisation hydrogéologique Master 2 HSE

Mesure des pressions de fluide en milieu poreux 

Dans le cas des nappes à surface libre la pression en un point du milieu poreux est donc positive au-dessous de la surface libre dans la zone saturée, négative au-dessus de la surface libre dans la frange capillaire et la zone non saturée. Ces pressions "négatives" sont souvent appelées tensions. Leur mesure (évidemment plus délicate que celle des pressions positives) s'effectue à l'aide de tensiomètres. Si on parle en terme de pression absolue on a : Pression absolue = Pression atmosphérique + P Unités Rappelons les principales unités et équivalences (dimension M L-1 T-2 ) - unités homogènes C.G.S. : la barye = 1 dyne/cm2 - unités dérivées le bar = 106 baryes (le bar est la seule unité légale en France) - unités courantes mais non homogènes FRANCE \ l'atmosphère = 7 60 mm de mercure le mm d'eau le kg/cm2 U,S.A. le p.s, i. = pound per square inch  Potentiel hydraulique Définition : on appelle potentiel hydraulique pour les milieux poreux ou charge. en un point le terme : p φ= + z ρg z = cote du point considéré par rapport à un plan horizontal de référence arbitrairement fixé p = pression de pore au point considéré

8

Modélisation hydrogéologique Master 2 HSE

Dans la zone saturée, le potentiel, comme la pression p0, se déduit d'une lecture de niveau dans un piézomètre dès lors que l'on a choisi le plan horizontal de référence. 2. LES EOUATIONS FONDAMENTALES DE LA PHYSIOUE DES MILIEUX POREUX SATURES Les équations du mouvement des fluides dans les milieux poreux dérivent de la combinaison de trois équations fondamentales de la physique des fluides en milieu poreux: 2.1 Equation de continuité Cette équation traduit le principe de la conservation de la masse fluide. Elle s'obtient dans le cas général en faisant la balance de matériel sur un élément infinitésimal cubique.

Figure. 2- Conservation de la masse dans un volume infinitésimal dans la direction X : ∂ ρu ∆ x ρu− Entrant : ∂x 2 Sortant : ∂ ρu ∆ x ρu+ ∂x 2 Dans la direction Y : 9

Modélisation hydrogéologique Master 2 HSE

Entrant :

ρv−

∂ ρv ∆ y ∂y 2

Dans la direction Z : ∂ ρw ∆ z ρw− Entrant : ∂z 2

sortant :

ρv+

sortant :

∂ ρv ∆ y ∂y 2

ρ w+

∂ ρw ∆ z ∂z 2

Masse spécifique ( ρ ¿ Définition ; la masse spécifique ("density" des anglo saxons) d'un corps est la masse de l'unité de volume de ce corps (dans des conditions précisées de température et de pression). Dimensions : M L 3 Unités : En Europe on emploiele plus souvent le g/cm3 (syst, C.G.S.O Aux USA la lb. sec2ft4 (déduite de l'unité de poids spécifique Ib/ft3) Pour l'eau. La masse spécifique de l'eau prise à la pression atmosphérique est assez peu variable en fonction de la température. Elle admet un maximum proche de 1 aux environs de 4 °. Pour cet élément le principe s'énonce physiquement ainsi ; - toute variation du flux massique total algébrique à travers les parois du cube se traduit par une variation égale de la masse emmagasinée dans le cube. Avec les notations suivantes complétées par la figure 2. ⃗v =u i⃗ + v ⃗j+ w k⃗

vitesse du fluide

ρ(rho) = masse spécifique du fluide

ω

(omega)= porosité du milieu poreux

Le principe s'exprime mathématiquement de la façon suivante : ±

Flux entrant

selon x :

[

ρu−

∂ ( ρu ) ∆ x [∆ y.∆ z] ∂x 2

selon y :

[

ρv−

∂ ( ρv ) ∆ y [ ∆x .∆ z] ∂y 2

selon z :

[

ρw−

Flux sortant

]

]

]

∂ ( ρw ) ∆ z [∆ x.∆ y ] ∂z 2

[ [

ρu+

ρv+

[

]

∂ ( ρv ) ∆ y [∆ x.∆ z] ∂y 2

ρw+

10

]

∂ ( ρu ) ∆ x [∆ y .∆ z] ∂x 2

]

∂ ( ρw ) ∆ z [∆ x.∆ y] ∂z 2

Modélisation hydrogéologique Master 2 HSE = variation de la masse d’eau du cube= ∂ ( ρω ) .∆ x.∆ y .∆ z ∂t

Si on admet flux entrant ¿ 0 on trouve : ∂ ( ρu ) ∂ ( ρv ) ∂ ( ρw ) ∂ ( ρω ) + + = ∂x ∂y ∂z ∂t

ρv −∂ ( ρω ) ou en écriture vectorielle : (¿)= ∂ t ¿ ⃗¿

Remarques 1. Soulignons que dans ce qui précède v désigne la vitesse macroscopique du fluide, c'est-àdire un flux par unité de surface, et non la vitesse microscopique réelle des particules fluides 2. On notera que dans certains cas particuliers, il est plus avantageux d'écrire l'équation de conservation en considérant un élément de volume non cubique. A retenir : Le principe de continuité exprime que dans un volume fermé fixe, la variation de la masse de fluide durant une unité de temps est égale à la somme algébrique des flux massiques traversant la surface du volume considéré. C’est le « rien ne se perd, rien ne se crée » de Lavoisier lorsqu’il énonce le principe fondamental de la conservation de la matière. L’équation de continuité s’écrit :

ρ

ε

Où est la masse volumique du fluide [ML-3] est la porosité du milieu poreux définie comme étant le rapport entre le volume des vides

par le volume total. Elle est sans dimension. ⃗ V la vitesse de filtration de l’écoulement (vitesse de Darcy) [LT-1] exprimant la vitesse fictive d'un fluide qui percolerait à travers le milieu en occupant tout l'espace (pores + grains) au lieu ⃗ de n’occuper que les vides. Le flux du vecteur V à travers une surface quelconque de milieu poreux représente ainsi le débit d'eau qui la traverse. En hydrogéologie, il est important de rajouter un terme source à l’équation de continuité Ce terme correspond aux prélèvements (ou apports) d’eau que l’on peut réaliser dans un milieu (forage par exemple). Ce terme source noté w représente le débit 11

Modélisation hydrogéologique Master 2 HSE d’eau prélevé (ou apporté) par unité de volume en chaque point. Le débit massique apporté sera donc ρ w. L’équation de continuité s’écrit alors:

2.2. Equation d'équilibre Elle traduit le principe fondamental de la mécanique qui s'écrit sous sa forme la plus générale : m⃗ Γ =∑ ⃗ F

m = masse ⃗ Γ = accélération(gamma)

∑ ⃗F

= forces agissant sur le système

En mécanique des fluides classiques, on démontre : pour les fluides visqueux compressibles, cette équation se traduit par :

ρ

d ⃗v 1 =−⃗ grad p+ ∑ ⃗ F e + μ Δ ⃗v + μ ⃗ grad θ dt 3

avec :

Fe ∑⃗

(Équation de Navier-Stocks)

: somme des forces extérieures

∂2 ∂2 ∂2 ∆=∇ = 2 + 2 + 2 ∂x ∂ y ∂ z 2

θ=¿ ⃗v

1 ⃗ μ grad θ (pour les fluides incompressibles le terme : 3 =0

Cette équation régit donc l'écoulement des fluides réels dans le cas général. Dans le cas des écoulements en milieu poreux, c'est elle qui régit l'écoule¬ ment mliçros^02iaM§ des fluides dans les canalicules du milieu poreux solide (voir remarque 1 précédemment). Malheureusement, on sait que même dans le cas d'écoulements présentant des conditions aux limites simples, l'équation de NAVIER STOCKES est d'un maniement très difficile. A fortiori, elle est pratiquement inapplicable de façon rigoureuse aux écoulements tortueux qui se forment dans le milieu poreux. Toutefois, certains auteurs (HUBBERT, N.K. BOSE, DE WLESI) ont essayé de déduire de la loi de NAVIER STOCKES assortie de considérations statistiques et de diverses

12

Modélisation hydrogéologique Master 2 HSE simplifications portant sur la nature des écoulements dans le milieu poreux, la loi générale régissant les écoulements macroscopiques dans les milieux poreux. Ces diverses démonstrations sont parfois citées dans les ouvrages fondamentaux (cf. DE WIEST pp. 176-177). Mais d'une part ces "démonstrations" sont critiquables et ne font pas l'unanimité, d'autre part l'aspect macroscopique des écoulements en milieux poreux peut être abordé en première approximation de façon plus simple et plus pratique En effet, on connaît expérimentalement la relation qui régit un écoulement macroscopique linéaire dans un milieu poreux : c'est la loi de Darcy qui sous sa forme originelle s'écrit :

Q=K . A .

H 1−H 2 L

Q ∆H =V =K . A L ou

avec : Q débit volumique A section de l'écoulement Q/A = V vitesse macroscopique (ou moyenne) h1 et h2 charges piézométriques amont et aval L longueur de l'écoulement Cette relation est strictement expérimentale et "vérifiée" uniquement dans les conditions de l'expérience de DARCY (c'est-à-dire, en particulier, pour une dimension). On peut toutefois admettre qu'elle traduit dans les faits une relation plus générale qui s'énoncerait ainsi : "la vitesse macroscopique d'un fluide circulant à travers un milieu poreux est directement proportionnelle aux forces agissant sur le fluide". Si on donne à cette relation la valeur d'un postulat et si on tient compte des relations de dimension, on peut écrire par extension la loi plus générale suivante (ou loi de DARCY généralisée ) On peut écrire l’équation de Darcy sous la forme (equation #): V x=

−k x ∂ P −F x μ ∂x

V y=

−k y ∂ P −F y μ ∂y

V z=

−k z ∂ P −F z μ ∂z

(

(

(

) )

) 13

Modélisation hydrogéologique Master 2 HSE dans laquelle v(vx,vy, vz)= vitesse macroscopique 0 0 =tenseur de perméabilité pour un terrain supposé anisotrope (les kz

‖ ‖

kx 0 k= 0 ky 0 0

axes principaux d 'anisotropie ayant été pris comme axe de référence

P = pression de pore →

F (Fx/ Fy / Fz ) = résultante des forces de volume μ

= viscosité dynamique du fluide

dans le cas particulier ou F dérive d'un potentiel V (F = - grad V), (Equation#) devient k ⃗v = gradϕ avec ϕ =( P+V ) μ Si donc les seules forces de volume qui interviennent sont les forces de pesanteur, on aura : −‖k‖ ⃗v = ρg gradφ μ signe+ si axe ↑ et−siaxe ↓ p avec : φ= ± z ¿ ρg φ

étant par définition le potentiel hydraulique (encore appelé charge hydraulique*)

C'est en général sous cette forme que l'on a à employer la loi de Darcy, Elle constitue pour les écoulements macroscopiques en milieu poreux l'équivalent de la loi de Navier Stockes pour les écoulements fluides en milieu libre. Elle traduit le principe fondamental de la dynamique. C'est elle que nous retiendrons comme 2ème loi fondamentale.

2.2. Equation d'état Pour les fluides, en général, l'équation d'état est une relation de la forme : F (p, ρ , T) = 0 Mais en fait les liquides, et l'eau en particulier, en présence de milieux poreux restent toujours dans des conditions quasi-isotherme et l'équation d'état se réduit à une relation entre les termes p et ρ . F ( p,

ρ ,T)=0

Cette relation s'écrit pour les liquides :

14

Modélisation hydrogéologique Master 2 HSE −β 0( p0− p)

ρ= ρ0 . e

ρ masse spécifique de la pression p ρ0 masse spécifique à la pression p de référence 0 β 0 coefficient de compressibilité à la pression p de référence 0

2.3 La loi de Darcy La formulation vectorielle courante de la loi de Darcy est la suivante :

´ K est la perméabilité du milieu poreux [LT-1] relativement à un fluide donné (l’eau en ´ l’occurrence). On admet que K

Où

est un tenseur dissymétrique :

Lorsque les directions d’anisotropie du milieu coïncident avec les axes ´ orthogonaux utilisés pour décrire le milieu, le tenseur K se réduit à ses composantes ⃗ diagonales. Dans ce cas, la direction du vecteur vitesse de Darcy V est parallèle au gradient de la charge hydraulique ⃗ grad (ℎ). Dans un milieu sédimentaire stratifié, il est évident que les directions parallèles et perpendiculaires à la stratification s'identifient à ces directions privilégiées de l'écoulement. Le modèle numérique tridimensionnel est basé sur l'équation de l'écoulement de l'eau dans un milieu poreux qui s'écrit :

Où : x, y et z sont les coordonnées cartésiennes [L]; 15

Modélisation hydrogéologique Master 2 HSE Kxx, Kyy et Kzz, sont les conductivités hydrauliques du milieu selon x, selon y et selon z [L/T] ∂ est le symbole de dérivée partielle; h est la charge hydraulique produisant l'écoulement, correspondant typiquement à l'élévation du toit de la nappe d'eau au-dessus d'un niveau de référence [L]; t est le temps [T]; P est un flux volumétrique unitaire représentant l'infiltration [1/T]; W est un flux volumétrique unitaire représentant le pompage [1/T]; Ss est le coefficient d'emmagasinement du milieu [1/L].

3. Introduction à la modélisation hydrogéologique 3.1 C’est quoi un Modèle Simplification de la réalité dont le but est de comprendre l’évolution d’un système réel pour pouvoir en prédire l’évolution. La plupart des modèles hydrogéologiques utilisés de nos jours sont des modèles mathématiques déterministes. Ils sont basés sur les principes de la conservation de la masse et leur mise en œuvre nécessite généralement la résolution d’équations différentielles partielles. Des solutions exactes peuvent être obtenues analytiquement (modèle analytique) alors que les méthodes numériques fournissent des solutions approchées à travers la discrétisation du système dans l’espace et le temps (modèle numérique). 3.2Discrétisation Le modèle numérique utilise une solution mathématique des équations d’écoulement et/ou de transport. Le domaine étudié est discrétisé dans l’espace et dans le temps et un schéma numérique est utilisé (Différences Finies, Eléments Finis, Volumes Finis, etc). Ces outils permettent de simuler des conditions hydrogéologiques hétérogènes, des gradients verticaux de charge ou de concentration ainsi que des phénomènes transitoires (pompages, variation d’émission de source de pollution dans l’espace et dans le temps). La Figure suivante montre un exemple de maillage régulier, le système est ici discrétisé en cellules de 5 m de côté. Ce type de modèle peut être employé pour simuler le transport de traceurs non réactifs. Le modèle peut également être utilisé pour simuler le transport réactif de substances organiques et inorganiques, soit en utilisant un modèle Kd, soit en le couplant à un modèle géochimique de spéciation (à ce sujet voir les réserves quant à l’usage de ces modèles dans le guide TRANSPOL « Recommandations pour la modélisation des transferts des éléments traces métalliques dans les sols et les eaux souterraines », Rollin et Quiot, 2006 ; ainsi que dans l’article « Importance

16

Modélisation hydrogéologique Master 2 HSE de la spéciation géochimique pour l’atténuation des pollutions métalliques » Claret et al., 2007). 3.2.1 Différences Finies (DF) Cette méthode de résolution est basée sur la recherche, au centre de chacune des mailles, d’une valeur numérique de la charge hydraulique, supposée représenter une valeur moyenne de la véritable charge. Les mailles considérées peuvent être des carrés, des rectangles ou des parallélépipèdes rectangulaires (en 3D). Cette méthode, facile à comprendre et à programmer, convient très bien aux problèmes régionaux d’écoulement de nappes, mais elle est plus limitée en ce qui concerne la résolution de l’équation de transport et le calcul des concentrations (d’autres méthodes existent comme le TVD, Total Variation Diminishing). Le calcul des charges se fait en établissant l’équation d’équilibre des débits entre la maille de calcul (C) et ses six voisines (haut, bas, sud, nord, est, ouest). 3.2.2 Eléments Finis (EF) Plus délicate à comprendre et à programmer cette approche est néanmoins plus flexible que celle des DF. La forme des mailles est moins limitée (triangle et quadrilatère en 2D, tétraèdre et parallélépipède en 3D). Cependant, la multiplication des données d’entrée induit un temps de calcul plus important et selon la connaissance du site, une incertitude accrue vis à vis des résultats. Pour les écoulements, cette méthode est généralement employée dans le cas de problèmes locaux de génie civil tels que le calcul du débit d’exhaure d’un fond de fouille ou la simulation des écoulements autour d’un barrage, elle s’applique aussi lorsque l’on observe des discontinuités géologiques telles que des fractures ou fissures. Pour le transport, les EF sont bien adaptés puisqu’ils permettent de stabiliser la dispersion numérique plus efficacement qu’avec les DF. Les principales différences entre ces deux méthodes sont donc le temps de calcul et l’échelle des phénomènes naturels à simuler. Dans le cas des DF le maillage est structuré, il est non structuré pour les EF..

17

Modélisation hydrogéologique Master 2 HSE

Représentation du maillage en élements finis et des différents types de limites du modèle(nappe de la plaine des Galets)

3.2.3 Principes de discrétisation On utilise généralement un découpage de l’aquifère à étudier en mailles homogènes de formes variables (figure ci-dessus) dans lesquelles les différents éléments caractéristiques ont des valeurs moyennes : - Transmissivité – T – - Coefficient d’emmagasinement – S – - Débit prélevé ou injecté – Q – - Infiltration par la pluie efficace – Inf – (ou par la maille de la nappe d’au-dessus) - Niveau piézométrique ou charge – H On applique à chaque maille les lois fondamentales de l’hydrodynamique. On mène les calculs d’une maille à l’autre par approximations successives de manière itérative. Partant d’un état initial des charges dans les mailles, on les recalcule les unes après les autres plusieurs fois avec les charges des mailles voisines et les conditions de débit, d’infiltration et de charges imposées dans certaines mailles situées en limite. On arrête les itérations lorsqu’on obtient une quasi stabilisation des charges calculées dans toutes les mailles. Il est possible, au cours du calcul, de tenir compte des variations de transmissivité avec la hauteur d’eau (nappe libre), de l’assèchement ou du débordement en d’autres points (rivières, gravières, etc…). 

Taille et forme des mailles

18

Modélisation hydrogéologique Master 2 HSE

La taille des mailles dépend de plusieurs facteurs : précision souhaitée sur les calculs, contours plus ou moins sinueux des limites, nombre et éloignement des singularités (puits) capacité de l’ordinateur Il est évident que pour une même précision il y aura un plus grand nombre de mailles carrées égales que de mailles inégales, choisies plus petites dans les zones sensibles. Mais un modèle à mailles égales est beaucoup plus facile à utiliser. Suivant le problème, la dimension des mailles peut être de l’ordre de 5 mètres (pour l’étude d’une digue ou d’un champ captant) à 5 kilomètres (pour l’étude d’une nappe régionale). En pratique, on essaiera de limiter le nombre de mailles à une valeur comprise entre 1000 et 10000 mailles (par couche). La géométrie de la nappe est définie (en chaque point) par : - Le substratum (parfois appelé « mur »), c’est à dire le fond de la ou des couche(s) aquifères, - Le toit, c’est à dire la limite supérieure de l’aquifère. Dans certains cas cependant si on est assuré que l’aquifère sera toujours en charge au cours des exploitations, ou bien si la variation de hauteur mouillée est faible, on peut se contenter de déterminer l’épaisseur de l’aquifère. - Les limites latérales peuvent être :  Des limites naturelles à POTENTIEL IMPOSE (c’est à dire à niveau imposé). En ces limites, le niveau est fixé (par un cours d’eau, par la mer, etc…) et le modèle calcule le débit d’échange qui permet d’imposer ce niveau, c’est à dire le débit que la rivière fournit à la nappe ou au contraire que la nappe fournit à la rivière ; le débit de la nappe qui se perd dans la mer, etc …  Des limites naturelles imperméables : collines, zones argileuses. En ces limites, l’aquifère se termine, mais un débit constant peut arriver (ou partir) par exemple par ruissellement sur un coteau. On appelle parfois ces limites à DEBIT IMPOSE, car le débit y est constant (ou nul).  Des limites arbitraires : il est indispensable d’étudier une nappe dans son ensemble, c’est à dire jusqu’à ses limites. Dans certains cas cependant, quand la nappe est très étendue, on pourra n’en étudier qu’une partie en prenant soin d’isoler une zone assez grande pour que les modifications ou les sollicitations (pompages) qu’on veut faire subir à la nappe aient des influences négligeables sur les limites introduites dans le modèle. On choisira donc une limite : - soit le long d’une ligne de courant : le débit d’échange par cette limite sera alors nul,

19

Modélisation hydrogéologique Master 2 HSE - soit le long d’une ligne équipotentielle sur laquelle on impose la charge mesurée sur le terrain. Le modèle donne alors le débit d’échange par cette limite. Si cette limite est suffisamment éloignée des exploitations prévues (captage), on peut supposer que ce débit ne changera pas. Pour simuler les exploitations prévues, on mettra alors une limite imperméable et on imposera le débit calculé précédemment. 3..3 Types de modèles Il y’a plusieurs façons de classer les modèles, ils peuvent être classés selon : *le régime permanent ou transitoire ; *le type de nappe (libre ou captive) ; *le nombre de dimension (domaine monodimensionnel, bidimensionnel ou tridimensionnel) ; *Model d’écoulement ou hydrodispersif (transport de soluté : variation de la concentration). 3.3.1. Le model bidimensionnel Quatre types de nappes peuvent être simulés : *Nappe captive : niveau piézométrique dépasse le toit ( T= K. m), avec m : épaisseur de la nappe captive : équation linéaire puisque m est constante * nappe libre : niveau piézométrique au -dessous du toit (T= K. n), avec n= épaisseur de la nappe libre : équation non linéaire *nappe avec drainance *nappe mixte : transformation d’une nappe libre en nappe captive. L’anisotropie peut être considérée selon les axes de direction x y.  Lors de la simulation de la nappe captive la transmissivité T et le coefficient d’emmagasinement S sont introduits dans chaque maille : T va dépendre soit de la variation de K soit de la variation de m Dans un écoulement bidimensionnel, on introduit T selon x et y donc introduire Tx et Ty.  La simulation de la nappe libre est basée sur la théorie de « Dupuit (1863) » : la composante verticale est négligée car elle est trop faible, donc les écoulements sont plus importants sur le plan horizontal que sur le plan vertical. Rappel Théorie de Dupuit Rabattement par puits en nappe libre Objectif : calculer le débit Q que l’on peut extraire d’un puits de rayon r atteignant le substratum imperméable, de sorte à maintenir une hauteur d’eau z o constante dans le puits, lorsque le régime permanent est atteint. L’influence du pompage se fait sentir sur une certaine distance R (rayon d'action) de l’axe du puits. Surface de la nappe = surface conique de révolution (cône de dépression). Isopièzes: cercles concentriques; les filets liquides convergent vers le puits.  nappe à filets convergents Débit Q à travers une surface cylindrique de rayon x et de hauteur z concentrique au puits: Q = q S et : S = 2 π x z

20

Modélisation hydrogéologique Master 2 HSE

Rabattement par un puits en nappe libre Q=K s

dH dz =K s dx dx

Q=2 πxz K x

dz dx R

h

dx =2 π K s∫ z .dz En intégrant : Q ∫ r x z

soit

Q(lnR−lnr )=

0

2 π Ks 2 2 (h −z 0) 2

h2−Z20 ¿ ¿ Q=π K s ¿ * signe positif car le flux q est opposé à la direction de l'axe x La critique de la théorie de DUPUIT a été faite par VIBERT et JAEGER , notons simplement que dans la mise en équation Dupuit à négligé la composante verticale de la vitesse de filtrations. Ceci n’est admissible que si la surface libre rabattue ne représente que de faibles pentes. Un fait plus gênant est que le cône de rabattement se raccorde au plan d’eau dans le puits, ce qui est nettement démenti par l’expérience.  Rabattement par puits en nappe libre Calcul de l’équation de la surface de la nappe z(x) Pour un point (x, z) quelconque, l'éq. de débit devient:

En égalant les 2 expressions fournissant le débit, on obtient: 21

Modélisation hydrogéologique Master 2 HSE

La surface de la nappe est convergente ou conique. z2 varie linéairement avec le logarithme de la distance x au puits. La surface de la nappe est donc représentée par une surface de révolution à écoulement convergent de type logarithmique.  Dans un système avec drainance l’effet de la drainance est représenté par le coefficient de la drainance A0, où K0 A 0= m0

Remarque : ne pas confondre ce coefficient avec le facteur de drainance B, B=

√

Km 0 K0

) 3.4 Conditions aux limites Pour pouvoir résoudre les équations précédentes, il est nécessaire de préciser les conditions aux limites qui peuvent être de différents types : a) Charge connue le long de la frontière (condition de Dirichlet) b) Flux (débit) connu le long de la frontière (condition de Neumann) c) combinaison de (a) et (b) (condition mixte ou de Fourier-Cauchy), le débit est calculé en fonction d'une charge.

22

Modélisation hydrogéologique Master 2 HSE

23

Modélisation hydrogéologique Master 2 HSE

24

Modélisation hydrogéologique Master 2 HSE

Résumé de conditions aux limites avec exemples  Charge imposée Les conditions aux limites en charge imposée Elles correspondent à la schématisation de conditions très usuelles, telles que : 

:

Un cours d’eau connecté à une nappe : la hauteur du cours d’eau se raccorde à une charge en nappe, et ce quelques soit le sens de l’écoulement (cours d’eau réalimentant la nappe ou nappe soutenant le débit du cours d’eau).

25

Modélisation hydrogéologique Master 2 HSE



Une source émergeant à une cote donnée. La charge en nappe se raccorde à la cote de la source



Une nappe subaffleurante : la charge est égale à la cote du sol.

Dans tous les cas, une condition à charge imposée doit être maniée avec beaucoup de précautions, car qui dit charge imposée dit absence de contrôle sur le flux. Dans le cas d’une modélisation, si les transmissivités dans la nappe ne sont pas représentatives de la réalité, les débits qui transitent au point de charge imposée seront tout à fait fantaisistes. Ainsi, il est dans le cas général plus judicieux de traiter une source ou une limite de recharge comme un lieu à débit imposé.  Flux imposé La condition à flux imposé la plus courante est la limite à flux nul, ou barrière hydraulique.

On emploiera également cette condition dans le cas d’une nappe rechargée par un taux d’infiltration connu ou bien pour représenter un débit soustrait à la nappe, par exemple au niveau d’un puits (on peut aussi traiter un puits en charge imposée ou en rabattement imposé, mais ce n’est guère prudent, ainsi qu’on l’a souligné précédemment.)

26

Modélisation hydrogéologique Master 2 HSE

 Conditions mixtes Les conditions aux limites en combinaison linéaire de la charge et du flux.

Elles correspondent par exemple au cas d’une rivière connectée à la nappe mais dont le fond est colmaté. Il y a donc perte de charge à la traversée du lit.

On peut alors écrire q=k*(hn-hr)/e=kn*grad(hn), d’où hn+kn/k*e*grad(hn)=-hr (avec hr= charge de la rivière, hn=charge en nappe, kn = perméabilité de la nappe, k= perméabilité du lit. 4.3.1 Introduction des conditions aux limites L’introduction des conditions aux limites dépend de leur localisation et du type de maillage adopté (point centré ou bloc centré). Les conditions aux limites sont directement imposées sur 27

Modélisation hydrogéologique Master 2 HSE le nœud pour le procédé des éléments finis et au centre des mailles pour le procédé des différences finis. Les conditions de potentiel sont imposées dans les mailles représentant les océans, la mer, les grands lacs et les rivières ayant une relation directe avec les grands aquifères. Les conditions de flux, représentant les flux sortant ou entrant à l’aquifère, sont simulés par des débits d’injection pour les débits entrants et des débits pompés pour les débits sortants, ces débits ont assigner au centre de la maille sur la face supérieure pour la méthode des différences finis. Pour la méthode des éléments finis les flux sont assignés sur la portion de la limite entre les nœuds, le programme les attribue ensuite aux nœuds.

Les conditions de flux nul sont imposées par : -les limites imperméables (faille), -la ligne de partage des eaux, -par les lignes de courants perpendiculaires à cette limite.

28

Modélisation hydrogéologique Master 2 HSE

La limite imperméable est simulée en donnant à T et K la valeur 0. Le choix d’imposer une condition de potentiel est du fait que c’est plus facile à mesurer un niveau piézométrique que de calculer un débit. 

Les conditions mixtes (Cauchy- Fourier) : Le débit traversant ce type de limite dépend de la différence de charge entre l’aquifère et la source qui sont séparés par une couche semi-perméable. Ce phénomène est appelé « drainance » et régit par l’équation suivante : K Qd = o ( h source−h ) . A m0

Avec : A : surface de la maille Qd : Débit de la drainance K0 : perméabilité de la couche semi-perméable m0 : épaisseur de la couche semi-perméable hsource : niveau d’une rivière (hr) ou niveau de l’aquifère adjacent h : niveau de l’aquifère qu’on simule Les conditions aux limites mixtes peuvent également représentées l’évapotranspiration, quand celle-ci dépend de la profondeur de la surface piézométrique et des arbres « phréatophytes » c.ad les racines utilisent les eaux de la nappe. Le débit évaporéest régit par la relation suivante : QET =Q ETM pour h>h s où : QET =Débit de l' évapotranspiration

29

Modélisation hydrogéologique Master 2 HSE QETM

' = Débit de l évapotranspirationmaximum

QET =Q ETM

h−(h s−d ) pour (h s−d )≤ h ≤ hs d

4. Modèles mono, multicouches et tridimensonnels Les modèles peuvent s’intéresser à une seule couche, à plusieurs ou même être tridimensionnels : - l’aquifère peut être considéré comme monocouche : les vitesses d’échange sont horizontales et la perméabilité en un point est considérée comme constante sur la verticale, - l’aquifère peut être multicouche, c’est à dire formé de plusieurs couches aquifères (sables, calcaire) séparées par des couches imperméables (argile), - l’aquifère peut être représenté par un modèle tridimensionnel si les vitesses verticales ne sont pas négligeables, par exemple au voisinage d’un puits ou une tranchée ou dans une galerie de mine. 5. La résolution En pratique les calculs se font par approximations successives. Partant d’un état initial des charges dans les mailles, on les recalcule toutes, les unes après les autres plusieurs fois, à 30

Modélisation hydrogéologique Master 2 HSE partir des équations précédentes avec les charges des mailles voisines et les conditions de débit, d’infiltration et de charge imposée aux limites. On arrête le calcul lorsqu’on obtient une quasi stabilité des charges recalculées. La méthode converge si on part de charges réalistes. On peut tenir compte au cours du calcul des variations de transmissivité avec la hauteur d’eau (eau de nappe libre), de l’assèchement ou du débordement en certains points. Les calculs en régime transitoire se font également par approximations successives avec un découpage en intervalles (subdivisions des pas de temps) en partant d’un état équilibré. Il est également possible de faire les calculs par inversion directe, mais cette méthode beaucoup plus lourde ne donne qu’une première solution qui doit être affinée par itérations pour tenir compte des éventuelles non-linéarités. Pour toutes les méthodes de résolution, on impose sur la périphérie du modèle étudié (qui forme les limites du modèle) : - soit les niveaux : le modèle calcule alors les débits qui transitent par cette limite, - soit une perméabilité nulle : c’est une limite imperméable, - soit (plus rarement) un débit d’échange imposé (entrant ou sortant).

6. Les données nécessaires et le calage du modèle Il faut connaître débit d’échange et d’alimentation, ainsi que perméabilité et coefficient d’emmagasinement : - Débit d’échange et d’alimentation. Il faut déterminer de façon aussi précise que possible les débits pompés et celui des sources ainsi que la recharge par la pluie corrigée de l’évapotranspiration et du ruissellement. Pour régler le modèle on cherche à déterminer les transmissivités produisant les charges observées quand la nappe est exploitée avec les débits constatés. - La perméabilité n’est connue qu’en quelques points par interprétation des pompages d’essais d’où on déduit, en fonction de la géologie une carte approchée des perméabilités. Le calage en corrigera les valeurs. - Coefficient d’emmagasinement. Sa connaissance n’est nécessaire que si on s’intéresse au régime transitoire. Il est mal connu sauf en certains points grâce aux pompages d’essais. On l’ajuste lors du calage notamment si on dispose de plusieurs champs piézométriques à différentes dates. En général, on connaît assez bien la géométrie de l’aquifère et mal la perméabilité, le colmatage des cours d’eau et la recharge par infiltration. Le calage consiste à ajuster ces valeurs pour reproduire au mieux : - Niveaux en tous points ; - Débits des sources et émergences. Si on dispose de plusieurs états de la nappe, à des périodes différentes, (basses et hautes eaux par exemple) on peut régler les perméabilités sur un état et les contrôler sur le second. Pour régler l’emmagasinement il faut disposer d’une variation de niveau au cours du temps, en quelques points, ou sous l’effet de variation des débits d’exploitation par exemple, ou encore de deux cartes piézométriques dont l’une pour un équilibre non atteint. Le calage du modèle permet d’ajuster les paramètres et également de vérifier la cohérence des données et du schéma d’interprétation choisi. Un premier modèle grossier permet de déterminer les zones où il faut compléter les mesures.

6.1 Singularités Si les charges varient lentement dans la maille ou bien si elles varient linéairement, la charge calculée correspond à peu près au centre de la maille. Au voisinage d’une singularité (un captage par exemple), la charge calculée ne correspond plus du tout à celle d’un puits au

31

Modélisation hydrogéologique Master 2 HSE centre de la maille, mais il faut appliquer une correction donnée par la formule suivante (valable pour des mailles carrées) :

avec a = Côté de la maille rp = Rayon du puits T = transmissivité de la maille

6.2 Aquifères complexes : modélisation quasi-tridimensionnelle Pour la modélisation de systèmes complexes, on utilise un modèle multicouche, ou un modèle quasi tridimensionnel. Un modèle quasi-tridimensionnel prend les hypothèses suivantes (figure Schéma d’un modèle multicouche ou quasi-tridimensionnel (d’après Hydroexpert)

Le système multicouche est représenté par un empilement d’aquifères séparés par des semiperméables ; - Les écoulements sont strictement bidimensionnels horizontaux dans les aquifères et monodimensionnels verticaux dans les semi-perméables. Dans le premier cas le tenseur de perméabilité se réduit à sa seule composante horizontale (Kx = Ky = Kh), dans le second à sa seule composante verticale (Kz = Ky) ; - Toute différence de charge entre deux aquifères séparés par une semi-perméable produit un débit d’échange par drainance ; Les aquifères sont géométriquement délimités par les cotes de leur toit et de leur mur ; - Les semi-perméables peuvent être réels (i.e. la molasse du Gatinais, les marnes vertes…) ou fictifs dès lors que l’on veut différencier deux réservoirs de comportement hydraulique distinct ; - L’extension des semi-perméables est automatiquement définie par le toit de l’aquifère sousjacent et le mur de l’aquifère sus-jacent ; - Les semi-perméables sont supposés non capacitifs ; le volume d’eau qu’ils peuvent libérer est négligeable par rapport au volume d’eau disponible dans les aquifères (ce qui est effectivement vérifié dans le cas présent compte tenu de la faible épaisseur des semiperméables) ; - Le rôle hydrodynamique des semi-perméables (échanges entre aquifères) fictifs ou réels dont l’utilité ou la présence n’est que locale mais la représentation généralisée à l’échelle du modèle (pour les raisons évoquées ci-dessus), est numériquement neutralisée partout où sa présence n’est pas nécessaire. La neutralisation s’opère en agissant sur la perméabilité verticale KV. Une difficulté importante est la détermination des débits drainés par les rivières du réseau hydrographique interne (y compris la localisation et l’évolution des assecs). La représentation qui en est faite par la modélisation est représentée en figure

32

Modélisation hydrogéologique Master 2 HSE

Représentation Hydroexpert)

des

émergences

par

modèle

quasi-tridimensionnel

(d’après

En situation réelle le mécanisme de débordement vers des cours d’eau superficiels écrêtant juste la surface libre de la nappe se traduit obligatoirement par l’apparition d’une singularité où l’écoulement, tout en restant majoritairement bi- dimensionnel horizontal, présente une composante verticale non négligeable. La prise en compte de cette structure locale particulière de l’écoulement (singularité) est nécessaire pour calculer correctement le débit. La solution consiste à introduire un horizon semi-perméable fictif, de faible épaisseur (1 m) au regard de l’épaisseur totale de l’aquifère, qui sépare celui-ci en deux parties : - une partie supérieure à la surface libre et à transmissivité linéairement variable en fonction de la hauteur piézométrique ; - une partie inférieure captive, à transmissivité fixe, sous le toit du semi-perméable. La condition du débordement est alors imposée sur l’aquifère supérieur uniquement. Notons que dans ce cas il convient d’affiner la discrétisation au niveau de la rivière afin d’accroître la précision des calculs. Ainsi en jouant à la fois sur la perméabilité verticale du semi-perméable et sur la perméabilité horizontale de « l’aquifère supérieur » il est possible d’ajuster très précisément le débit drainé en respectant l’effet tridimensionnel local. L’imposition d’une condition de débordement sur une maille ayant l’épaisseur totale de l’aquifère constituerait une grave erreur, dans la mesure où elle ne rendrait pas compte de 33

Modélisation hydrogéologique Master 2 HSE l’effet tridimensionnel et donc obligerait, pour reproduire les débits, à diminuer exagérément la perméabilité horizontale de l’aquifère avec pour conséquence un impact régional trop important sur la piézométrie.

34