Solutions analytiques du problème

On étudie ici le problème d'intrusion saline sous un îlot. Pour simplifier le problème on prendra le cas d'une bande de terre rectiligne afin d'établir une symérie plane du problème. La bande de terre sera soumise à une recharge par le dessus (symbolisant la pluie). Dans un premier temps, on va détailler les équations à résoudre et les hypothèses faites pour simplifier les calculs. Ensuite on représentra les solutions analytiques grâce à un programme Matlab écrit par Mr Rachid Ababou (professeur à l'ENSEEIHT et chercheur à l'Institut de Mécanique des Fluides de Toulouse (IMFT), groupe de recherche Milieu Poreux).

 

Schéma du problème

 

Représentation schématique du problème et des axes

 

Présentation des paramètres et des inconnues

 

  • ρsalé : masse volumique de l'eau salée (ρsalé= 1025 kg/m3)
  • ρdouce : masse volumique de l'eau douce (ρdouce=1000 kg/m3)
  • ε : contraste de densité entre eau douce et eau salée ( ε=(ρdsalé - ρdouce)/ρdouce d'ou ε=0,025 )
  • I : taux de recharge nette de la nappe d'eau douce (on a choisi ici de prendre une moyenne annuelle des précipitations soit I=700 mm/an c'est à dire I=8.7 x 10-8 m/s)
  • Lo : distance entre la mer et le centre de l'ile (Lo=1000 m)
  • L : longueur de pénétration du coin salé sous l'eau.
  • Q(X) : densité de flux de l’eau douce vers la mer (verticalement intégrée) en [m²/s] ou débit spécifique en [m3/s/mètre de côte].
  • T : transmissivité hydraulique de la nappe d'eau douce (m2/s)
  • K : conductivité hydraulique du massif (pour évaluer K, nous avons choisi d'utiliser la formule de Kozeny-Karman (2), avec une porosité estimée à 25 % (sols essentiellement sableux), une constante de Kozeny égale à 5 (empilages de grains isométriques avec des porosités n'excédant pas 0,8) et un diamètre moyen du milieu dmoyen=35 micromètres (3), on trouve la perméablilité k= 5*10-9 m²
    et que pour l'eau à 15.5°C, K=0.873 x 10e7 x k
    d'où K = 4.356 x 10-3 m/s.
  • Ho : profondeur du plancher imperméable s'il existe (m). Puisque l'on ne connait pas la profondeur ni l'existence du plancher on imaginera deux cas, le cas où il est profond et le cas où il ne l'est pas.
  • H(X,Z) : profondeur de l'interface par rapport à la mer (m)
  • h(X,Z) : hauteur de la surface d'eau douce par rapport à la mer (m)


Hypothèses de travail

 

  • interface abrupte entre eau douce et eau de mer, c'est à dire que l'on néglige la diffusion de sel.
  • les lignes de courants suivent le front salé et ne peuvent pas le traverser sauf près de la surface de la mer (point 0)
  • on prend le cas d'un aquifère homogène
  • substratum imperméable plan et horizontal à une profondeur H0.
  • vitesses quasi horizontales dans l'eau douce, on utilisera les équations verticalement intégrées d'écoulement plan (approximation de Dupuit-Boussinesq) d'où les pressions verticalement hydrostatiques.
  • on suppose que l'eau de mer (le coin salé) est totalement hydrostatique.
  • la nappe d'eau douce est à surface libre (pas de supersubstratum)

Présentation des équations à résoudre

 

  • Bilan de masse

On fait un bilan de masse sur la moitié droite de l’ilôt, on obtient par bilan de masse (égalité des flux entrant et sortant en régime permanent) avec Q(x) la composante suivant X du "débit spécifique".

 

  • Loi de Darcy (intégré verticalement)

On écrit la loi de Darcy sur le demi-domaine là où se trouve l'île.

 

On a maintenant deux équations (Darcy, bilan de masse) mais trois inconnues ( h(x,z), H(x,z) et L), on a donc besoin d'une troisième équation.

  • Relation de Ghyben-Herzberg sur l’interface salée

On rappelle que le niveau de la mer est en z=0 et que l'on a fait les approximations de Ghyben-Herberz (eau douce verticalement hydrostatique, eau salée hydrostatique) et la continuité de pression à l'interface. On note que Zd la cote de la surface libre d'eau douce et Zs la cote de l'interface salée. En suivant un circuit hydraulique en forme de « tube en U », on obtient :


Pression en un point B de la surface libre (coté eau douce) : Pb = Patm + ρdouce x g x (H + h)
Pression au même point B de la surface libre (coté eau salée) : Pb' = Patm + ρsalé x g x H

Par continuité de pression Pb=Pb', on obtient :

 

Résolution du problème et solutions analytiques

  • Solution dans la zone du coin salé, c'est à dire pour X < L on a que h= ε x H , d'ou

Il reste maintenant à trouver les constantes A et B, on fait pour cela une hypothèse simplificatrice supplémentaire: on admet que la totalité du débit d'eau douce Qo passe par le point triple interface/côte/mer (point 0). On est en régime permanent donc tout ce qui rentre par infiltration doit aussi sortir vers la mer, on obtient le débit Qo qui sort vers la mer : Q(0)= I x Lo.

or Q(X)= IX + A d'où Q(0) = A = -Qo (débit sortant) c'est à dire A = - I x Lo et ainsi Q(X) = I (X - Lo)

De plus, en X=0, la surface libre de la nappe d’eau douce atteint le niveau de la mer, soit h(0)=0 ce qui nous donne B=0. Les constantes déterminées, il ne reste plus qu'a réecrire les solutions et on trouve :

  • Calcul de la longueur du coin salé

Si le calcul montre que la profondeur H(X) est supérieure à Ho pour X < Lo, alors l’interface salée est interceptée par le substratum imperméable à une certaine distance de la côte (on a noté cette longueur L). Plusieurs cas peuvent se présenter et L peut exister ou pas si Ho est trop grand et que dans ce cas l'eau salée pénètre complétement sous l'île. Pour le savoir, il faut résoudre l'équation H(L) = Ho. En insérant l’expression précédente de H(x) dans cette relation, on obtient une équation du 2nd degré en L, soit

I x L2+ (2 x Qo) x L − (K x ε x (ε+1) x Ho)=0

  • Solution hors du coin salé, c'est à dire pour L < X < Lo

Tout d'abord, on sait que H(X)=Ho. Ensuite il reste à trouver l'équation de h(X). Pour cela, il faut réutiliser la solution de h(X) donnée avant la détermination des constantes A et B et trouver leurs nouvelles valeurs. La première condition limite qui nous donnait A est toujours valable d'où A= -I x Lo. On sait que H(L)=Ho en utilisant l'équation de Ghyben-Herzberg, on obtient que h(L)= ε x Lo ce qui est la seconde condition limite qui nous permet de trouver la constante B. On pourrait vérifier que flux Q(X) est nul au centre de l'île Q(Lo)=0 par symétrie du problème. On trouve alors que

 

Représentation avec Matlab des solutions analytiques

Les valeurs paramètres choisies pour les simulations analytiques sont celles données et détaillées au début de cette page. Les équations des solutions utilisées dans Matlab dans le cas d'une bande de terre sont données dans le paragraphe ci-dessus.

Aussi, puisque l'on ne connaît pas la profondeur du substrat imperméable, on fera une simulation où l'interface eau salée, eau douce intercepte le substrat et une deuxième où elle ne l'intercepte pas. Dans un troisième temps, on a représenté les solutions analytiques dans le cas d'un îlot rond (les équations des solutions ne sont pas détaillées sur ce site internet).

  • Substrat imperméable profond Ho=30 mètres


Solution analytique du problème d'intrusion saline représenté avec Matlab en coupe et en 3D

Sur la figure précédente, on peut observer la lentille d'eau douce qui se forme sur le dessus de la nappe, avec un hmax d'environ 1 mètre dans ce cas.

  • Substrat imperméable peu profond Ho=20 mètres

Solution analytique du problème d'intrusion saline représenté avec Matlab en coupe et en 3D
 
  • Cas d'une île circulaire



Dans le cas d'une île circulaire, on a aussi représenté le cas où le substrat est profond (Ho= 30m) et celui où il ne l'est pas  (Ho=15m)

Solution analytique du problème d'intrusion saline représenté avec Matlab en coupe et en 3D (Ho = 30m)

Solution analytique du problème d'intrusion saline représenté avec Matlab en coupe et en 3D (Ho = 15m)


Conclusion
 

On remarque que la résolution de ce problème simple est déja un peu complexe. La résolution analytique de ce problème en rajoutant un pompage existe aussi mais suppose d'avoir un puits profond jusqu'au plancher imperméable et situé au milieu de l'île (ce qui n'est absolument pas le cas de notre situation) afin de pouvoir conserver les hypothèses de vitesses quasi horizontales dans l'aquifère et de résoudre ainsi facilement le problème. De plus, tous les calculs ont été faits en régime permanent ce qui est très limitant pour calculer un temps caractéristique de pollution du puits par l'eau salée.

Pour toutes ces raisons, nous avons choisi de faire une simulation numérique du problème avec le logiciel Comsol.

 


Bibliographie
1. L'ensemble des solutions analytiques sont tirés du polycopié suivant (Rachid Ababou, Elements d'hydrologie souterraine, 2007)
HSOUT : cours de 3ème année option SEE/MFN/Master H2SE

2. Formule de Kozeny-Karman (J.C.Charpentier, Techniques de l'ingénieur, 2003)

3. Rapport de l'institut des facultés des sciences et techniques de l'université Cheikh Anta Diop (Jean Henri Sene, Evaluation du potentiel agro-écologique de l'estuaire du Siné-Saloum, 2005)