0% ont trouvé ce document utile (0 vote)

148 vues132 pages

Cours GMRES

Transféré par

Nous prenons très au sérieux les droits relatifs au contenu. Si vous pensez qu’il s’agit de votre contenu, signalez une atteinte au droit d’auteur ici.

Formats disponibles

Téléchargez aux formats PDF, TXT ou lisez en ligne sur Scribd

0% ont trouvé ce document utile (0 vote)

148 vues132 pages

Cours GMRES

Transféré par

manalmoussaoui41

Nous prenons très au sérieux les droits relatifs au contenu. Si vous pensez qu’il s’agit de votre contenu, signalez une atteinte au droit d’auteur ici.

Formats disponibles

Téléchargez aux formats PDF, TXT ou lisez en ligne sur Scribd

Institut Agronomique et Vétérinaire Hassan II

Département du Génie Rural

Cours
Modélisation et Gestion
des Ressources en
Eau Souterraine
Ali HAMMANI

3ème année Génie Rural

Décembre 2021
Table des matiÃ¨res

1 Chapitre 2: Modélisation numérique des écoulements souterrains. . . . . . . . . . . . . . . . . . . . 1

1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
Définitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
Différents types de modèles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Choix d'une technique numérique et d'un code . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Méthode de différences finie en aquifère homogène et isotrope . . . . . . . . . . . . . . . 4
Ecoulement en régime permanent . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
Ecoulement en régime transitoire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.4 Méthodes de différences finies en aquifère hétérogène et anisotrope . . . . . . . . . . . 11
Equation l'écoulement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
Approximation de l'équation d'écoulement par les différences finies . . . . . . . . . . . . . 12
Terme source W(x,y,t) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.5 Validité des solutions numériques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.6 Calage des modèles numériques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.7 Résolution du problème inverse de l'écoulement souterrain . . . . . . . . . . . . . . . . . . . 16
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
Paramétrisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
Identification des paramètres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.8 Analyse de sensibilité d'un modèle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2 Chapitre 2: Modélisation numérique des écoulements souterrains. . . . . . . . . . . . . . . . . . . . 18

2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
Définitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
Différents types de modèles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2 Choix d'une technique numérique et d'un code . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.3 Méthode de différences finie en aquifère homogène et isotrope . . . . . . . . . . . . . . . 21
Ecoulement en régime permanent . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
Ecoulement en régime transitoire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.4 Méthodes de différences finies en aquifère hétérogène et anisotrope . . . . . . . . . . . 28
Equation l'écoulement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
Approximation de l'équation d'écoulement par les différences finies . . . . . . . . . . . . . 29
Terme source W(x,y,t) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.5 Validité des solutions numériques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.6 Calage des modèles numériques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.7 Résolution du problème inverse de l'écoulement souterrain . . . . . . . . . . . . . . . . . . . 33
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
Paramétrisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
Identification des paramètres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.8 Analyse de sensibilité d'un modèle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

0
Etudedesnappessouterraines
Chapitre 1

parlebilanhydrogéologique
Dans la gestion des ressources en eau souterraines, l'homme intervient dans le cycle hydrologique
dans le but de tirer profit. Cette intervention prend la forme de modifications imposées sur les
variations des composantes du bilan hydrique des nappes.
Une des parties la plus importante de l'hydrogéologie est, donc, la connaissance du bilan hy-
drogéologique. Déterminer la part de l'eau qui circule souterrainement est un problème majeur,
puisque c'est cette part qui fournit les ressources utilisables. C'est un problème difficile à résoudre.
Le bilan se réalise dans le but de déterminer la contribution de chaque terme dans l'alimentation et le
rabattement de la nappe et de connaître et évaluer sa réserve dans le temps (notamment pendant
la période de sécheresse ou les périodes à fortes précipitations). Le bilan peut servir comme outil
de prise de décisions quant à la gestion de la ressource en eau souterraine.
Le bilan peut être calculé pour toute une région, par exemple pour un bassin fluvial. C'est facile.
Mais il englobe toutes les nappes du bassin, les terrains perméables et imperméables, les zones
saturées et non saturées.

Les différents éléments du bilan

Examinons le bilan d'un horizon aquifère pour un intervalle du temps (figure ??). Les entrées dans le
système aquifère peuvent se composer des éléments suivants:
1. Entrée d'eau d'autres horizons aquifères voisins ou sous-jacents;
2. Entrée d'eau directement des eaux de surface composée de précipitation P et/ou d'apport
artificiel par irrigation I;
3. Entrée d'eau des réserves d'eaux souterraines (abaissement du niveau de la nappe −∆S).
Les sorties d'eau de l'horizon aquifère peuvent se composer des éléments suivants:
1. Départ d'eau vers d'autres horizons aquifères −L;
2. Départ d'eau par évaporation E et par évapotranspiration des cultures ET ;
3. Départ d'eau par des sources à l'air libre; par drainage naturel Dn (vers les oueds et/ou
d'autres cours d'eau naturels) ou par drainage artificiel vers le réseau de drainage.
4. Départ d'eau par exploitation de la nappe (pompages...).
5. Quantité d'eau passée en réserve des eaux souterraines+∆S (une remontée de la nappe
phréatique).

1
Chapitre 1 : Etude des nappes souterraines par le bilan hydrogéologique 2

Equation générale du bilan

Selon que le système défini, l'équation du bilan s'écrira différemment. Deux cas de schéma sont
envisagés. Dans le premier cas le système se compose de tout le bassin. Dans le deuxième le système
se compose seulement de la couche aquifère; on s'intéresse donc à ce qui entre et sort du système
aquifère. C'est le bilan de ce dernier qu'il faut calculer.

Equation du bilan hydrologique globale d’un bassin versant

AP + AI + L = I + E + R + Ex + ±∆S (1.1)

Où :
• Ap et Ai sont les apports respectivement par les précipitations et par irrigation;
• L : Apports occultes;
• I : Percolation profonde des eaux de pluie et d'irrigation;
• E : Composée de l'évapotranspiration des cultures et de l'évapotranspiration directe à partir
de la nappe dans le cas où la nappe est très proche de la surface du sol;
• R : Terme qui englobe le ruissellement de surface le drainage à partir de la nappe.
• Ex : Terme qui englobe les quantités d'eau prélevées de la nappe pour l'exploitation dans le
cas où cette eau est utilisée à l'extérieur du bassin.
• ∆S est la variation du stock d'eau dans le sol. Pour une période longue de l'ordre de l'année
∆S peut être négligée.

Equation du bilan hydrogéologique dans une nappe phréatique

La nappe phréatique est supposée comme étant un réservoir inerte; dans un tel cas l'équation du
bilan s'écrit:
IP + II + L = E + D + Ex + ±∆S (1.2)

Où :

Ali HAMMANI, 2021

Chapitre 1 : Etude des nappes souterraines par le bilan hydrogéologique 3

• Ip et Ii sont respectivement les infiltrations efficaces des eaux de précipitation et d'irrigation,

c'est-à-dire les volumes d'eau qui atteindront réellement à la nappe;
• E est l'évapotranspiration directe à partir de la nappe;
• D est le volume d'eau drainée à partir de la nappe;
• Ex étant les volumes prélevés à partir de la nappe à des fins agricoles, industriels ou domes-
tiques.
• ∆S est la variation de la réserve de la nappe.

Evaluation des différents termes du bi-

lan
Les entrées
L'alimentation d'une nappe phréatique peut se réaliser par l'infiltration des eaux de précipitations,
la percolation des eaux d'irrigation et les apports naturels constitués par les apports latéraux et/ou
verticaux par drainance.

7.3.1.1. Evaluation des apports par les précipitations

a. Calcul de la lame d’eau tombée sur le périmètre
La lame d'eau moyenne tombée sur un bassin est égale au quotient de la précipitation en millimètre
par la superficie en m2 . Elle peut ôtre calculée par plusieurs méthodes dont les principales sont les
suivantes:
• La moyenne arithmétique des précipitations;
• L'utilisation des isohyètes;
• La méthode de Theissen.
i) Calcul de la lame d’eau par la moyenne arithmétique
La méthode la plus simple consiste à calculer la moyenne arithmétique des hauteurs de précipitation
relevées dans le même intervalle de temps aux diverses stations du bassin.
ii) Calcul de la lame d’eau par l’utilisation des isohyètes
C'est le procédé le plus rationnel et le plus précis (figure ??). Le calcul de la lame d'eau se fait de
la façon suivante :
• Cartographie des courbes isohyètes (Courbes de même hauteur pluviométrique).
• Planimétrage de la surface partielle Si de la nappe comprise entre deux isohyètes successives
Pi et Pi+1 .
• Calcul de la pluviométrie moyenne par la formule suivante:

∑ Pi + Pi+1 Si
P = Pi (1.3)
i
2 S

iii) Calcul de la lame d’eau par la méthode de Theissen

La méthode de calcul avec les cartes isohyètes est longue. C'est pourquoi elle est remplacée souvent
par celle de Theissen qui est plus rapide (figure ??). Les stations pluviométriques sont reportées sur
une carte. On relie les stations adjacentes par des droites au milieu de chacune desquelles est
élevée une perpendiculaire. Les intersections de ces perpendiculaires déterminent des polygones.
Dans chaque polygone, la hauteur de précipitation choisie est celle relevée à la station située à
l'intérieur. La surface de chaque polygone élémentaire est calculée en mètre carré est exprimée

Ali HAMMANI, 2021

Chapitre 1 : Etude des nappes souterraines par le bilan hydrogéologique 4

en pour-cent de la surface totale du bassin. Ce pourcentage est utilisé comme coefficient de

pondération propre à chaque station. Donc, si P est la hauteur de précipitation recherchée sur
tout l'étendue de la nappe de superficie A, Pi celle relevée sur le polygone de surface Ai , nous
obtenons:

∑ Ai
P = Pi (1.4)
A

b. Calcul de la lame d’eau efficace

La lame d'eau efficace est définie comme étant la proportion d'eau de pluie qui alimente réellement
la nappe. Pour l'évaluer, on utilise la formule suivante:

Ip = P iS (1.5)

Où :
• Ip est l'infiltration efficace des eaux de pluie; P est la hauteur de la pluviométrie moyenne
calculée sur tout le domaine de la nappe.
• S est la superficie totale du domaine de la nappe.
• i est le coefficient d'infiltration. Il peut ôter estimé par plusieurs méthodes.
c. Détermination du coefficient d’infiltration

Ali HAMMANI, 2021

Chapitre 1 : Etude des nappes souterraines par le bilan hydrogéologique 5

L'eau de précipitation et/ou d'irrigation, après avoir subi des pertes par évapotranspiration et
absorption par la couverture végétale parvient à la surface du sol où elle se répartit en deux
fractions : Le ruissellement de surface et l'infiltration.
L'infiltration représente la quantité d'eau qui pénètre dans le sol et le sous-sol où elle alimente les eaux
souterraines: eau de rétention, écoulement hypodermique, écoulement souterrain et reconstitution
des réserves aquifères.
L'évaporation dépense une part importante du stock d'humidité. Ainsi une fraction seulement des
précipitations alimente les nappes. C'est l’infiltration efficace.
Le taux d'infiltration est exprimé en millimètres de hauteur d'eau. Les hydrologues utilisent très souvent
le terme de coefficient d'infiltration qui est le quotient, exprimé en pourcentage, du taux d'infiltration
par la hauteur de précipitation.
L'infiltration est, en hydrologie, le facteur le plus important du cycle de l'eau et aussi le plus difficile
à évaluer car il échappe aux mesures directes par des processus simples.
Les méthodes de mesure de l'infiltration sont nombreuses et variées en relation avec les phénomènes
étudiés. Elles peuvent ôter directes ou indirectes.
i) Méthodes directes
Elles permettent d'évaluer la quantité d'eau infiltrée sur une surface de sol déterminée. Elles com-
prennent:
• La méthode expérimentale par alimentation en eau à l'aide de fossés, rigoles, tranchées ou
tubes crépinés enfoncés dans le terrain;
• Le débit des sources dans un bassin hydrogéologique individualisé en fonction de la plu-
viométrie;
• Les lysimètres;
• Les mesures de gradient vertical d'humidité dans le sol;
• L'étude des variations de la surface piézométrique dans un bassin hydrogéologique et des
oscillations du niveau piézométrique dans les puits. Cette méthode sera développée par la
suite.
ii) Méthodes indirectes
Dans les méthodes indirectes l'infiltration est calculée par différence des entrées et des sorties. Les
autres termes du bilan hydrologique, précipitation, ruissellement et évapotranspiration étant connus.
Pour une longue période (de l'ordre de l'année par exemple), on peut négliger la variation de la
réserve et donc le bilan peut s'écrire :

P =E+R+I (1.6)

Où :
• P : Précipitation;
• E : Evapotranspiration;
• R : Ruissellement;
• I : Infiltration.
On en déduit l'infiltration par différence :

I =P −E−R (1.7)

Et le coefficient d'infiltration i :
I
i= (1.8)
P

Ali HAMMANI, 2021

Chapitre 1 : Etude des nappes souterraines par le bilan hydrogéologique 6

7.3.1.2. Evaluation des apports par irrigation

a. Calcul du coefficient d’infiltration
Contrairement aux apports par la pluviométrie, l'infiltration des eaux d'irrigation peut se faire à deux
niveaux :
• Au niveau des canaux d'irrigation d'ordre supérieur; dans ce cas une partie du volume d'eau
qui transite par le réseau, alimente la nappe d'une façon systématique. Par exemple, les pertes
au niveau du réseau sont estimées à 10% du volume d'eau distribué dans le périmètre de
Béni-Moussa du Tadla).
• Au niveau de la parcelle, la percolation des eaux d'irrigation varie avec la variation de
l'évapotranspiration des cultures. A ce niveau le coefficient d'infiltration par la méthode du
bilan annuel.
b. Formule de calcul des apports par irrigation
Comme dans le cas des apports météoriques, le calcul des apports nets à la nappe s'effectue
moyennant la formule suivante:
AI = V (i + i′ ) (1.9)

Où:
• V = Volume d'eau provenant du barrage et distribué par le réseau d'irrigation;
• i = Coefficient d'infiltration au niveau de la parcelle;
• i′ = Coefficient d'infiltration au niveau du réseau d'irrigation.

7.3.1.3. Evaluation des apports latéraux

a. Source des apports latéraux
On désigne par apports latéraux tout flux d'eau horizontal, dirigé vers la nappe et permettant son
alimentation naturelle. Pour les caractériser, on se base sur les cartes des courbes piézométriques qui
peuvent indiquer le sens d'écoulement, et on peut donc localiser les zones où les apports latéraux
ont lieu (figure ??).

b. Méthode d’approche
L'estimation des flux latéraux peut ôter faite en utilisant la loi de Darcy.
Pour ce faire, on doit avoir les données piézométriques, et hydrodynamiques pour un ensemble de
points situés à l'extrémité où les apports se réalisent. Le calcul se fait moyennant la formule suivante:
∆H
Q(m3 /s) = T L (1.10)
l

Ali HAMMANI, 2021

Chapitre 1 : Etude des nappes souterraines par le bilan hydrogéologique 7

Dans cette formule :

• T est la transmissivité.
• L est la longueur du front de la nappe, c'est la longueur de la section à travers laquelle se
fait l'écoulement. Cette longueur peut ôtre mesurée directement sur les cartes piézométriques.
• ∆Hoverl est le gradient hydraulique moyen.

7.3.1.4 Les apports par drainance

La drainance qv (volume d'eau par unité de surface et par unité de temps) à travers une couche
semi-perméable (vers le haut ou vers le bas) d'une nappe de piézométrie hext à une autre nappe
de hauteur piézométrique h, est donné par.
hext − h
Qv = K ′ (1.11)
B′

Où K ′ et B ′ sont, respectivement, la conductivité hydraulique et l'épaisseur de la couche semi-

perméable.

Les sorties:
Les prélèvements à partir de la nappe peuvent ôter constitués des pompages, du drainage naturel
et artificiel, de la drainance et de l'évapotranspiration directe à partir de la nappe.

7.3.2.1. Evaluation des prélèvements par pompages et par drainage artificiel

Le drainage se réalise, généralement, pour rabattre le niveau de la nappe. Le réseau de drainage
peut ôter soit des fossés ou des drains enterrés. Les volumes drainés sont connus si des stations
de mesures sont installées au niveau des principaux drains et collecteurs de drainage. Dans le cas
contraire ils peuvent ôter calculé moyennant un bilan hydrique à l'échelle du bassin versant.
Quant aux pompages, il est effectué, généralement, dans le but de l'exploitation de la nappe pour
des fins agricoles, industriels ou domestiques. Les volumes pompés sont connus si les débits des
stations de pompage sont connus. Sinon, ils peuvent ôter estimés par les besoins en eau des cultures
et par les besoins en eau domestiques et industrielles.

7.3.2.2. Le drainage naturel et la drainance vers le bas

Le drainage naturel peut ôter considéré comme étant un flux latéral sortant de la nappe. Son calcul
se fait par la formule de Darcy appliquée près de la barrière jouant le rôle de drain. De momie la
drainance, peut ôter calculée en utilisant l'équation (??).

7.3.2.3. Evaporation et évapotranspiration

Les nappes phréatiques sont soumises à l'action de l'évapotranspiration qui épuise les réserves en
eau souterraines. Mais cette action est limitée à une profondeur relativement faible et à la zone
d'évapotranspiration (figure ??). Elle est, toutefois, importante dans les régions semi-arides et arides
où elle augmente la teneur en sel des eaux et des sols.
White, 1932 a proposé la variation de l'évapotranspiration à partir d'une nappe, en fonction de la
profondeur de la surface de l'eau, et exprimée en pourcentage de l'évaporation bac.
D'une manière générale, l'évapotranspiration à partir de la nappe cesse dès que la frange capillaire
atteint une profondeur donnée. Dans les régions arides, il est facile de reconnaître cette profondeur.
C'est la profondeur ∆o dite critique au-delà de laquelle il ne se forme plus de dépôts salins.
D'après Kovda cette profondeur serait en climat continental:

∆ = 170 + 8T ± 15cm (1.12)

Ali HAMMANI, 2021

Chapitre 1 : Etude des nappes souterraines par le bilan hydrogéologique 8

T étant la température moyenne annuelle en o C.

Variation de la réserve
La différence entre les entrées et les sorties pour une période où le bilan a été établi donne la
variation de la réserve de la nappe. Pour une nappe de surface A, un volume d'eau ∆S stocké
dans un aquifère cause une remontée du niveau de la nappe de ∆h:

∆S = SA∆h (1.13)

Où : S est le coefficient d'emmagasinement de l'aquifère. Si la remontée n'est pas uniforme et

l'emmagasinement varie d'une zone à l'autre, on peut subdiviser l'aquifère en N parties telles que:

∑
N
∆S = Sj Aj ∆hj (1.14)
j=1

Exemple de bilan hydrogéologique :

Nappe de Berrechid
[2]: #importer les modules requis
import os, json
import geopandas as gpd
from ipyleaflet import Map, GeoData, basemaps, LayersControl, ScaleControl, \
FullScreenControl, basemap_to_tiles, Choropleth, MeasureControl

# objet map
m = Map(center=(33.25, -7.57), zoom = 10, basemap = basemaps.Esri.WorldImagery)
m.layout.width = '100%'
m.layout.height = '700px'

m.add_layer(basemaps.Esri.WorldStreetMap)
m.add_control(ScaleControl(position='bottomleft'))
m.add_control(FullScreenControl(position='topright'))
m.add_control(LayersControl())
m.add_control(MeasureControl(position='bottomleft',active_color =␣
,→'orange',primary_length_unit = 'kilometers'))

Ali HAMMANI, 2021

Chapitre 1 : Etude des nappes souterraines par le bilan hydrogéologique 9

display(m)

selDf = gpd.read_file('data/NappeBerchidXY.shp')
geoDfA = selDf.to_crs(4326)
geoDataAquifer = GeoData(geo_dataframe=geoDfA, name='NappeBerchidXY')
m.add_layer(geoDataAquifer)

geoDfB = gpd.read_file('data/BassinsAmont.shp')
geoDataBasins = GeoData(geo_dataframe=geoDfB, name='Bassins Amont')
m.add_layer(geoDataBasins)

selDfP1016 = gpd.read_file('data/Piezo1016a.shp')
geoDfP1016 = selDfP1016.to_crs(4326)
geoDataPiezo1016 = GeoData(geo_dataframe=geoDfP1016, name='Piézo 10/2016')
m.add_layer(geoDataPiezo1016)

/Users/alihammani/opt/miniconda3/lib/python3.9/site-
packages/geopandas/_compat.py:111: UserWarning: The Shapely GEOS version
(3.9.1-CAPI-1.14.2) is incompatible with the GEOS version PyGEOS was compiled
with (3.10.1-CAPI-1.16.0). Conversions between both will be slow.
warnings.warn(
Map(center=[33.25, -7.57], controls=(ZoomControl(options=['position', 'zoom_in_text',␣
,→'zoom_in_title', 'zoom_o…

[3]: import numpy as np

import pandas as pd
import geopandas as gpd
import glob
from pykrige.ok import OrdinaryKriging
from pykrige.kriging_tools import write_asc_grid
import pykrige.kriging_tools as kt
import matplotlib.pyplot as plt
#from mpl_toolkits.basemap import Basemap
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.patches import Path, PathPatch
from shapely.geometry import Point, Polygon

df1 = pd.read_excel('data/PiezoBerchid.xlsx',sheet_name='06_2016')
x1 = df1['X']
y1 = df1['Y']
z1 = df1['Z']
h1 = df1['H']
df2 = pd.read_excel('data/PiezoBerchid.xlsx',sheet_name='10_2016')
x2 = df2['X']
y2 = df2['Y']
z2 = df2['Z']
h2 = df2['H']

Boundary_Shapefile = './data/NappeBerchidXY.shp'
boundary = gpd.read_file(Boundary_Shapefile)

hmin = h1.min();hmax = h1.max()

print(hmin,hmax)

xmin, ymin, xmax, ymax=275500, 275000, 335500, 325000

Ali HAMMANI, 2021

Chapitre 1 : Etude des nappes souterraines par le bilan hydrogéologique 10

# Générer une grille pour x et y pour 100 points dans chaque direction
grid_x = np.linspace(xmin, xmax, 100)
grid_y = np.linspace(ymin, ymax, 100)

# Exécuter la méthode de krigeage ordinare, variogram Gauss

OK = OrdinaryKriging(x1, y1, h1, variogram_model='linear', verbose=True,␣
,→enable_plotting=False,nlags=20) #gaussian

h1, ss1 = OK.execute('grid', grid_x, grid_y)

OK = OrdinaryKriging(x2, y2, h2, variogram_model='linear', verbose=True,␣

,→enable_plotting=False,nlags=20) #gaussian

h2, ss2 = OK.execute('grid', grid_x, grid_y)

# Résultat de l'interpolation spatiale

xintrp, yintrp = np.meshgrid(grid_x, grid_y)
fig, ax = plt.subplots(figsize=(30,30))

boundarygeom = boundary.geometry

dx=(xmax-xmin)/100
dy=(ymax-ymin)/100

dh=h2-h1
xx=xmin
V=0
S=0.05
A=dx*dy
N=0
for i in range(100):
yy=ymin
for j in range(100):
point= Point(xx,yy)
v=int(boundarygeom.contains(point))
#print(len(x),len(y))
if v==0:
h1[j][i]=np.NaN
h2[j][i]=np.NaN
dh[j][i]=np.NaN
else:
V += dh[j][i]
N += 1
yy+=dy
xx+=dx

V = A*S*V/1000000
print('La variation de la réserve est : %5.3f'% V, 'Mm3')

contour3 = plt.contour(xintrp, yintrp, dh,20)

plt.clabel(contour3, fmt = '%4.1f', colors = 'k', fontsize=18)
boundary.plot(ax=ax, color='white', alpha = 0.2, linewidth=5.5, edgecolor='black',␣
,→zorder = 5)

npts = len(x1)
plt.yticks(fontsize = 20)

# Charges hydraulique

Ali HAMMANI, 2021

Chapitre 1 : Etude des nappes souterraines par le bilan hydrogéologique 11

plt.title('Interpolation spatiale de la charqe hydraulique (%d points)' %␣

,→npts,fontsize = 30)

plt.show()

110.13999999999999 282.52
Adjusting data for anisotropy…
Initializing variogram model…
Coordinates type: 'euclidean'

Using 'linear' Variogram Model

Slope: 0.025846182547349147
Nugget: 234.81706617996386

Calculating statistics on variogram model fit…

Executing Ordinary Kriging…

Adjusting data for anisotropy…

Initializing variogram model…
Coordinates type: 'euclidean'

Using 'linear' Variogram Model

Slope: 0.026652772130290412
Nugget: 230.26343841083747

Calculating statistics on variogram model fit…

Executing Ordinary Kriging…

La variation de la réserve est : -24.400 Mm3

Ali HAMMANI, 2021

Chapitre 1 : Etude des nappes souterraines par le bilan hydrogéologique 12

[4]: import plotly

import plotly.graph_objects as go
import oceansdb

colorscale=[[0.0, 'rgb(20,29,67)'],
[0.1, 'rgb(28,76,96)'],
[0.2, 'rgb(16,125,121)'],
[0.3, 'rgb(92,166,133)'],
[0.4, 'rgb(182,202,175)'],
[0.5, 'rgb(253,245,243)'],
[0.6, 'rgb(230,183,162)'],
[0.7, 'rgb(211,118,105)'],
[0.8, 'rgb(174,63,95)'],
[0.9, 'rgb(116,25,93)'],
[1.0, 'rgb(51,13,53)']]

fig = go.Figure(data=[
go.Surface(z=h1, x=xintrp, y=yintrp, opacity=0.8, colorscale='Electric',␣
,→showscale=True),

go.Surface(z=h2, x=xintrp, y=yintrp, opacity=1.0, colorscale='Blues',␣

,→showscale=True),

#go.Surface(z=dh, x=xintrp, y=yintrp, opacity=1.0, colorscale='Blues',␣

,→showscale=True),

])
"""
go.Surface(
x=xintrp,
y=yintrp,
z=dh,
#surfacecolor=z2,
opacity=0.5,
cmin=-1,
cmax=1,
colorscale='RdBu',
colorbar=dict(len=0.3, x=0.8, y=0.5),
#figsize=(20,10),
),
"""

fig.update_traces(contours_z=dict(show=True, usecolormap=True,
highlightcolor="limegreen", project_z=True))

fig.update_layout(
autosize=True,
scene = dict(
zaxis = dict(nticks=4, range=[0,320],),
),
width=1000,
margin=dict(r=20, l=10, b=10, t=10)
)

fig.show()

Ali HAMMANI, 2021

Chapitre 1 : Etude des nappes souterraines par le bilan hydrogéologique 13

Bilan de la nappe de Béni Moussa

(Tadla)
[5]: import matplotlib
import matplotlib.pyplot as plt
import pandas as pd
from ipywidgets import *

def WaterBalance(iir,iip,ip,a,S):
DSc = (iir+iip)*Ai + 10*ip*Ap*a/1000000 + Al - Dr - NDr - DW - Agri
n=len(Period)
CumC.clear()
CumM.clear()
SumC=0
SumM=0
for i in range(0,n):
SumC += DSc[i]
CumC.append(SumC)
SumM += DSm[i]*S*a/100
CumM.append(SumM)

fig = plt.figure(figsize=(18,9))
ax = fig.add_subplot()
ax.plot(Period,CumM)
ax.plot(Period,CumC)
ax.set(xlabel='Période', ylabel='Variation de la réserve',title='Evolution de␣
,→lréserve mesurée et la réserve calculée')

plt.show()

df = pd.read_excel('data/BilanBM.xlsx')
Period=df['Period']
Ai = df['Irrigation']
Ap = df['Rainfall']
Al = df['LateralInflow']
Dr = df['ArtificialDrainage']
NDr = df['NaturalDrainage']
DW = df['DWPumping']
Agri = df['AgriPumping']
DSm = df['DS']

CumC=[]
CumM=[]

interact(WaterBalance,
iir=widgets.FloatSlider(value=0.1,min=0, max=0.5, step=0.005),
iip=widgets.FloatSlider(value=0.28,min=0, max=0.5, step=0.005),
ip=widgets.FloatSlider(value=0.2,min=0, max=0.5, step=0.005),
a=widgets.FloatSlider(value=69500,min=40000, max=100000, step=500),
S=widgets.FloatSlider(value=0.009,min=0, max=0.2, step=0.001,␣
,→layout=Layout(width='70%'))

interactive(children=(FloatSlider(value=0.1, description='iir', max=0.5, step=0.005),␣

,→FloatSlider(value=0.28, …

Ali HAMMANI, 2021

Chapitre 1 : Etude des nappes souterraines par le bilan hydrogéologique 14

[5]: <function main.WaterBalance(iir, iip, ip, a, S)>

Ali HAMMANI, 2021

Modélisationnumériquedes
Chapitre 2

écoulementssouterrains

Introduction
Définitions
a. Modèle
Un modèle est un outil désigné pour représenter une version simplifiée de la réalité.
En hydraulique souterraine, tout système pouvant reproduire la réponse ou le comportement d'un
réservoir souterrain est appelé modèle de celui-ci.
b. Simulation
La simulation est l'utilisation du modèle dans les applications réelles.
Si le modèle est bien élaboré, il peut être un outil de prédiction pour la gestion des ressources
en eau souterraine. Par exemple, en utilisant le modèle d'écoulement des eaux souterraines, il est
possible de tester différents schémas de gestion et de prédire l'effet de certaines actions. La validité
de la prédiction dépend des conditions d'approximation du modèle. Une bonne masse de données
est nécessaire pour une bonne prédiction par le modèle.

Différents types de modèles

Plusieurs types de modèles ont été utilisés pour l'étude des écoulements souterrains. Ils peuvent être
subdivisés en trois catégories:

Les modèles physiques

Le modèle le plus utilisé est le modèle de cuve de sable qui consiste en un réservoir rempli d'un
milieu poreux à travers lequel l'eau s'écoule. Ce sont des modèles de laboratoires. Cependant,
les phénomènes mesurés à l'échelle de la cuve de sable sont souvent différents des conditions
observées au champ, et les conclusions tirées de ces modèles doivent être corrigées quand on les
translate dans la réalité du terrain.

Les modèles analogiques

On dit qu'il y a analogie entre deux phénomènes physiques lorsqu'il y a une correspondance terme
à terme entre les paramètres physiques et les équations de base qui régissent les comportements
des deux systèmes physiques.

15
Chapitre 2: Modélisation numérique des écoulements souterrains 16

Il y a deux types d'analogie: analogie entre l'écoulement souterrain et le transfert de chaleur, et

l'analogie entre l'écoulement souterrain et le transfert du courant électrique dans un milieu conduc-
teur.

Variable Loi régissant le Equation générale

phénomène
Ecoulement Souterrain Charge h Loi de Darcy ∇2 h = 0
qx = −K ∆h
∆x

Transfert de chaleur Température T Loi de Fourrier ∇2 T = 0

qx = −K ∆T
∆x

Electricité Potentiel V Loi d’Ohm ∇2 V = 0

Ii = −K ∆V
∆x

L'analogie éléctrique est la plus utilisée.

Modèles mathématiques
Ils consistent à résoudre l'équation différentielle qui gouverne l'écoulement souterrain. La fiabilité
des prédictions par les modèles mathématiques dépend de la façon dont le modèle approche la
situation réelle. Des hypothèses simplificatrices doivent être faites pour construire le modèle du fait
que les situations sur le terrain sont plus complexes à être simulées exactement.
Il y a deux types de modèles mathématiques: les modèles analytiques et les modèles numériques.
• modèles analytiques se limitent aux cas simples et nécessitent beaucoup de simplifications dans
les cas complexes. Exemples: Formules de Dupuit, Jacob, Theis, ...etc.
• Pour traiter les situations réelles, il est souvent nécessaire de résoudre les modèles mathématiques
en utilisant les techniques numériques vu le développement des moyens de calcul. Les méthodes
les plus couramment utilisés pour la modélisation des problèmes d'écoulement souterrains sont
les méthodes de différences finies et les méthodes des éléments finies. La caractéristique de
ces deux méthodes est que le domaine d'écoulement est limité. Leur avantage est que les
paramètres hydrodynamiques peuvent être facilement variés à travers le système aquifère,
la formulation est convenable pour la modélisation des écoulements transitoires et que la
comparaison des résultats peut se faire directement.

Choix d’une technique numérique et

d’un code
On peut se voir obligé d'utiliser une solution numérique des équations de l'écoulement plutôt que
des solutions analytiques pour une ou plusieurs raisons:
1. Le domaine de l'écoulement est délimité par des limites complexes jouant un rôle pendant
le temps qui intervient dans la solution recherchée. Les solutions analytiques disponibles
s'appliquent à des milieux infinis ou semi-infinis; la méthode des images ne peut être utilisée, ou
bien il devient trop compliqué de représenter le rôle des limites par cette méthode.
2. Le problème est non linéaire (par exemple, la transmissivité varie avec la charge dans une
nappe libre), et in n'existe pas de solution analytique.
3. Les propriétés du milieu varient dans l'espace, tandis que les solutions analytiques supposent
que le milieu est homogène ou que la géométrie des hétérogénéités est très simple.
4. La géométrie et la grandeur du terme source sont trop complexes pour qu'elles puissent être
représentées par une source ponctuelle ou une ligne de source, ou une intégrale des deux,
suivant un tracé simple.

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 17

5. On peut trouver une solution analytique, mais son expression est si complexe (par exemple,
la somme d'une série infinie, l'intégrale de fonctions complexes) que les calculs numériques de
sa valeur demandent beaucoup plus d'effort (de programmation et de temps) que l'utilisation
directe d'une solution numérique du problème originel.
Dans les cas où des solutions numériques sont demandées, il faut d'abord décider: (i) quelle méthode
numérique choisir (essentiellement, différences finies, éléments finis, éléments limites), et (ii) comment
obtenir un code (le programmer ou obtenir l'accès à un code qui existe déjà ).
Il n'y a pas de réponse unanimement acceptée à la première question; au sujet des trois méthodes
décrites ci-dessus, on peut faire les remarques suivantes:
1. Différences finies : Cette méthode est facile à comprendre et à programmer. Elle convient très
bien à la résolution des problèmes régionaux d'écoulement des nappes, en une ou deux dimensions,
dans des systèmes multicouches ou en trois dimensions. Bien qu'elle soit, en principe, capable de
traiter des mailles de n'importe quelle forme et taille, elle est, en pratique, limitée à des mailles simples:
des carrés réguliers, des carrés gigognes, des rectangles ou des parallélépipèdes rectangulaires
en trois dimensions (figure 2.1). Elle peut très bien représenter les hétérogénéités des propriétés du
milieu, pourvu que la forme de ces hétérogénéités puisse être décrite de façon adéquate par la
forme des mailles; dans la pratique, l'anisotropie doit être limitée aux directions parallèles aux côtés
des mailles.

2. Eléments finis : Cette méthode est moins facile à expliquer et beaucoup moins facile à program-
mer que la précédente. Comme cette approche est plus flexible que celle des différences finies, un
programme d'éléments finis peut être plus compliqué à utiliser (davantage de données d'entrée, par
exemple, sur la géométrie des mailles, donc plus de possibilités d'erreurs) et peut demander davan-
tage de temps ordinateur. Cependant, la forme des mailles est beaucoup moins limitée: en pratique,
on prend des triangles et des quadrilatères en deux dimensions et, en trois dimensions, des tétraèdres
ou des parallélépipèdes de n'importe quel angle (figure 2.1). Ceci permet de décrire d'une manière
beaucoup plus satisfaisante la forme des limites du milieu ainsi que celle des hétérogénéités ou les
fonctions source, ce qui rend également la méthode éléments idéale pour résoudre les problèmes
à limites mobiles, par exemple, ceux ayant une surface libre et une interface abrupte entre eau
douce et eau salée ou entre deux fluides immiscibles. Elle est capable de traiter toutes les directions
d'anisotropie, et ces directions peuvent même changer d'un élément à un autre ou avec le temps.
En pratique, dans les problèmes d'écoulement, la méthode des éléments finis peut être employée
pour des études régionales, mais elle est pratiquement efficace dans les problèmes locaux de génie
civil tels que l'exhaure de l'eau d'une excavation, le drainage d'une mine et l'écoulement autour d'un
barrage, où les formes des limites et des hétérogénéités doivent être représentés avec précision.
Remarquez que quand il faut calculer les poussées en tant qu'entrée d'un modèle mécanique, il est
souvent nécessaire de le faire sur le même réseau que celui utilisé pour les calculs de structure, et
presque tous ceux-ci utilisent des éléments finis. Pour résoudre l'équation de transfert, la méthode
des éléments finis est bien supérieure à celle des différences finies.
3. Des approches d’éléments aux limites ou d’intégrales de limites ont été proposées récemment pour
résoudre l'équation d'écoulement. L'avantage principal de celles-ci est que la précision du calcul
ne dépend pas de la taille des éléments finis. Ainsi, on peut se servir de quelques éléments très

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 18

grands (ou même infinis), ce qui rend la méthode très efficace du point de vue du temps de calcul.
Dans un premier temps, la solution numérique se calcule uniquement le long des limites des éléments;
si l'on demande explicitement la solution à l'intérieur d'un élément, sa valeur est calculée dans un
second temps par une intégration numérique dans cet élément. La restriction principale est que les
propriétés du milieu, dans un élément donné, sont supposées constantes: si les hétérogénéités du
milieu sont telles qu'il faut inclure un grand nombre d'éléments afin de les décrire convenablement,
alors la méthode de l'intégrale des limites perd de sa supériorité, et l'on peut aussi bien se servir
de celle des différences finies ou de celle des éléments finis. Par conséquent, cette méthode est
beaucoup moins flexible et générale que les précédentes.
La deuxième question : comment obtenir un code relève plutôt du jugement personnel et on ne
peut qu'offrir quelques suggestions. La programmation d'un code simple de différences finies pour
un problème en une ou deux dimensions avec des mailles simples (carrés ou rectangle) et des limites
simples peut se faire en quelques jours. Toutefois, ce code ne sera pas facile à utiliser; si on veut
le rendre d'utilisation aisée (c'est-à -dire lui donner, par exemple, des données d'entrée simples, une
recherche d'erreurs et des massages d'erreur, des sorties graphiques), cela peut demander quelques
mois. Un code plus complexe en différences finies, multicouches ou à trois dimensions, ayant plusieurs
options (par exemple, des non linéarités) peut exiger de six à un an de travail, tout comme un code
en éléments finis à deux dimensions qui soit facilement maniable. Un code très complexe de transfert
en éléments finis à deux dimensions peut demander un effort allant d'un à deux ans, et un modèle de
gisement pétrolier à composantes et phases multiples à trois dimensions peut représenter un effort
de cinq à dix années (ou plus). Il ne faut pas oublier qu'un nouveau code doit être soigneusement
vérifié et validé à l'aide de solutions analytiques (ou numériques) connues avant qu'il puisse être
utilisé dans un contexte sérieux. Ces tests risquent d'être très longs.
Dans ce qui suit, nous décrivons brièvement la méthode des différences finies pour résoudre les
équations d'écoulement en nappe homogène et isotrope, et pour résoudre la même équation dans
le cas général d'écoulement transitoire en nappe homogène et anisotrope.

Méthode de différences finie en

aquifère homogène et isotrope
Ecoulement en régime permanent
Equation d’écoulement
L'équation différentielle pour un écoulement permanent dans un aquifère homogène et isotrope est:
∂2Φ ∂2Φ
+ 2 = −N (x, y) (2.1)
∂2x ∂ y

Où:

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 19

• Φ est le potentielle de charge:

– Φ = K H h pour une nappe en charge;
– Φ = 12 K h2 pour une nappe phréatique.
• N est e taux d'infiltration.

Discrétisation du domaine d’étude

La dicrétisation permet de remplacer l'équation aux dérivées partielles par un système d'équations
algébriques linéaires dans le cas des nappes captives (la transmissivité est constante), non linéaires
dans le cas des nappes phréatiques (la transmissivité est une fonction de la charge hydraulique).
Le domaine d'écoulement est découpé en mailles carrées (figure 2.3).

Les noeuds sont définis par les indices (i, j) où:

• x = (i − 1)a; y = (j − 1)a
• a = ∆x = ∆y
∂Φ
Aux points A et B ∂x peut être estimé par:
[ ∂Φ ] Φi,j −Φi−1,j
∂x A = a
(2.2)
[ ∂Φ ] Φi+1,j −Φi,j
∂x B = a

L'expression de la dérivée seconde de Φ par rapport à x au niveau du noeud (i, j) est obtenue
comme suit: [ 2 ]
∂ Φ Φi+1,j − 2Φi,j + Φi−1,j
2
= (2.3)
∂x i,j a2

∂2Φ
∂x2 est exprimée de la même façon et on obtient l'expression suivante pour le Laplacien:
[ 2 ] Φi−1,j + Φi+1,j + Φi,j−1 + Φi,j−1 − 4Φi,j
∆ Φ i,j = (2.4)
a2

En utilisant l'expression de l'équation différentielle (2.1) et en représentant N au noeud (i, j) par Ni,j ,
on obtient:
Φi−1,j + Φi+1,j + Φi,j−1 + Φi,j+1 − 4Φi,j = −a2 Ni,j (2.5)

On résout cette équation pour Φi,j :

1
Φi,j = [Φi−1,j + Φi+1,j + Φi,j−1 + Φi,j+1 + a2 Ni,j ] (2.6)
4

Cette équation est appliquée à chaque noeud et est utilisée avec les conditions aux limites pour
obtenir le potentiel hydraulique dans tous les noeuds.

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 20

Méthodes de résolution
La méthode simple pour la résolution du système d'équations (2.6) est la méthode de Jacobi, où
chaque valeur Φi,j est estimée en premier lieu à chaque noeud intérieur, et les valeurs des itérations
suivantes sont obtenues par l'application de l'équation (2.6) à chaque noeud en utilisant la valeur
calculée de l'itération précédente.
1 m
Φm+1 = [Φ + Φm m m 2
i+1,j + Φi,j−1 + Φi,j+1 + a Ni,j ] (2.7)
i,j
4 i−1,j

Où m est l'indice des itérations.

La méthode qui converge plus rapidement est la méthode de Gauss-Seidel. Dans cette méthode,
on commence par i = 2 et j = 2 et on balaye de gauche à droite et ligne par ligne comme si on lit
dans une page. De cette manière on peut utiliser les deux nouvelles valeurs calculées à l'itération
en cours, ainsi:
1
Φm+1
i,j = [Φm+1 + Φm m+1 m 2
i+1,j + Φi,j−1 + Φi,j+1 + a Ni,j ] (2.8)
4 i−1,j

Conditions aux limites

Deux types de conditions aux limites peuvent être appliquées à ce genre de problème :
i) Conditions aux limites de Derichlet: En terme de Φ
Dans ce cas de problème le potentiel est connu le long de la limite de l'aquifère, et par conséquent
à chaque noeud extrême.
∂Φ
ii) Conditions aux limites de Neumann: En terme de ∂n
∂Φ
Dans ces cas de problèmes, le flux est constant dans toute la limite du domaine. Afin d'estimer ∂n
aux noeuds extrêmes on ajoute des mailles fictives (figure 2.4).

Les valeurs de i et j sont nulles pour les noeuds extrêmes, qui sont appelés noeuds imaginaires. Pour
le cas ci-dessus, ∂Φ ∂Φ
∂n est donnée le long de la limite x = 0, i = 1 et ∂n pour le noeud (1, j) estimé
par: [ ] [ ]
∂Φ ∂Φ Φ2,j − Φ0,j
= = (2.9)
∂n 1,j ∂x 1,j 2a

Donc:
[ ]
∂Φ
Φ0,j = Φ2,j − 2a (2.10)
∂n 1,j

L'équation (2.10) est utilisée pour calculer la valeur de ∇2 Φ aux noeuds limites, l'expression (2.8)
devient: [ ] ]
1[ m ∂Φ
Φ1,j = Φ2,j + Φ1,j−1 + Φ1,j+1 + a N1,j − 2a
m+1 m m+1 2
(2.11)
4 ∂x 1,j

∂Φ
Les valeurs limites de ∂n sont ainsi appliquées par l'utilisation de (2.10) au lieu de (2.8).

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 21

Programme Python : Ecoulement unidimenionnelle en nappe captive

[10]: import numpy as np

import matplotlib.pyplot as plt
from ipywidgets import *

# DÃ©finition de la fonction
def Confined1D(N,h1,hn,L,kmax):
x = np.arange(0, L+L/N,L/N)
# Solution Analytique
phi_ana = (hn-h1)/L*x + h1
# Solution numÃ©rique
phi0_num = np.ones(N+1)*h1
phi0_num[0] = h1; phi0_num[N] = hn;
phi_num = phi0_num
Er=np.full(N+1, 9999)
k=0
while k<kmax: #Er>0.00001:
k += 1
for i in range(1,N):
phi_num[i] = (phi0_num[i+1] + phi_num[i-1]) / 2
phi0_num = phi_num
#Er=abs(phi_num[1:3]-phi0_num[1:3])
Er=abs(phi_num-phi_ana)
print("Erreur moyenne introduite : %10.5f" % Er.mean())

fig = plt.figure(figsize=(18,9))
ax = fig.add_subplot()
ax.plot(x,phi_ana)
ax.plot(x,phi_num)
ax.set(xlabel='x', ylabel='Charge hydraulique',title='Ecoulement␣
,→unidimentionnel en nappe captive')

ax.set_xticks(np.arange(0, L+L/N, L/N))

hmax=phi_ana.max()
plt.yticks(np.arange(10,hmax,1))
ax.grid()
plt.show()

Confined1D(10,15,10,100,20)

Erreur moyenne introduite : 0.30652

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 22

Programme Python : Ecoulement unidimensionnelle en nappe libre avec infiltration

[4]: import math

N=10;h1=15.;h2=10.;L=100.;I=0.00001;K=0.0001
phi = math.sqrt( -I/K*(-40**2 - L**2/4) - (h1**2-h2**2)/L*(-40) + (h1**2 + h2**2)/
,→2 )

print(phi)

24.9499498997493

[11]: import numpy as np

import matplotlib.pyplot as plt
from ipywidgets import *
import math

# DÃ©finition de la fonction
def Unconfined1D(N,h1,h2,L,I,K,kmax):
DeltaX=L/N
x = np.arange(-L/2, L/2+DeltaX,DeltaX)

# Solution Analytique
phi_ana=np.full(N+1, h1)
phi_ana[N] = h2
for i in range(1,N):
phi_ana[i] = math.sqrt( -I/K*(x[i]**2 - L**2/4) - (h1**2-h2**2)/L*x[i] +␣
,→(h1**2 + h2**2)/2 )

# Solution numÃ©rique
phi0_num = np.full(N+1, h1)
phi0_num[0] = h1; phi0_num[N] = h2;
phi_num = phi0_num
Er=np.full(N+1, 9999)
k=0
while k<kmax:
k += 1
for i in range(1,N):
Phi = (K/2*phi0_num[i+1]**2 + K/2*phi_num[i-1]**2 + I*DeltaX**2 ) / 2
phi_num[i] = math.sqrt(2*Phi/K)

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 23

phi0_num = phi_num
#Er=abs(phi_num[1:3]-phi0_num[1:3])
Er=abs(phi_num-phi_ana)
print("Erreur moyenne introduite : %10.5f" % Er.mean())
xd = -K/I/L/2*(h1**2-h2**2)
print("Point de partage des eaux x_d : %6.3f" % xd)

fig = plt.figure(figsize=(18,9))
ax = fig.add_subplot()
ax.plot(x,phi_ana, marker = '+')
ax.plot(x,phi_num)
ax.set(xlabel='x', ylabel='Charge hydraulique',title='Ecoulement␣
,→unidimentionnel en nappe captive')

ax.set_xticks(np.arange(-L/2, L/2+DeltaX, DeltaX))

hmax=phi_ana.max()
plt.yticks(np.arange(10,hmax,(hmax-10)/10))
ax.grid()
plt.show()

Unconfined1D(10,15.,15.,100.,0.00001,0.0001,20)

Erreur moyenne introduite : 0.50912

Point de partage des eaux x_d : -0.000

Ecoulement en régime transitoire

On discutera l'application de la méthode de différences finies aux problèmes d'écoulement transitoire
dans une nappe à surface libre. La formulation du problème dans un aquifère en charge est plus
simple, et il est facile d'adopter le schéma discuté ci-après aux problèmes d'écoulement en charge.

Equation de l’écoulement
L'équation différentielle pour un écoulement transitoire dans un aquifère à surface libre homogène
et isotrope est donnée par:
∂h
∆2 Φ = Sp −N (2.12)
∂t

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 24

Discrétisation du domaine d’étude

Il y a deux méthodes différentes pour discrétiser l'équation (2.12) par rapport au temps:
• La formulation explicite;
• La formulation implicite.
i. Formulation explicite
Dans ce cas la dicrétisation dans le temps est comme suit:
[ ] ht+∆t − ht
∆2 Φ t = Sp −N (2.13)
∆t

Cette expression est appelée explicite car le Laplacien est calculé au temps t de telle sorte que
l'expression de la charge au temps t + ∆t est obtenue comme suit:
∆t 2
ht+∆t = ht + [∆ Φ + N ] (2.14)
Sp

L'avantage de cette formulation est que la linéarisation de l'équation différentielle n'est pas néces-
saire. La charge à chaque noeud au temps t+∆t est obtenue par le calcul de ∇2 Φ par l'application
de l'équation (2.4) dans laquelle les valeurs de Φ au temps t sont utilisées.
Le schéma explicite a l'avantage d'être flexible et simple, mais il ne donne des résultats significatifs
que si ∆t est inférieur à une valeur critique ∆tc . (∆t < ∆tc ).
La valeur du pas de temps critique peut être déterminé par des tests d'erreurs: si ∆t excède ∆tc la
solution devient instable, et donne des résultats erronés.
La valeur de ∆tc peut être estimée par l'utilisation de l'expression suivante:
2
Sp (∆x)
∆tc = (2.15)
2Khmax

Où hmax est la charge maximale qui peut avoir lieu.

En pratique, il est intéressant de prendre un pas de temps inférieur à la valeur critique, par exemple
∆t = 21 ∆tc .
ii.Formulation implicite
Dans la formulation implicite on prend la valeur de ∇2 Φ entre t et t + ∆t. Le système résultant est
implicite en terme des inconnues et doit être résolu à chaque pas de temps.
Du fait de la nature implicite de la formulation, l'équation différentielle d'écoulement dans une nappe
à surface libre doit être linéarisée. La forme linéarisée de l'équation différentielle est donnée par:
Sp ∂Φ
∆2 Φ = −N (2.16)
kh ∂t

On représente le Laplacien comme moyenne entre les approximations aux temps t + n ∆t et to + (n +

1)∆t, où le temps to correspond au début de l'écoulement transitoire.

∆2 Φ = ω∆2 Φn+1 + (1 − ω)∆2 Φn (0 < ω < 1) (2.17)

On introduit la fonction Ωi,j :

1 n
Ωni,j = [Φ + Φni+1,j + Φni,j−1 + Φni,j+1 ] (2.18)
4 i−1,j

Donc: ∇2 Φ peut être écrite comme suit:

4
∆2 Φ = [ω(Ωni,j − Ωn+1
i,j ) + (1 − ω)(Ωi,j − Φi,j )]
n n
(2.19)
a2

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 25

La discrétisation de l'équation différentielle donne:

a2 Sp (Φi,j − Φni,j ) a2 n
n+1

i,j − Φi,j ) + (1 − w)(Ωi,j − Φi,j ) =

ω(Ωn+1 n+1 n n
N (2.20)
4Khni,j ∆t 4 i,j

Où h est prise égale à hni,j . Cette équation est résolue pour la valeur inconnue Φn+1
i,j , qui donne:

4Khni,j ∆t
Φn+1 i,j + (1 − ω)(Ωi,j − Φi,j )+
[ωΩn+1 n n
i,j =
a2 Sp + 4khni,j ω∆t
a2 Sp a2 n
n Φni,j + Ni,j ] (2.21)
4khi,j ∆t 4

n
Des améliorations peuvent être apportées en remplaçant Ni,j par une moyenne sur tout le pas du
temps, en utilisant la forme d'expression (2.17). L'équation (2.21) peut être résolue par l'utilisation de
la méthode itérative de Gauss-Seidel. Le schéma implicite est rendu explicite pour le cas où w = 0.
Pour w ̸= 0 le schéma est implicite; la formulation en différences finies est appelée régressive, au
contraire de la formulation progressive pour w = 0. Pour w = 1 le schéma est complètement implicite;
w = 0.5 la méthode est appelée schéma de différences finies centrale ou schéma de Crank-Nickolson.
L'avantage de la formulation implicite est que le pas de temps peut être choisi plus grand que celui
de la formulation explicite. Cependant, elle est complexe est la méthode itérative de Gauss-Seidel
doit être appliquée à chaque pas de temps.

Les conductions aux limites et initiales

Les conditions aux limites peuvent être du même type que celles citées précédemment. On aura
besoin en plus des conditions initiales qui peuvent être spécifiées de la façon suivante: Connaissons
la piézométrie en un ensemble de points de la nappe au temps initial (t = 0), on génère les charges
dans les noeuds restants à l'aide du régime permanent. Une fois l'état initial est connu, la simulation
du régime transitoire démarre.

Méthodes de différences finies en

aquifère hétérogène et anisotrope
Equation l’écoulement
L'équation différentielle partielle de l'écoulement transitoire dans une nappe captive hétérogène et
anisotrope peut être exprimée par:
∂ ∂h ∂ ∂h ∂ ∂h ∂ ∂h ∂h
(Txx ) + (Txy ) + (Tyx ) + (Tyy ) = S + W (x, y, t) (2.22)
∂x ∂x ∂x ∂y ∂y ∂x ∂y ∂y ∂t

Où:
• Txx , Txy , Tyx et Tyy sont les composantes du tenseur de transmissivité;
• Φ est Charge hydraulique;
• S est Coefficient d'emmagasinement;
• W (x, y, t) est le flux volumique de recharge ou de prélèvements par unité de surface de
l'aquifère.
Dans les modèles de simulation l'équation (2.22) est simplifiée en supposant que les axes des coor-
données Cartésiens x et y sont colinéaires avec les composantes principales des composantes de
transmissivité Txx et Tyy , ce qui donne:
∂ ∂h ∂ ∂h ∂h
(Txx ) + (Tyy ) = S + W (x, y, t) (2.23)
∂x ∂x ∂y ∂y ∂t

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 26

Pour une nappe phréatique la transmissivité est fonction de la charge. L'équation fondamentale
d'écoulement peut être exprimée comme:
∂ ∂h ∂ ∂h ∂h
(Kxx b ) + (Kyy b ) = Sy + W (x, y, t) (2.24)
∂x ∂x ∂y ∂y ∂t

Où:
• Kxx et Kyy sont les principales composantes du tenseur de la conductivité hydraulique;
• ωd est la porosité de drainage de la nappe libre;
• b est l'épaisseur saturée de l'aquifère.

Approximation de l’équation d’écoulement par les différences finies

Dans le but de résoudre les équations précédentes pour un aquifère hétérogène avec des limites
irrégulières, la région de l'aquifère est subdivisée en blocs réguliers dans lesquels les propriétés de
l'aquifère sont supposées uniformes. A chaque bloc de l'aquifère l'équation aux dérivées partielles
sera donc remplacée par une équation algébrique. On se trouvera donc avec N équations à N
inconnues (valeur de la charge aux noeuds) où N est le nombre de blocs représentant l'aquifère.
En utilisant discrétsation en blocs centrés (figure 2.5) dans laquelle les espacements sont variables
en fonction de la position dans l'aquifère, l'équation (2.23) peut être écrite comme suit:
( ) ( ) ( ) ( )
1 [ ∂h ∂h ] 1 [ ∂h ∂h ]
Txx − Txx + Tyy − Tyy
∆xi ∂x i,j+ 1 ∂x i,j− 1 ∆yj ∂y i+ 1 ,j ∂y i− 1 ,j
2 2 2 2

Si,j ( k )
= hi,j − hk−1
i,j + Wi,j
∆t
(2.25)

Dans laquelle:
• ∆xj est le pas d'espace dans la direction x pour la colonne j;
• ∆yi est le pas d'espace dans la direction y pour la ligne i;
• ∆t est le pas de temps t;
• i est l'indice dans la direction y;
• j est l'indice dans la direction x;
• k est le compteur de temps.

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 27

L'équation (2.25) peut être écrite aussi sous la forme:

1 [ hki,j+1 − hki,j hki,j − hki,j−1 ]

Txxi,j+ 1 − Txxi,j− 1
∆xi 2 ∆xj+ 21 2 ∆xj− 21
1 [ hki+1,j − hki,j hki,j − hki−1,j ]
+ Tyyi+ 1 ,j − Tyyi− 1 ,j
∆yj 2 ∆y i+ 12 2 ∆y i− 12
(
Si,j k )
= hi,j − hk−1
i,j + Wi,j
∆t
(2.26)

Dans laquelle:
• T xxi,j+ 12 est la transmissivité entre les noeuds (i, j) et (i, j + 1);
• ∆xi+ 12 est la distance entre les noeuds (i, j) et (i, j + 1).

L'équation (2.26)est écrite implicitement. En utilisant les notations F , D, H, et B l'équation (2.26)

peut être simplifiée en écrivant:
Si,j k
Fi,j (hki,j+1 − hki,j ) − Di,j (hki,j − hki,j−1 ) + Hi,j (hki+1,j − hki,j ) − Bi,j (hki,j − hki−1,j ) = (h − hk−1 (2.27)
∆t i,j i,j

Dans laquelle: [ ]
2Tyyi,j Tyyi−1,j
Tyyi,j ∆yi−1 +Tyyi−1,j ∆yi
Bi,j = (2.28)
∆yi
Tyyi,j Tyy[i−1,j] Tyyi−1,j
Le terme entre crochets est la moyenne harmonique de ∆yi , ∆yi−1 et de ∆yi−1

∆xi + ∆xi+1
Tyy 1
= ∆xi ∆xi+1
(2.29)
i+
2, j Ti,j + Ti+1,j

De la même façon: [ ]
2Txxi,j Txxi,j−1
Txxi,j ∆xj−1 +Txxi,j−1 ∆xj
Di,j = (2.30)
∆xj
[ 2Txxi,j Txxi,j+1
]
Txxi,j ∆xj+1 +Txxi,j+1 ∆xj
Fi,j = (2.31)
∆xj
[ 2Tyyi,j Tyyi+1,j
]
Tyyi,j ∆yi+1 +Tyyi+1,j ∆yi
Hi,j = (2.32)
∆yi

L'équation (2.27) est utilisée aussi pour l'approximation de l'équation (2.24) en définissant les trans-
missivités dans les équations (2.28), (2.30), (2.31) et (2.32) comme fonction de la charge de l'itération
précédente.
n
Txx i,j
= Kxxi,j bn−1k
i,j (2.33)

Après réarrangement l'équation (2.27) peut être écrite sous la forme suivante:

Bhki−1,j + Dhki,j−1 + Ehki,j + F hki,j+1 + Hhki+1,j = Q (2.34)

Dans laquelle:
( )
E =− B+D+F +H + S
∆t ;
(2.35)
Q= − ∆t
S k−1
h + Wi,j

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 28

Terme source W(x,y,t)

Le terme source W (x, y, t) comprend les débits de pompage, l'évapotranspiration, la recharge par
la pluviométrie ou/et par irrigation et la drainance à partir d'un aquifère soujacent.
Qkwi,j
− q ′ i,j + qet
k
k
Wi,j = − qre
k k
(2.36)
∆xj ∆yi i,j i,j

Où:
• Qw ki,j est le débit de pompage au niveau du puits situé au noeud (i, j, [L3 T −1 ];
• qre ki,j est le flux de recharge par unité de surface, [LT −1 ];

• q ′ i,j est le flux qui remonte de l'aquifère soujacent par unité de surface, [LT −1 ];
k

• qet ki,j est le taux d'évapotranspiration par unit de surface, [LT −1 ];

La drainance
La drainance à partir d'une nappe en charge sous-jacente peut être estimée par: (Bredefoeft and
Pinder, 1970).
( ) K ′ i,j
q ′ i,j = h0i,j − hki,j (
k
)1
πK ′ i,j t
3m2i,j Ssi,j
mi,j
  
∞
∑ ( )
  −n 2
 K ′ i,j
1 + 2 exp  ′ t
Ki,j  + ĥ0i,j − h0i,j (2.37)
n=1
mi,j
2m2i,j Ssi,j

Où:
• h0i,j est la charge hydraulique au noeud (i, j) au début de la période de pompage, [L];
• h0i,j est la charge au niveau de la nappe en charge drainant, [L];

• K ′ i,j est la conductivité hydraulique de la couche séparant les deux aquifère au niveau du
0

noeud (i, j), [L/T ];

• mi,j est l'épaisseur de cette couche au niveau du noeud (i, j), [L];
• Ssi,j est le coefficient d'emmagasinement spécifique au niveau de la nappe en charge, [L−1 ];
′
Ki,j t
• m2i,j Ssi,j
est un terme adimenssionnel;

• t est écoulé de la période de pompage.

Il faut noter que la drainance est la somme de deux termes, le premier considère l'effet transitoire;
le second correspond à la drainance permanente due au gradient initial à travers la couche
semi-perméable.

L’évapotranspiration
L'évapotranspiration est une fonction linéaire de la profondeur de la nappe elle est estimée par :

Qet hki,j ≥ Gi,j

k
qet i,j
= Qet − Qet
ET z (Gi,j − hki,j ) ET z > (Gi,j − hki,j ); hki,j > Gi,j (2.38)

0 [ET z ≤ (Gi,j − hki,j )]

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 29

Où:
• Qet est le taux maximal d'évapotranspiration, [LT −1 ];
• ETz est la profondeur au-delà de laquelle l'évapotranspiration cesse, [L];
• Gi,j est l'élévation de la surface du sol, [L].

Calcul de la charge au niveau d’un puits de pompage

Le débit au niveau du puits est donné par la formule de Theim (1906).

Qkwi,j πTi,j ∆h
= (2.39)
4 2 ln rrwe

Où:
k
Qwi,j
• 4 est le débit traversant une face de la maille;

Qkwi,j ∆h
= ∆xj Ti,j (2.40)
4 ∆y

• ∆h = hkw,j − hki,j ;
• Ti,j = T xxi,j = T yyi,j ;
r1
• re est un rayon fictif du puits: re = 4.81 où r1 = ∆xj ∆yi ;
• rw est la distance séparant le puits du point où la charge hw sera calculée qui peut être
écrite sous la forme;

Validité des solutions numériques

Dans la plupart des exemples, on vérifie la validité de la solution numérique en comparant les
valeurs de la charge générées par la méthode numérique à celles calculées par les solutions
analytiques. Cependant, les solutions analytiques ne sont disponibles que pour un nombre limité
de problèmes. Les méthodes numériques nous permettent de résoudre l'équation fondamentale qui
décrit l'écoulement souterrain à trois dimensions, pour des conditions aux limites complexes et pour les
aquifères hétérogènes et anisotropes, tandis que les solutions analytiques considèrent que l'aquifère
est homogène et isotrope. Pour les cas où la solution analytique est non disponible, les solutions
numériques peuvent être vérifiées par l'élaboration d'un bilan de masse. Cependant, les modèles
basés sur le bilan ne donnent pas beaucoup de précision à leur tours.
Le but principal des simulations numériques est de prédire les effets des schémas de gestion proposés
pour un système souterrain particulier. Le test final d'un modèle numérique est de vérifier s'il simule
réellement les situations réelles observées sur le terrain. Un tel modèle doit être calé et vérifié. Le
processus de calage nécessite à une identification des paramètres du modèle et/ou des conditions
aux limites, ce processus sera discuté au paragraphe suivant.

Calage des modèles numériques

La mise en place d'un modèle numérique est basée sur le détail géométrique de l'aquifère, les
informations concernant ses paramètres physiques, les limites, les entrées et les sorties...etc. Toutes
ces informations sont obtenues à partir des études géologiques et à partir des observations faites
sur le système aquifère. Quand une information est non disponible elle doit être estimée et vérifiée
après par le calage.

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 30

Le calage d'un modèle est le processus par lequel les paramètres de l'aquifère (qui peuvent aussi
inclure sa géométrie, les entrées, ...etc) sont déterminés, s'ils sont non disponibles, ou vérifiés s'ils sont
disponibles. Le calage est basé sur les données observées sur le terrain pendant un période donnée
du passé. Ces données sont, généralement, la piézométrie, les taux de pompage et de recharge, la
qualité de l'eau, les flux latéraux et verticaux ...etc. Le calage, ou l'identification des paramètres, est
connue sous le terme de ``problème inverse''.
Le calage consiste à faire des simulations d'un aquifère pour des périodes du passé pour lesquelles
les données sont disponibles dans le système aquifère réel (la piézométrie par exemple). Les résultats
calculés par cette simulation sont comparés à ceux observées. Le modèle est dit calé quand la
différence entre ces résultats est inférieure à une valeur spécifiée par le modélisateur. A la fin
de la phase de calage, on aura un modèle bien défini pour le système aquifère sous certaines
considérations. Tous les paramètres sont biens définis et le modèle pourra être utilisé pour la simulation
des opérations futures.

Résolution du problème inverse de

l’écoulement souterrain
Introduction
Les paramètres hydrauliques du modèle d'écoulement souterrain ne sont pas des grandeurs physiques
directement mesurables. L'élément majeur dans le processus de modélisation est l'estimation, à partir
des observations historiques de la charge hydraulique, des paramètres inconnus du système. His-
toriquement, les paramètres d'aquifères se déterminaient par des techniques basées sur des solutions
analytiques des équations gouvernant le système (Theis 1935). Ces méthodes sont généralement non
applicables aux paramètres distribués sur un système souvent caractérisé par son non homogénéité
et son anisotropie. En plus, les résultats obtenus ne sont que ponctuels, d'où un grand risque de non
représentativité du système en question .

Paramétrisation
Pour un aquifère non homogène, la dimension des paramètres est théoriquement infinie. Cependant,
dans le domaine d'application, les variables spatiales sont approchées par des méthodes de dif-
férences finies ou éléments finis. La réduction des paramètres d'une forme de dimension infinie à une
forme de dimension finie est appelée paramétrisation.
Il y a deux approches pour accomplir la paramétrisation : la première par zonation, et la deuxième
par interpolation .
a. Approche par zonation
C'est une approche, par laquelle le domaine d'écoulement est divisé en un certain nombre de sous-
régions ou zones, et un paramètre constant en valeur est utilisé pour caractériser chaque zone. Nous
obtenons un vecteur dont la dimension est le nombre total de zones.
b. Approche par interpolation
C'est une technique par laquelle, les paramètres inconnus dans les noeuds où aucune mesure n'a
été effectuée, sont générées par interpolation.

Identification des paramètres

L'identification des paramètres découle des propriétés physiques du milieu aquifère qui cumule les
fonctions de conducteurs et de réservoirs d'eau. Deux paramètres sont généralement considérés :
la transmissivité et le coefficient d'emmagasinement.

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 31

Analyse de sensibilité d’un modèle

L'analyse de sensibilité d'un modèle a pour but d'étudier les répercussions de faibles variations des
paramètres du modèle sur les résultats. De façon plus précise, elle se limite essentiellement à l'étude
des variations sur chacun des coefficients ou paramètres qui n'entraà®nent pas des changements
importants.
Cette analyse est très utile, notamment, lorsque les données du problème ne sont pas connues avec
précision, et permet éventuellement de repérer les coefficients qu'il y a lieu de déterminer avec plus
de précision et ceux qui, au contraire, peuvent être estimés plus grossièrement .
La solution optimale, trouvée par la résolution inverse peut être très sensible aux erreurs d'estimation
des paramètres servant comme données de base du modèle. Donc, avant de passer à l'utilisation,
il faut chercher dans quelle mesure la solution optimale est stable .
Pour chaque paramètre, on fixe un pas de variation, soit pour sa surestimation, soit pour sa sous-
estimation, et on analyse la variation de la fonction objectif .
Une surestimation (ou sous-estimation) d'un paramètre, ou deux paramètres peut affecter le résultat.
C'est pour cette raison que nous n'allons pas nous contenter de varier seul chaque paramètre
pour voir le degré de stabilité de la solution optimale, mais nous aurons à considérer toutes les
combinaisons possibles de variation des différents paramètres .

Ali HAMMANI, 2021

Modélisationnumériquedes
Chapitre 3

écoulementssouterrains

Différents types de modèles

Plusieurs types de modèles ont été utilisés pour l'étude des écoulements souterrains. Ils peuvent être
subdivisés en trois catégories:

Les modèles physiques

Les modèles analogiques

32
Chapitre 2: Modélisation numérique des écoulements souterrains 33

Il y a deux types d'analogie: analogie entre l'écoulement souterrain et le transfert de chaleur, et

l'analogie entre l'écoulement souterrain et le transfert du courant électrique dans un milieu conduc-
teur.

Variable Loi régissant le Equation générale

phénomène
Ecoulement Souterrain Charge h Loi de Darcy ∇2 h = 0
qx = −K ∆h
∆x

Transfert de chaleur Température T Loi de Fourrier ∇2 T = 0

qx = −K ∆T
∆x

Electricité Potentiel V Loi d’Ohm ∇2 V = 0

Ii = −K ∆V
∆x

L'analogie éléctrique est la plus utilisée.

Choix d’une technique numérique et

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 34

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 35

Méthode de différences finie en

Où:

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 36

• Φ est le potentielle de charge:

– Φ = K H h pour une nappe en charge;
– Φ = 12 K h2 pour une nappe phréatique.
• N est e taux d'infiltration.

Discrétisation du domaine d’étude

Les noeuds sont définis par les indices (i, j) où:

• x = (i − 1)a; y = (j − 1)a
• a = ∆x = ∆y
∂Φ
Aux points A et B ∂x peut être estimé par:
[ ∂Φ ] Φi,j −Φi−1,j
∂x A = a
(3.2)
[ ∂Φ ] Φi+1,j −Φi,j
∂x B = a

L'expression de la dérivée seconde de Φ par rapport à x au niveau du noeud (i, j) est obtenue
comme suit: [ 2 ]
∂ Φ Φi+1,j − 2Φi,j + Φi−1,j
2
= (3.3)
∂x i,j a2

∂2Φ
∂x2 est exprimée de la même façon et on obtient l'expression suivante pour le Laplacien:
[ 2 ] Φi−1,j + Φi+1,j + Φi,j−1 + Φi,j−1 − 4Φi,j
∆ Φ i,j = (3.4)
a2

En utilisant l'expression de l'équation différentielle (2.1) et en représentant N au noeud (i, j) par Ni,j ,
on obtient:
Φi−1,j + Φi+1,j + Φi,j−1 + Φi,j+1 − 4Φi,j = −a2 Ni,j (3.5)

On résout cette équation pour Φi,j :

1
Φi,j = [Φi−1,j + Φi+1,j + Φi,j−1 + Φi,j+1 + a2 Ni,j ] (3.6)
4

Cette équation est appliquée à chaque noeud et est utilisée avec les conditions aux limites pour
obtenir le potentiel hydraulique dans tous les noeuds.

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 37

Où m est l'indice des itérations.

Conditions aux limites

Donc:
[ ]
∂Φ
Φ0,j = Φ2,j − 2a (3.10)
∂n 1,j

L'équation (2.10) est utilisée pour calculer la valeur de ∇2 Φ aux noeuds limites, l'expression (2.8)
devient: [ ] ]
1[ m ∂Φ
Φ1,j = Φ2,j + Φ1,j−1 + Φ1,j+1 + a N1,j − 2a
m+1 m m+1 2
(3.11)
4 ∂x 1,j

∂Φ
Les valeurs limites de ∂n sont ainsi appliquées par l'utilisation de (2.10) au lieu de (2.8).

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 38

Programme Python : Ecoulement unidimenionnelle en nappe captive

[10]: import numpy as np

import matplotlib.pyplot as plt
from ipywidgets import *

ax.set_xticks(np.arange(0, L+L/N, L/N))

hmax=phi_ana.max()
plt.yticks(np.arange(10,hmax,1))
ax.grid()
plt.show()

Confined1D(10,15,10,100,20)

Erreur moyenne introduite : 0.30652

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 39

Programme Python : Ecoulement unidimensionnelle en nappe libre avec infiltration

[4]: import math

N=10;h1=15.;h2=10.;L=100.;I=0.00001;K=0.0001
phi = math.sqrt( -I/K*(-40**2 - L**2/4) - (h1**2-h2**2)/L*(-40) + (h1**2 + h2**2)/
,→2 )

print(phi)

24.9499498997493

[11]: import numpy as np

import matplotlib.pyplot as plt
from ipywidgets import *
import math

# DÃ©finition de la fonction
def Unconfined1D(N,h1,h2,L,I,K,kmax):
DeltaX=L/N
x = np.arange(-L/2, L/2+DeltaX,DeltaX)

# Solution Analytique
phi_ana=np.full(N+1, h1)
phi_ana[N] = h2
for i in range(1,N):
phi_ana[i] = math.sqrt( -I/K*(x[i]**2 - L**2/4) - (h1**2-h2**2)/L*x[i] +␣
,→(h1**2 + h2**2)/2 )

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 40

ax.set_xticks(np.arange(-L/2, L/2+DeltaX, DeltaX))

hmax=phi_ana.max()
plt.yticks(np.arange(10,hmax,(hmax-10)/10))
ax.grid()
plt.show()

Unconfined1D(10,15.,15.,100.,0.00001,0.0001,20)

Erreur moyenne introduite : 0.50912

Point de partage des eaux x_d : -0.000

Ecoulement en régime transitoire

Equation de l’écoulement
L'équation différentielle pour un écoulement transitoire dans un aquifère à surface libre homogène
et isotrope est donnée par:
∂h
∆2 Φ = Sp −N (3.12)
∂t

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 41

Discrétisation du domaine d’étude

Où hmax est la charge maximale qui peut avoir lieu.

On représente le Laplacien comme moyenne entre les approximations aux temps t + n ∆t et to + (n +

1)∆t, où le temps to correspond au début de l'écoulement transitoire.

∆2 Φ = ω∆2 Φn+1 + (1 − ω)∆2 Φn (0 < ω < 1) (3.17)

On introduit la fonction Ωi,j :

1 n
Ωni,j = [Φ + Φni+1,j + Φni,j−1 + Φni,j+1 ] (3.18)
4 i−1,j

Donc: ∇2 Φ peut être écrite comme suit:

4
∆2 Φ = [ω(Ωni,j − Ωn+1
i,j ) + (1 − ω)(Ωi,j − Φi,j )]
n n
(3.19)
a2

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 42

La discrétisation de l'équation différentielle donne:

a2 Sp (Φi,j − Φni,j ) a2 n
n+1

i,j − Φi,j ) + (1 − w)(Ωi,j − Φi,j ) =

ω(Ωn+1 n+1 n n
N (3.20)
4Khni,j ∆t 4 i,j

Où h est prise égale à hni,j . Cette équation est résolue pour la valeur inconnue Φn+1
i,j , qui donne:

4Khni,j ∆t
Φn+1 i,j + (1 − ω)(Ωi,j − Φi,j )+
[ωΩn+1 n n
i,j =
a2 Sp + 4khni,j ω∆t
a2 Sp a2 n
n Φni,j + Ni,j ] (3.21)
4khi,j ∆t 4

Les conductions aux limites et initiales

Méthodes de différences finies en

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 43

Approximation de l’équation d’écoulement par les différences finies

Si,j ( k )
= hi,j − hk−1
i,j + Wi,j
∆t
(3.25)

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 44

L'équation (2.25) peut être écrite aussi sous la forme:

1 [ hki,j+1 − hki,j hki,j − hki,j−1 ]

Txxi,j+ 1 − Txxi,j− 1
∆xi 2 ∆xj+ 21 2 ∆xj− 21
1 [ hki+1,j − hki,j hki,j − hki−1,j ]
+ Tyyi+ 1 ,j − Tyyi− 1 ,j
∆yj 2 ∆y i+ 12 2 ∆y i− 12
(
Si,j k )
= hi,j − hk−1
i,j + Wi,j
∆t
(3.26)

Dans laquelle:
• T xxi,j+ 12 est la transmissivité entre les noeuds (i, j) et (i, j + 1);
• ∆xi+ 12 est la distance entre les noeuds (i, j) et (i, j + 1).

L'équation (2.26)est écrite implicitement. En utilisant les notations F , D, H, et B l'équation (2.26)

peut être simplifiée en écrivant:
Si,j k
Fi,j (hki,j+1 − hki,j ) − Di,j (hki,j − hki,j−1 ) + Hi,j (hki+1,j − hki,j ) − Bi,j (hki,j − hki−1,j ) = (h − hk−1 (3.27)
∆t i,j i,j

Dans laquelle: [ ]
2Tyyi,j Tyyi−1,j
Tyyi,j ∆yi−1 +Tyyi−1,j ∆yi
Bi,j = (3.28)
∆yi
Tyyi,j Tyy[i−1,j] Tyyi−1,j
Le terme entre crochets est la moyenne harmonique de ∆yi , ∆yi−1 et de ∆yi−1

∆xi + ∆xi+1
Tyy 1
= ∆xi ∆xi+1
(3.29)
i+
2, j Ti,j + Ti+1,j

De la même façon: [ ]
2Txxi,j Txxi,j−1
Txxi,j ∆xj−1 +Txxi,j−1 ∆xj
Di,j = (3.30)
∆xj
[ 2Txxi,j Txxi,j+1
]
Txxi,j ∆xj+1 +Txxi,j+1 ∆xj
Fi,j = (3.31)
∆xj
[ 2Tyyi,j Tyyi+1,j
]
Tyyi,j ∆yi+1 +Tyyi+1,j ∆yi
Hi,j = (3.32)
∆yi

Après réarrangement l'équation (2.27) peut être écrite sous la forme suivante:

Bhki−1,j + Dhki,j−1 + Ehki,j + F hki,j+1 + Hhki+1,j = Q (3.34)

Dans laquelle:
( )
E =− B+D+F +H + S
∆t ;
(3.35)
Q= − ∆t
S k−1
h + Wi,j

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 45

Terme source W(x,y,t)

Où:
• Qw ki,j est le débit de pompage au niveau du puits situé au noeud (i, j, [L3 T −1 ];
• qre ki,j est le flux de recharge par unité de surface, [LT −1 ];

• q ′ i,j est le flux qui remonte de l'aquifère soujacent par unité de surface, [LT −1 ];
k

• qet ki,j est le taux d'évapotranspiration par unit de surface, [LT −1 ];

Où:
• h0i,j est la charge hydraulique au noeud (i, j) au début de la période de pompage, [L];
• h0i,j est la charge au niveau de la nappe en charge drainant, [L];

• K ′ i,j est la conductivité hydraulique de la couche séparant les deux aquifère au niveau du
0

noeud (i, j), [L/T ];

• t est écoulé de la période de pompage.

L’évapotranspiration
L'évapotranspiration est une fonction linéaire de la profondeur de la nappe elle est estimée par :

Qet hki,j ≥ Gi,j

k
qet i,j
= Qet − Qet
ET z (Gi,j − hki,j ) ET z > (Gi,j − hki,j ); hki,j > Gi,j (3.38)

0 [ET z ≤ (Gi,j − hki,j )]

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 46

Calcul de la charge au niveau d’un puits de pompage

Le débit au niveau du puits est donné par la formule de Theim (1906).

Qkwi,j πTi,j ∆h
= (3.39)
4 2 ln rrwe

Où:
k
Qwi,j
• 4 est le débit traversant une face de la maille;

Qkwi,j ∆h
= ∆xj Ti,j (3.40)
4 ∆y

Validité des solutions numériques

Calage des modèles numériques

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 47

Résolution du problème inverse de

Identification des paramètres

Ali HAMMANI, 2021

Chapitre 2: Modélisation numérique des écoulements souterrains 48

Analyse de sensibilité d’un modèle

Institut
Agronomique et Vétérinaire
Hassan II

Institut Agronomique et Vétérinaire Hassan II

Département du Génie Rural

Ali HAMMANI, 2021

Chapitre 4

Initiationà GMS

MODFLOW – L’approche par maillage

Deux approches peuvent être utilisées pour construire des simulations Modflow cans GMS : l'approche
par maillage st l'approche conceptuel. L'approche par maillage implique le travail avec des mailles
3D en appliquant les termes puits/sources et d'autres modèles bloc par bloc. L'approche con-
ceptuelle implique l'utilisation d'un outil SIG pour élaborer une carte pour développer le modèle
conceptuel du site à modéliser. Les données sont copiées après dans la grille du maillage.
L'approche par maillage est décrite dans se présent chapitre. Dans la plus part des cas, l'approche
conceptuelle eat plus efficace que l'approche par maillage. Cependant, cette dernière est très
pratique pour les problèmes simples ou pour des exercices académiques où l'édition des données
par maille est nécessaire.

Description du problème
Le problème à résoudre dans cet exercice est donné dans la figure 1. Trois aquifères vont être
simulés en utilisant trais couches dans une grille conceptuelle de blocs. na grille couvre une région
carrée de 75 km sur 75 km. Le maillage consiste en 15 lignes et 15 colonnes, les mailles sont carrées
et mesurent 5000 m de coté dans le plan horizontal. Pour un souci de simplicité la topographie
de chaque couche est considéré plate. Les valeurs re conductivité hydraulique sont données dons
la direction horizontale. Pour la direction verticale on utilisera une fonction de la conductivité
hydraulique horizontale.

49
Initiation à GMS 50

Pour chaque couche les données sont comme suit :

• Couche 1: K = 15 m/j, côte du niveau du sol = 60 m, côte de substratum = -45 m;
• Couche 2: K = 0.9 m/j, côte du toit = -45 m, côte de mur = -120m;
• Couche 3: K = 2 m/j, côte du toit = -120 m, côte de mur = -215.
L'écoulement dans le système est dû à l'infiltration des précipitations qui sera défini comme une
recharge de la nappe. Les prélèvements du système sont ceux qui se font vers un drain, un puits de
pompage (non visible sur la figure 1) et vers un lac représenté par une limite à charge hydraulique
constante à gauche de la figure. Les charges initiales sont considérées égales à 0 et on cherche
à trouver la solution mu régime permanent.

Mise en route
1. SF nécessaire, lancer GMS. Si GMS est déjà lancé, sélectionnée File | New dans le menu pour
assurer que la configuration du programme revient a son état initial.

Unités
A ce stade, on peut définir les unités utilisées dans le modèle. Les unités choisies seront modifiées en
éditant les champs de l'interface GMS.
1. Sélectionner la commande du menu Edit | Units ;
2. Pour Length, enter m (pour mètre). Poor Time, entrer d (pour jours). On ignore pour l'instant les
autres unités (elles ne sont pas utilisées dans les simulations de l'écoulement).
3. Cliquer sur le bouton OK.

Création du maillage
La première étape de résolution du problème est de créer un maillage de différences finies en 3D.
1. Dans l'espace Explorateur de Projet (Project Explorer) cliquer par le bouton droit de souris dans
une zone vide puis dans le menu flottant choisir New | 3D Grid.
2. Dans lm section intitulée X-dimension, entrer la valeur 22860 pour la longueur (Length value)
et 15 pour le numéro de mailles (Number cells) ;
3. Dans la section intitulée Y-dimension, entrer la valeur 22860 pour la longueur (Length value)
et 15 pour le numéro de mailles (Number cells) ;
4. Dans la section intitulée Z-dimension, entrer 3 pour le nombre de mailles (Number cells) ;
Plus tard, nous allons entrer les valeurs des côtes du toit et mur de chaque couche. Ainsi, l'épaisseur
des mailles dans la direction z introduite à ce niveau ne vont pat affecter les calculs du modèle
MODFLOW.
1. Cliquer sur le bouton OK.
Le maillage va apparaître sur l'écran de l'ordinateur. Une représentation simplifiée du maillage doit

également apparaître dans Mini-Grid Plot. .

Création d’une Simulation MODFuOW

L'étape suivante dans la mise en oeuvre du modèle est s'initialiser la simulation de MODFLOW.
1. Dans l'espace Explorateur de Projet (Project Explorer) cliquer par le bouton droit sur le dossier
3D Grid cata puis sélectionner New MODFLlW dans le menu flottant.

Ali HAMMANI, 2021

Initiation à GMS 51

Le Package global
Les données d'entrée de MODFLOW sont subdivisées en packages. Certains packages sont option-
nels et d'autres sont obligatoires. L'un des packages obligatoire est le Package Global par lequel
nous allons commencer.
En premier lieu nous allons choisir le package.
1. Cliquer sur le bouton Packages.

Les packages
La boîte de dialogue packages est utilisés pour spécifier quel package va être utilisé pour la mise
el œuvre du modèle. Le package de baae est déjà utilisé et né peut pas être arrête. Pour choisir
d'autres packages il faut passer par les étapes suivantes:
1. Dans la section Boundary Conditions cocher les options Drain (DRN1) et pour la prise en compte
du drainage, Well (WEL1) pour la prise en compte des pompages.
2. Pour la prise en compte de la recharge choisir l'option Recharge (RCH1).
3. Dans la section Solver, choisir le package Strongly Implicite Procedure (SIP1) pour choisir la
méthode de résolution des équations de différences finies.
4. Cliquer sur OK pour quitter la boîte de dialogue Packages.

Le tableau IBOUND
Dans l'étape suivante nous allons spécifier le tableau IBOUND (variable qui prend des valeurs <0,
=0 ou >1 en fonction du type de la condition aux limites). Quand IeOUND>0 la maille est dite
active, quand IBOUrD=0 la maille est inactive et quand IBOUND<0 la maille est à charge imposée.
Pour notre problème, toutes les mailles seront actives excepté pour la colonne de gauche des deux
premières couches dans les mailles où la charge hydraulique est imposée.
1. Cliquer sur le bouton IBOUND .
La boîte ds dialogue IBOUND affiche les valeurs du tableau IBOUND pour chaque louche à la fois.
Le champ Layer dans le coin haut à gauche est utilisé pour changer la couche en cours. Pour notre
problème, nous allons spécifier que toutes les valeurs du tableau sont supérieures à zéro, excepté
pour la colonne gauche des couches 1 et 2 où les valeurs sont inférieures à 0. Par défaut, les
valeurs du tableau sont supérieures à zéro. Ainsi, ce dont on a besoin est de changer les valeurs
des mailles où la charge hydraulique est imposée. Ceci peut être accompli en introduisant la valeur
-1 pour chaque maille à charge imposée. Cependant, il existe une autre façon pour éditer le
tableau IBOUND qui est plus simple dans ce cas. cette méthode vu être décrite plus tard. Pour
l'instant on va garder toutes les mailles active (IBOUND=1).
1. Cliquer sur bouton OK pour quitter la boîte de dialogue IBOUND.

Charges initiales
Dans l'étape suivante on va spécifier les valeurs des charges initiales.
1. Cliquer sur le bouton Strating Heads.
Le tableau des charges hydrauliques est utilisé pour établir une valeur initiale en exécutant une
simulation en régime transitoire. Puisque nous allons exécuter une simulation en régime permanent,
ces charges initiales ne vont pas changer la solution. Cependant, plus les charges initiales sont près
des valeurs de la solution plus MODFLOW converge rapidement. En outre, pour certains types de
couches, si les valeurs principales de départ sont très faibles MODEFLOW peut l'interpréter par des
mailles qui deviennent sèches. Pour le problème qu'on est entrain de traiter des valeurs de la charge
hydraulique égales à 0 sont suffisantes.
Le tableau des charges hydrauliques initiales est aussi utilisé pour spécifier les valeurs de la charge
hydraulique des mailles à charge imposées. Pour notre problème, les charges hydrauliques imposées

Ali HAMMANI, 2021

Initiation à GMS 52

sont égales à 0. On a n'a pas donc besoin de changer ces valeurs puisqu'elles sont déjà égale
à 0.
1. Cliquer sur OK pour quitter la boîte de dialogHe Starting heads.

Les côtes du toit et du mur de chaque couche

Dans l'étape suivante nous allons spécifier les côtes du toit et du mur de chaque couche.
1. Cliquer sur le bouton Top Elevation.
2. S'assurer que la couche en cours en 1.
3. Cliquer sur le bouton Constante Layer.
4. Entrer la valeur 60 et cliquer sur OK.
5. Cliquer sur OK pout quitter la boîte de dialogue Top Elevations.
GMS force les côtes du toit de la couche d'être à la même côte que le mur de la couche en dessus.
Ainsi, nous avons juste besoin d'entrer les côtes des murs des autres couches et les toits vont être
spécifiés automatiquement.
1. Cliquer sur le bouton Bottom Elevation.
2. S'assurer que la couche en cours est 1.
3. Cliquer sur le bouton Constante Layer.
4. Entrer la valeur -45 et clique sur OK.
5. Changer Ca coche Layer à 2.
6. Cliquer sur le bouton Constante Layer.
7. Entrer la valeur -120 et cliquer sur OK.
8. Change the Layer to 3.
9. Cliquer sur le bouton Constante layer.
10. Entrer la valeur -215 et clique sur OK.
11. Cliquer sur OK pour sortir de la boîte de dialogue Bottom Elevation.
12. Cliquer sur OK pour sortir de la boîte de dialoogue MODFLOW Global Package.

Affectation des valeurs à IBOUND directement aux mailles

Comme mentionner plus haut, les valeurs de IBOUND peuvent être entrées à travers la boîte de
dialogue du Tableau IBOUND. Dans certains cas, il est plus simple d'affecter les valeurs directement
aux mailles. Ceci peut être accompli en utilisant la commande propriétés de la maille (Cell Properties).
Avant d'utiliser la commande, nous devons d'abord sélectionner les mailles de la colonne gauche
des deux couches de dessus.

Affichage de la colonne gauche

Pour simplifier la sélection des mailles, nous allons changer l'affichage de telle manière à voir la
couche gauche du maillage.

1. Cliquer sur le bouton Side View

Le maillage apparait très fin. Pour faciliter les choses, nous allons amplifier la valeur de Z de telle
manière à faire apparaître les mailles dana la directioo Z.

1. Clique sur le bouton Display .

2. Changer Z magnification à 15.

Ali HAMMANI, 2021

Initiation à GMS 53

3. Cliquer sur OK.

Sélection des mailles

Pour sélectionner une maille :

1. Cliquer sur l'outil Select lells .

2. Changer le numéro de colonne à 1 dans l'affichage Miai-Grid.
Il faut noter que nous affichons la colonne 1 à gauche du domaine.
1. Dessiner un rectangle autour de toutes les mailles des deux premières colonnes en haut du
maillage.

Changement des valeurs de IBOUND

Pour éditer les valeurs de IBOUND:
1. Cliquer par le bouton droit sur l'une des mailles sélectionnées.
2. Choisir Properties dans le menu flottant.
3. Changer l'option IBOUND à Specified head.
4. Cliquer sur OK poor sortir de la boîte de dialogue 3D Cell Properties.

5. Cliquer Plan View .

Noter qu'un symbole est affiché dans les mailles que nous avons éditées, indiquant que ces mailles
sont des mailles à charge hydraulique imposée.

Vérification des laveurs

Pour s'assurer que les valeurs de IBOUND ont été entrées correctement:
1. Exécuter la commande Du menu MOdFLOW | Global Options.
2. Cliquer sur lu bouton IBOUND.
3. Changer la valeur de la couche (Layer) en haut à gauche de la boîte de dialogue.
Notez que les mailles de la colonne gauche ont toutes une valeur de IBOU N D = −1 dans les deux
couches de dessus. Toutes ses données d'entrée de MODFLOW peuvent être éditées dans GMS en
utilisant à la fois le tableau de données ou par sélection d'une série de maille directement.
1. Cliquer sur OK pour quitter la boîte de dialogue du tableau IBOUND.
2. Cliquer sur OK pour quitter la boîte de dialogue MODFLOW Global Package.

Le package LPF
La prochaine étape dans la mise en œuvre du modère est d'entrer les données sur le package
des propriétés de l'écoulement dans chaque couches (Layer Property Flow ou LPF). Le package LPF
calcule la conductance entre les mailles et configure les équations de différences finies décrivant
l'écoulement d'une maille à l'autre. Pour entrer les données LPF:
1. Lancer le menu MODFLOW | LPF Package.

Ali HAMMANI, 2021

Initiation à GMS 54

Types de couches
Les options dans la section Layer Data de la boîte de dialogue sont utilisées pour définir le type
de la couche et la conductivité hydraulique pour chaque couche. Pour notre problème, nous avons
trois couches. La couche supérieure est une nappe phréatique les autres couches sont des nappes
captives. Le couche par défaut dans GMS est `convertible' c'est à dire elle peut être captive ou
libre. Ainsi nous n'avons pas besoin de changer le type de la couche.

Les paramètres des couches

Les boutons dans la section Layer Data sont utilisés pour entrer les paramètres nécessaires pour
calculer les conductances d'une maille à l'autre. MODFLOW requière un jeu de paramètres pour
chaque couche dépendant du type de la couche.

La couche supérieure
En premier lieu, nous allons entrer les données pour la couche supérieure:
1. Cliquez sur le bouton Horizontal Hydraulic Conductivity.
2. Cliquez sur le bouton Constante Layer.
3. Entrer la valeur 15.
4. Cliquer uur OK.
5. Cliquer sur OK pour sortir de la boîte de dialogue Horizontal Hydraulic Conductivity.
6. Répéter ce processus pour entrer la valeur 10 pour l'anisotropie verticale.

Couche du milieu
Dans la suite nous allons introduire les données pour la couche du milieu:
1. Changer la couche à 2 pour le champ layer pour pouvoir éditer les données de la couche 2.
2. Entrer les valeurs suivantes pour la couche 2:

paramètre Value
Conductivité Hydraulique Horizontale 0.9 m/j
Anisotropie verticale 5

Couche inférieure
Nous introduisant maintenant les données pour la couche inférieure:
1. Changer la couche layer à 3 et entrer les valeurs suivantes:

Paramètre Value
Conductivité Hydraulique Horizontale 2 m/j
Anisotropie verticale 5
1. Cliquer sud OK pour quitter boîte de dialogue LPF Package.

Le package recharge
Dans la suite nous allons entrer lea données sur la recharge en utilisant le package recharge. Le
package recharge est utilisé pour simuler la recharge d'un aquifère due à la pluie et l'infiltration.
Pour entrer les données sur la recharge:

Ali HAMMANI, 2021

Initiation à GMS 55

1. Exécuter le menu MODFLOW | Source/Sink Packages | Recharge Package.

2. Cliquer sur Constante Array.
3. Entrer la valeur 0.0009 et cliquer sur OK.
4. Cliquer sur OK pour sortir de la boîte de dialogue Recharge Package.

Le package drain
Nous allons maintenant définir la ligne du drain dans la couche supérieure du modèle. Pour définir
le drain, il faut tout d'abord sélectionner les maillis où le drain est localisé et exécuter après la
commande Point Sources/Sinks.

Sélection des mailles

Le drain est localise dans la couche supérieure (layer 1). Nous avons besoin de sélectionner les
mailles 2 à 10 de la ligne 8. Pour cela :

1. Choisir i'outil sélection des mailles .

2. Noter que quand on déplace le curseur dans la grille, les indices i, j, k de la maille
au-dessous du curseur sont affichés dans la barre d'état en bas de l'écran comme s'est
montré dans la figure 2.

1. Sélectionner les mailles à i=8, j=2, k=1.

2. Appuyer sur la touche Shift pour permettre le mode multi-sélection et sélectionner les mailles
dans les colonnes 3-10 dans la même ligne (Figure 3).

Affectation des drains aux mailles

1. Cliquer par le bouton droit sur les mailles sélectionnées et exécuter la commande
Sources/Sinks à partir du menu flottant.
2. Cliquer sur l'onglet Drain.
3. Cliquer sur le bouton New. Ceci ajoute une nouvelle instance de drain à chacune des
mailles sélectionnée.
A ce point, nous devons entrer une côte et une conductance pour les mailles sélectionnée. Les
mailles drains ont toutes la même conductance mais la côte n'est pas la même.

Ali HAMMANI, 2021

Initiation à GMS 56

1. Entrer les valeurs suivantes our la côte et pour la conductance des drains :

ID Elevation Codnuctance

107 0 7430
108 0 7430
109 3 7430
110 6 7430
111 9 7430
112 15 7430
113 20 7430
114 27 7430
115 30 7430

1. Cliquer sur OK.

2. Désélectionner les mailles en cliquant n'importe où sur le graphique.

Le package puits
Dans la suite nous allons définir plusieurs puits en sélectionnant les mailles où les puits sont localisés
en utilisant la commande Point Sources/Sicks.

Le puits de la couche supérieure

La plus part des puits pompent dans la couche supérieure mais certains d'entre eux pompent dans
la couche du milieu ou la couche inférieure. On va tout d'abocd définir les puits dans la couche
supérieure.
1. En maintenant le bouton de la souris et la touche Shift enfoncés, sélectionner les mailles
comme c'est montré dans la figure 4.
1. Cliquer par le bouton droit dans l'une des mailles sélectionnées et exécuter la commande
Sources/Sinks du menu flottant.
2. Cliquer sur l'onglet Well.
3. Cliquer sur le bouton New.

Ali HAMMANI, 2021

Initiation à GMS 57

4. Entrer la valeur -12230 pour tous les puits (une valeur négative signifie que c'est un
pompage).
5. Cliquer sur le bouton OK.
6. Désélectionner las mailles en cliquant sur n'importe quel point de la grille.

La couche du milieu
1. Cliquer sur la flèche vers le bas dans Mini-Grid Plot.
Pour sélectionner les mailles:
1. Appuyer simultanément sur la touche Shift et la souris et sélectionner les mailles comme c'est
montré dans la figure 5.

1. Cliquer sur le bouton droit le la souris dans l'une des mailles sélectionnées et exécuter da
commande Sources/Sinks.
2. Cliquer sur l'onglet Well.
3. Cliquer sur le bouton New.
4. Entrer la valeur -12230 pour les deux mailles.
5. Cliquer sur OK.
6. Désélectionner les mailles en cliquant sur an espace en dehors de ces dent mailles.

Ali HAMMANI, 2021

Initiation à GMS 58

La couche inférieure
Finalement nous allons définir un seul puits dans la couche inférieure. Pour afficher la couche
inférieure:
1. Cliquer sur la flèche vers le bas de Mine-Grid Plot.
2. Cliquer sur la maille comme c'est montré dans la figure 6.

1. Cliquer sur le bouton droit dans la maille sélectionnée et execute la commande du menu
Sources/Sinks.
2. Cliquer sur l'onglet Well.
3. Cliquer sur New.
4. Entrer une valeur de -0.15.
5. Clique sur OK.
6. Désélectionner la maille.
A ce stade tous les puits ont été définis et on peut revenir à la première couche.
1. Clique sur le flèche vars le haut de Mini-Grid Plot.

Vérification de la Simulation
Nous avons maintenant complètement défini les données de MODFLOW et nods sommes prêts pour
lancer une simulation. Cependant, avant la sauvegarde de la simulation et l'exécution de
MODFLOW, nous devons exécuter le vérificateur du modèle MODFLOW et vérifier les erreurs. En
raison de la quantité importante de données nécessitées par las simulations de MODFLOW, il est
souvent plus facile d'oublier certaines données requises ou de trouver une inconsistance et une
incompatibilité des options ou des paramètres. De telles erreurs peuvent causer un plantage de
MODFLOW ou la production de résultats erronés. L'objectif du vérificateur du modèle (Model
Checker) est d'analyser lee données d'entrée qui sont actuellement définies pour la simulation par
MODFLOW et faire un rapport sur les problèmes potentiels. Il faut noter que l'exécution avec
succès du vérificateur du modèle ne garantit pas que la solution soit correcte. Il sert simplement à
une vérification initiale des données d'entrée pouvant économiser un temps considéreble.
Pour exécuter le vérificateur du modèle (Model Checker):
1. Exécuter la commaee du menu MODFLOW | Check Simulation.
2. Cliquer sur le bouton Run Check.

Ali HAMMANI, 2021

Initiation à GMS 59

Une liste de messages est affichée pour chaque package d'entrée de données de MODFLOW. Si
tout a été fait correctement, il n'y aura aucune erreur. Quand une erreur existe, si vous sélectionner
une erreur, GMS sélecte la maille ou la couche associé au problème.
1. Cliquer sur le bouton Done pour quitter le vérificateur de modèle.

Sauvegarde de la simulation
Maintenant on est prêt pour sauvegarder la simulation et exécuter le modèle MODFLOW.
1. Exécuter la commande du menu File | Save.
2. Choisir le répertoire.
3. Sauver le projet sous le nom gridmod.gpr.

Exécution de MODFLOW
1. Exécuter la commande du menu MOEFLOW 'Run MODFLOW'.
A ce stade MODFLOW est lancé dans une nouvelle fenêtre. Le fichier de données est passé à
MODFLOW comme une ligne de commande DOS. MODFLOW ouvre le fichier et commence la
simulation. Quand la solution est atteinte, on voit du text défiler indiquant que la simulation est en
progression.
1. Quand MOFDoOW termine la simulation cliquer sur le bouton Close.

Affichage de la solution
GMS lié automatiquement la solution quand on ferme la fenêtre MODFLOW. A ce stade on peut
voir les courbes équipotentielles des charges hydrauliques pour la couche supérieure. On peut
aussi voir que quelques mailles contiennent un triangle bleu. Ces mailles sont ``innondées'', c'est à
dire que la charge hydraulique simulée est au-dessus du toit de la couche.

Changement de la couche
Pour afficher la solution de la couche du milieu:
1. Cliquer sur la flèche vers le bas dans Mini-Grid plot.
Pour afficher la solution de la couche du milieu:
1. Cliquer sur la flèche vers le bas dans Mini-Grid Plot.
Pour retourner à la couche supérieure :
1. Cliquer deux fois our la couche vers le haut.
14.2 Remplissage den courbes de niveau
On peut afficher les courbes de niveau en utilisant l'option couleur de remplissage color fill).
1. Exécuter la Commande du menu Data | contour Options.
2. Changer Contour Method à Color Fill.
3. Clique sur OK.

Couleur de la légende
1. exécuter la commande du menu Data | Color Ramp Options.
2. Cocher l'option Legend.
3. CliqCe sur OK.

Ali HAMMANI, 2021

Initiation à GMS 60

Zone du bilan
La zone du bilan (Zone Budget) est un programme développé par l'USGS (Harbaugh 1990) qui est
utilisé pour calculer un bilan hydrique sous-régional pour MODFLoW. GMS a incorporé un modèle
de bilan similaire. Dans GMS, l'utilisateur définit des zones en affectant à leurs mailles un identifiant
(Zone Budget ID). Une fois ces zones définies, un rapport peut être généré et qui affiche le bilan
dans la zone. Le rapport inclut également une composante qui montre les entrées et les sorties
dans les zones adjacentes.

Délimitation les zones de bilan

Dans ce modèle nous allons considérer que chaque couche est une zone.

1. Cliquer sur l'outil de sélection des mailles .

2. Si nécessaire changer l'affichage en Plan View .

3. S'assurer qu'on est entrain d'afficher la deuxième couche. Ajuster la couche dans Mini-Grid
Plot sur 2.
4. Dessiner un rectangle autour de toutes les maille de la couche 2.
5. Cliquer par le bouton droit sur les mailles sélectionnées.
6. Sélectionner la commande Properties dans le menu flottant.
7. Entrer 2 pour Zone budget ID et cliquer sur OK.
8. Changer la couche à 3 dans Mini-Grid Plot.
9. Répéter les étapes de 4 à 7 sauf qu'il faut introduire 3 pour Zone budget ID.

Affichage du rapport de la zone du bilan

1. Exécuter la commande du menu Data | Flow Budget.
2. Cliquer sur l'onglet Zones.
On est entrain d'affiches le rapport de la première zone (couche supérieure). Le rapport est divisé
en deux section : écoulements à l'intérieur de la zone et écoulements vers l'extérieur de la zone.
Chaque terme puits/sources pris en compte dans le modèle est lir-sté dans le rapport. Les
échanges de flux entre zones sont également affichés.
1. Dans la Zone drop-down box choisir 2.
On peut aussi afficher le rapport de la zone 3. Après avoir fini clique sur OK.

MODFLOW – Interpolation des don-

nées d’une Couche
Pour des sites avec une stratigraphie complexe et avec un écoulement tridimensionnel, un modèle
multi-couches peut être plus précis qu'un modèle à une couche bidimensionnel. Pendant la
création d'un modèle multi-couches, la définition des données du modèle est un vrai défi. GMS
contient une suite d'outils pour interpoler et manipuler les côtes topographiques des différentes
couches. Cette section, décrit l'utilisation des points dispersés, de l'interpolation et de MODFLOW.

Modules/Interfaces requis
On aura besoin des composants suivants pour compléter le présent excercice :
• Grid ;

Ali HAMMANI, 2021

Initiation à GMS 61

• Geostatistcis ;
• Map.
Ou peut savoir si ces modules sont disponibles en exécutant la commande du menu File 'Registed'.

Mise en route
1. Si nécessaire lancer GMS. Si GMS est déjà lancé cliquer sur File | New pour assurer que la
configuration du programme soit à son état par défaut.

Interpolation dans une couche MODFLOW

Une des façons de créer les côtes du mur et du toit de chaque couche du MODFLOW est de les
interpoler à partir d'un certain nombre d'un nuage de points. L'interpolation peut être réalisée
directement en utilisant la commande du menu Interpolation | Interpolate 2D Grid lorsque on a un
nuage de points en 2D. Les points peuvent être importés à partir d'une table de valeur dans en
fichier. Pour des régions ayant une stratigraphie complexe les points doivent être définis avec plus
d'attention. L'interpolation peut être vérifiée par le vérificateur de modèle (Model Checker).

Exemple de problèmes
Pour illustrer le processus d'interpolation des côtes et de correction des erreurs, on verra une série
d'exemples de problèmes. Chaque exemple illustre un problème different et décrit une approche
simple pour une modélisation correcte de la stratigraphie.

1er cas – Couche Complexe

Le premier cas qu'on va examines est montré dans la Figure ??. C'est le cas le plus simple dans cette
section. Le site à trois couches et toutes les trois couches s'étendent sur tout le domaine du modèle.

Importation un jeu de nuage de points

La première étape dans la définition des données d'une couche est de créer une grille. L'étape
suivante est de créer une série de points à différentes localisations x,y dans le modèle. Chaque
point a une côte pour le toit et le mur pour chaque couche. Dans un problème réel ces données
peuvent provenir des forages d'exploration. Les données sont entrées dans un tableau en fichier
texte, en utilisant par exemple Excel. On importe ensuite la table de données dans un GIS. Quatre
différentes séries de nuages de points ont été préparées préalablement en utilisant la technique
décrite plus haut. Le fichier a été déjà importé dans GMS et sauver dans un fichier projet GMS. Pour
lire se projet :

1. Exécuter la commande Open .

2. Localiser et ouvrir le dossier C:\GMS\.
3. choisir le ficher nommé points.gpr puis cliquer sur Open.

Ali HAMMANI, 2021

Initiation à GMS 62

Changement en vue de face(Front View)

Avant d'interpoler, on va changer la vue de telle manière à avoir une section transversale de la
grille. De cette façon on va immédiatement voir le résultat de l'interpolation. Avant de charger de
vue, on va sélectionner une maille à l'intérieer du domaine. Si une maille est sélectionner quand on
change de vue, la nouvelle vue changera à la ligne ou à la colonne passant à travers la maille
sélectionnée.
1. Cliquer sur l'outil de sélection des mailles .
2. Sélectionner une maille près du centre de la grille.

3. Clique sur le bouton Front View .

Interpolation des valeurs de l’élevation

Le nuage de points ont chacun 4 données sur la côte topographique, top1, bot1, bot2, et bot3.
L'étape suivante consiste à interpoler chacune de ces données dans la couche appropriée dans
le maillage de MODFLOW. Avant d'interpoler ces points il faut initialiser les données de MODFLOW.
1. Exécuter la commande du menu MODFLOW | New Simulation.
2. Cliquer sur OK.
Pour interpoler les valeurs de des côtes:
1. Pour afficher les séries de données, dans Project Explorer, étendre Case 1 scatter point si
nécessaire en cliquant sur le signe +.
2. Cliquer par le bouton droit sur Cas 1 dans Project Explorer et exécuter interpolate To | MOD-
FyOW layers dans le menu flottant.
La boîte de dialogue qui apparaît est utilisée pour définir quelles données du nuage de points
sont interpolées. Les données à interpoler sont listées en haut à gauche de la boîte de dialogue,
et les variables de MODFLOW dans la partie en haut à droite. Il faut donc chercher à faire
correspondre les données aux variables en utilisant les boutons Map et Unmap. GMS essaie de faire
ces correspondances d'une manière automatique en se basant sur les noms des données. Dans le
présent cas, toutes les relations ont été correctement affectées et on peut continuer.
1. Cliquer sur OK.
Noter que les couches interpolées représentent correctement la section transversale de la Figure ??.
On peut voir l'ensemble des sections verticales en appuyant sur les cuoches Mini-Grid Display.

2ème Cas – Couche encastrée

Le cas suivant qu'on va examiner est illustré dans la figure 6-1. Dans ce cas, la couche du milieu est
n'existe que la partie gauche (ouest) du modèle. Ce cas de couche est plus difficile à modéliser
avec MODFLOW du fait un tableau entier de K doit être défini pour chaque couche.

MODFLOW – Approche du modèle

conceptual
Introduction
L'approche du modèle conceptuel consiste en l'utilisation du SIG pour développer un model con-
ceptuel sur le site à modéliser. La localisation des puits/sources, les paramètres des couches comme
la conductivité hydraulique, les limites du domaine et toute autre donnée nécessaire pour la simu-
lation peuvent être définis au niveau du modèle conceptuel. Une fois ce modèle est complété, le

Ali HAMMANI, 2021

Initiation à GMS 63

maillage est généré et le modèle conceptuel est converti à modèle par maillage et les toutes les
affectations sont faites automatique maille par maille.

Description du Problème
Le problème qu'on va résoudre est illustrer dans Figure ??.

On va modéliser l'écoulement souterrain dans les sédiments de vallée limitée par des collines au
Nord et de deux rivières au Sud. Une coupe Nord-Sud du site est représentée dans la Figure ??. Le
substratum imperméable (mur) est composé de calcaire qui affleure au Nord. Le site présente deux
couches de sédiments. La couche supérieure va être modélisée comme une nappe phréatique et la
couche inférieure sera considérée comme un aquifère captif.
La limite du Nord sera considérée à flux nul et les autres limites seront considérées à potentiel imposé
correspondant à la cote du lit des rivières. On supposera également la recharge du système se
fait par infiltration des eaux de pluie. Les ruisseaux à l'intérieur du domaine sont parfois secs mais
de temps en temps ils s'écoulent en raison du drainage de la nappe libre. On représentera ces
ruisseaux par des drains. Il existe aussi deux puits de pompage qu'on va introduire dans le modèle.

Mise en route
Si nécessaire lancer GMS. Si GMS est déjà lance, exécuter la commande du menu File | New pour
assurer que GMS restaure sa configuration initiale.

Importation de l’image de fond

La première étape dans mise en œuvre de la simulation est d'importer une image numérique du site
à modéliser. L'image a déjà été importée dans GMS et enregistrée dans un fichier de projet GMS.
Pour lire l'image on va ouvrir le fichier projet. Une fois l'image est importée dans GMS, elle peut être
affichée en arrière plan pour aidé à la conception du modèle.

Ali HAMMANI, 2021

Initiation à GMS 64

Lecture de l’image
Pour importer l'image:
1. Cliquer sur le bouton Open de la barre d'outils.
2. Localiser et ouvrir le répertoire C:\GMS.
3. Ouvrir le fichier start.gpr.
Tous les autres objets de GMS sont dessinés sur l'image. L'image apparaît uniquement en vue plan.

Enregistrer le projet
Avant de faire des changements, il faut enregistrer le projet sous un autre nom.
1. Exécuter la commande du menu File | Save As.
2. Enregistrer le projet sous le nom conceptuel.gpr.
On peut maintenant cliquer à chaque fois sur le bouton Save On peut maintenant cliquer à chaque
fois sur le bouton Save ou le raccourcis clavier ctrl s pour enregistrer le projet à chaque instant.

Définition des unités

A ce niveau, on peut aussi définir les unités utilisées dans le modèle conceptuel. Les unités qu'on
choisi vont être appliquées pour éditer les champs dans l'interface de GMS pour se souvenir des
unités propre pour chaque paramètre.
1. Exécuter la commande du menu Edit | Units.
2. Pour Length choisir m (pour mètre). Pour Time, choisir d (pour jour). On ignore les autres unités
puisqu'elles ne sont pas utilisées dans cet exercice.
3. Cliquer sur OK.

Defining the Boundary

La première étape consiste à définir la limite extérieure du domaine. On fera cela par la création
d'un arc qui va former une limite fermée autour du site.

Créer la couverture
1. Dans Project Explorer cliquer par le bouton droit dans un espace vide puis dans le menu flottant
exécuter le menu New | Conceptual Model.
2. Pour Name, entrer ExempleMC. Pour le modèle choisir MODFLOW.
3. Cliquer sur OK.

4. Clique par le bouton droit sur le modèle conceptuel ExempleMC et exécuter New Coverage
(nouvelle couverture) dans le menu flottant.
5. Changer Coverage name à Limites. Changer Default elevation à 213. Changer Default layer
range pour aller de 1 à 2.
6. Cliquer sur OK.

Ali HAMMANI, 2021

Initiation à GMS 65

Création de l’Arc

1. Cliquer sur l'outil Create Arc .

2. Commencer le tracer d'une polyligne à gauche (ouest) du domaine dans la localisation montrée
dans Figure ??.
1. Créer une polyligne en allant vet l'Est et en commençant par le Sud. Ne pas faire attention à
l'espacement ou à la localisation exacte des points. Il faut juste un nombre suffisant de points
pour définir une localisation approximative de la limite. Les limites Sud et Est coïncide avec les
rivières. La limite Nord coïncide avec la limite calcaire au Nord du Site.
2. Pour terminer la polyligne, cliquer sur le point ou on a commencé.

Construction d’une couverture local Puits/Sources

L'étape suivante consiste en la construction d'une couverture locale puits/sources. Cette couverture
définit la limite de la région modélisée ainsi que les termes puits/sources comprenant les puits, les
rivières, les drains et les charges hydrauliques générales.
Les proprieties qui peuvent être assignees aux caractéristiques des objets dans la couverture de-
pendent du modèle et des options qu'on peut configurer dans la boîte de dialogue Coverage
Setup. Avant de créer les caractéristiques des objets (feature objects), on va changer les options
dans la boîte de dialogue Coverage Setup.
1. Cliquer par le bouton droit sur la couverture Limite et choisir la commande Duplicate dans le
menu flottant. Changer le nom de la nouvelle couverture en Puits & Sources (cela peut être
fait en cliquant par le bouton droit de la souris sur le nouvelle couverture et en choisissant
Rename dans le menu flottant).
2. Cliquer par le bouton droit sur la couverture Puits & Sources et exécuter la commande
Coverage Setup dans le menu flottant.
3. Dans la liste Sources/Sinks/BCs, cocher les options qu'on va utiliser dans cette exemple.
• Layer range (Etendue de la couche).
• Wells (Puits).
• Refine points (affiner tes points).
• Specific head (charge hydraulique imposée).
• Drain.

Ali HAMMANI, 2021

Initiation à GMS 66

4. S'assurer que la case Use to define model boundary (active area) est cochée.
5. Cliquer sur OK.

Définition de la limite à potentiel imposé

L'étape suivante consiste à définir la limite à charge imposée (specific head) le lond des limites Sud
et Est du modèle. Avant cela, on doit tout d'abord couper la polyligne en trois autres polylignes.
Une polyligne va définir la limite à flux nul au Nord et les deux autres limites vont définir les deux
rivières auSud et à l'Est. Une polyligne est coupée en sélectionnant un ou plusieurs vertex dans la
plolyligne et en convertissant la les vertex en nœuds.
1. Cliquer sur l'outil Select Vertices .
2. Clique sur deux vertex comme c'est indiqué dans la Figure ??. Le vertex 1 est localisé à la
jonction entre les deux rivières. Vertex 2 est localisé en haut de la rivière Est. Pour sélectionner
les deux vertex en une seule fois, sélectionner le premier puis appuyer sur la touche Majuscule
(Shift) tout en cliquant sur l'autre vertex.

1. Cliquer par le bouton droit sur l'une des vertex et exécuter la commande Vertex --> Node dans
le menu flottant.
Maintenant que nous avons les trois polyligne, nous allons spécifier les deux polylignes des rivières
comme des limites à charge imposée.

1. Cliquer sur l'outil Select Arcs .

2. Cliquer sur les polylignes Sud et Est en sélectionnant la première polyligne et en cliquant sur la
deuxième tout en appuyant sur la touche Majuscule (Shift) du clavier.
3. Cliquer par le bouton droit sur l'une des polyligne sélectionnées puis exécuter Attribute Table
dans le menu flottant.
4. Chercher la ligne All et la colonne Type du tableau. Dans cette cellule choisir spec. Head.
Ceci va changer le type pour les deux polyligne.
5. Cliquer sur OK.
6. Clique quelque part sur l'image pour désélectionner les polylignes.
Noter que la couleur des deux polyligne a change indiquant le type de condition aux limites.
L'étape suivante consiste à définir les charges hydrauliques aux bouts des deux polylignes. La
charge le long des deux polylignes sera supposée variée linéairement le long de la polyligne.

Ali HAMMANI, 2021

Initiation à GMS 67

1. Cliquer sur l'outil Select Points/Nodes .

2. Double cliquer sur le noeud à l'ouest du domaine au bout de la limite Sud.
3. Entrer la valeur 212 pour Head-Stage.
4. Cliquer sur OK.
5. D'une manière similaire affecter la valeur 208 au noeud de la jonction des deux rivières et la
valeur 214 au neoud en haut de la limite Est.

Définition des lignes du drain

A ce stade, on va entrer les polylignes définissant les ruisseaux à l'intérieur du domaine constituant
des drains.
1. Cliquer sur l'outil Create Arc .
2. Créer la plolilygne nommée arc 1 dans Figure ??. Commencer par cliquer sur la polyligne Sud,
créer la polyligne en cliquant sur des points le long du ruisseau et terminer la polyligne en
double cliquant sur la polyligne Nord.

Noter que lorsqu'on clique sur les environs d'un vertex dans une polyligne existante, GMS fait inter-
cepter automatiquement la polyligne crée à la polyligne existante et crée un nœud à la jonction
des deux polylignes.
1. Créer les polylignes arc 1 et arc 2 comme c'est montré dans Figure ?? et de la même manière
que la polyline arc 1.
Dans la suite, on va définir les polylignes comme des drains et on va leur affecter des conductances
et des cotes topographiques.

1. Cliquer sur l'outil Select Arcs .

2. Sélectionner l'ensemble des polylignes correspondant aux drains en cliquant dessus tous en
maintenant la touche Shift enfoncé.
3. Cliquer par le bouton droit sur l'une des polylignes sélectionnées et exécuter la commande du
menu Attribute Table dans le menu flottant.
4. Dans la ligne All et la colonne Type choisir l'option drain.
5. Entrer une conductance de 555 dans la ligne All et la colonne Cond.. Ceci représente une
conductance par unité de logueur. GMS calcule automatiquement calcul automatiquement la

Ali HAMMANI, 2021

Initiation à GMS 68

conductance dans les mailles du domaine quant la les drains sont affectés à au maillage du
modèle.
6. Changer les propriétés From layer et To layer à 1 pour chaque polygne. Cela signifie sont
uniquement affectés à la couche 1 du modèle.
7. Cliquer sur OK.
Les cotes des drains sont spécifiées dans les extrémités des polylignes. Elles sont supposées variées
linéairement le long de chaque polyligne.

1. Cliquer sur l'outil Select Points/Nodes .

2. Double cliquer sur le noeud 2 in dans la Figure ??. Noter que ce nœud a 2 propriétés qui lui
sont associées puisqu'il est attacher à deux limites.

1. Entrer 212 pour Bot. elev. De la propriété drain. Ne rien changer dans la propriété spec. head.
Cliquer OK.
2. Répéter cette procédure pour affecter les cotes des noeuds du bout de chaque drain en
utilisant les valeurs de la Figure ??.

Construction des polygones

Avec le type de couverture puits/sources locale, l'ensemble de la région à modéliser doit être cou-
verte par polygones qui ne se chevauchent pas. Ceci définir la région active du maillage du
modèle. Dans la plus parts des cas, tous les polygones vont être des polygones à charge variable.
Cependant, d'autres polygones peuvent être utilisés. Par exemple, pour modéliser un la, un poly-
gone de charge général peut être utilisé. La façon la plus simple de définir les polygones est de
créer en premier lieu toutes les polylignes utilisées dans la couverture et utiliser ensuite la commande
construction de polygones (Build polygons). Cette commande cherche à travers les polylignes et
crée un polygone pour chaque boucle fermée définie par les polylignes. Ces polygones sont par
défaut de type "NONE" mai peuvent être converti à d'autres types en sélectionnant le polygone et
en utilisant la commande propriétés.
Maintenant que toutes les polylignes ont été crées dans la couverture, on est prsè pour construire
les polygons. Tous les polygones auront une charge hydraulique variable.
1. Dans le menu principal, exécuter la commande Feature Objects | Build Polygons.
Noter que le polygone n'est pas rempli. On peut changer l'affichage du polygone si l'on veut en
exécutant la commande Display | Display Options.

Ali HAMMANI, 2021

Initiation à GMS 69

Création de puits
L'étape finale dans la création de la couverture du terme Puits/Sources local est de définir les puits
de pompage. Les puits sont définis comme des types d'objet point. Deux puits vont être créés.

1. Cliquer sur l'outil Create Point .

2. Déplacer le curseur la localisation approximative Puits 1 montré dans la Figure 7 et cliquer par
la souri pour créer le point.
3. Quand le nouveau point est sélectionner introduire les coordonnées (835, 1425) dans les
champs d'édition X et Y en haut de la fenêtre GMS.

4. Cliquer sur le bouton Properties .

5. Pour Type, choisir l'option well.
6. Pour Flow rate, entrer la valeur constante -680.
7. Changer les proprieties From layer et To layer pour être égales à 1. Cela signifie que le puits
1 pompe dans la couche 1 uniquement.
8. Clique sur OK.
9. D'une manière similaire, créer l'autre puits à la localisation (3220, 1000) et lui affecter un débit
de pompage de --2830. Cependant, pour ce puits, changer From layer et To layer à 2 pour
que ce puits soit affecté à la couche 2.

Affinement du maillage
Un puits de pompage représente un point de convergence des écoulements souterrains et cause
un gradient hydraulique très raide à son voisinage. Dans le but modéliser d'une manière précise
l'écoulement près du puits, le maillage est affiné au voisinage du puits. Ce type d'affinage peut être
réalisé automatiquement par GMS en affectant des données d'affinage directement au puits dans
le modèle conceptuel.

1. Cliquer sur l'outil Select Points/Nodes .

2. Sélectionner les deux puits e cliquant sur la souris pendant que la touche Shift du clavier est
enfoncée.

3. Cliquer sur le bouton Properties .

4. Chercher la colonne Refine et dans la ligne All cocher la case. Ceci permettra au maillage
d'être affiner pour au voisinage des deux puits.
5. Changer la taille de base (Base size) à 25, Bias à 1.1 et Max size à 150 pour les deux points.
6. Cliquer sur OK.

Délimitation des zones de recharge

L'étape suivante dans la conception du modèle conceptuel est de construire une couverture qui
définie les zones de recharge. On va supposer que la recharge sur l'ensemble de la région à
modéliser est uniforme excepter pour une décharge situé à l'intérieur de cette région. La recharge
dans la surface de la décharge sera réduite en raison d'un système d'étanchéisation.

Copie de la limite
On va créer une couverture de recharge en copiant la limite :
1. Cliquer par le bouton droit de la souris sur la couverture Limite et exécuter la commande
Duplicate dans le menu folottant.

Ali HAMMANI, 2021

Initiation à GMS 70

2. Change le nom de la couverture de new coverage à Recharge.

3. Cliquer par le bouton droit sur la couverture Recharge et exécuter la commande Coverage
Setup dans le menu flottant.
4. Dans la liste Areal Propertiest, cocher la propriété Recharge rate.
5. Cliquer sur OK.

Création de la limite de la décharge

Dans la suite on va créer une polyligne délimitant la décharge. Pour cela, on va tout d'abord créer
une boucle fermée de points sous forme de rectangle dans la zone de la décharge (Figure 7). On
va ensuite éditer les nœuds et les vertex de telle façon à ce que la polyligne coïcide exactement
avec les limite de la décharge.

1. Cliquer sur l'outil Create Arc tool .

2. Crée un polygone rectangulaire dans une zone proche de la localization de la décharge
comme c'est montré dans la Figure Figure ??.

Maintenant que la polyligne est créée dans une localization approximative, on va éditer les coor-
données des vertex et des noeuds pour definer des coordonnées précises.

1. Cliquer sur l'outil Select Vertices .

2. Dessiner un rectangle autour de l'ensemble de la décharge ce qui permettra à tous lex vertex
et les nœuds d'être sélectionnés.
3. Clique par le bouton droit sur l'un des vertex sélectionné puis exécuter la commande Vertex ->
Node du menu flottant.

4. Clique sur l'outil Select Points/Nodes .

5. Cliquer sur l'un des noeuds des coins du rectangle..

Ali HAMMANI, 2021

Initiation à GMS 71

6. Pendant que le nœud est sélectionner, entrer les coordonnées exactes du nœud dans les
champs d'édition X et Y en haut de la fenêtre GMS.
7. Répéter ce processus pour les 3 nœuds restant.

Construction de polygones
Maintenant que les polylignes sont définies on peut construire les polygones.
1. Dans le menu principal exécuter la commande Feature Objects | Build Polygons.

Affectation des valeurs de la recharge

Les zones de recharge étant définies, on peut affecter les valeurs de la recharge. On va affecter
une valeur de la recharge à la décharge et une autre valeur au reste du polygone.

1. Clique sur l'outil Select Polygons tool .

2. Double cliquer sur le polygone de la décharge.
3. Changer Recharge rate 0.00006.
NB: Le taux de recharge rate est faible que celui affecté à d'autre polygone. La décharge sera
revêtue et aura un impact faible sur la recharge.
1. Clique sur OK.
2. Double cliquer sur l'autre polyligne.
3. Changer Recharge rate à 0.00695.
4. Cliquer sur OK.

Definition de la conductivité Hydraulique

Dans la suite on va entrer la conductivité hydraulique pour chaque couche. Dans beaucoup de
cas, on souhaite définir plusieurs polygones correspondant à différentes zones de conductivité hy-
draulique. Dans un souci de simplification on utilisera une seule valeur par couche.

Copie des limites

On va créer une autre couverture en copiant la couverture Limite.
1. Cliquer par le bouton droit sur la couverture Limite et exécuter la commande Duplicate du
menu dans le menu flottant.
2. Changer le nom new coverage à Couche 1.
3. Cliquer par le bouton droit sur la couverture Couche 1 et exécuter la commande Coverage
Setup dans le menu flottant.
4. Dans la liste Areal Properties, cocher les options suivantes :
• Horizontal K.
• Vertical anis..
5. Changer Default layer range de 1 à 1.
6. Cliquer OK.
7. Cliquer par le bouton droit sur la couverture Couche 1 puis exécuter la commande Duplicate
du menu flottant. Changer le nom de la couverture à Couche 2.

Ali HAMMANI, 2021

Initiation à GMS 72

8. Cliquer par le bouton droit sur la couverture Couche 2 et exécuter la commande Coverage
Setup du menu flottant.
9. Changer Default layer range de 2 à 2.
10. Cliquer OK.

Couche supérieure
Tout d'abord, on va affecter les valeurs à K pour la couche supérieure.
1. Cliquer sur la couverture Couche 1 dans Project Explorer.
2. Dans le menu principal exécuter la commande Feature Objects | Build Polygons.

3. Cliquer sur l'outil Select Polygons , double cliquer sur le polygone.

4. Changer la valeur de Horizonal K à 5.5.
5. Change la valeur de Vertical anis. à 4.
6. Cliquer OK.

Couche inférieure
Pour la couche inférieure :
1. Cliquer sur la couverture Couche 2 dans Project Explorer.
2. Dans le menu principal exécuter la commande Feature Objects | Build Polygons.

3. Cliquer sur l'outil Select Polygons , double cliquer sur le polygone.

4. Changer la valeur de Horizontal K à 10.
5. Change la valeur de Vertical anis. à 4.
6. Cliquer OK.
Ceci complète la définition des couvertures dans le modèle conceptuel. Avant de continuer la
création du maillage, on va faire de Puits & Sources la couverture active.
1. Cliquer sur la couverture Puits & Sources dans Project Explorer.

Localisation du cadre du maillage

L'ensemble des couvertures du modèle conceptuel étant complétées, on est près pour créer un le
maillage. La première étape dans la création du maillage est de définir la localisation et l'orientation
du maillage en utilisant le cadre du maillage. Celui-ci représente la limite externe du maillage. Il
peut être posisitoner au-dessus de carte géographique du site à modéliser.
1. Dans Project Explorer cliquer par le bouton droit sur un espace vide et exécuter la commande
New | Grid Frame dans le menu flottant.
2. Dans Project Explorer cliquer par le bouton droit sur Grid Frame et exécuter la commande Fit
to Active Coverage dans le menu flottant.
3. Double-cliquer sur Grid frame dans Project Explorer pour lancer la boîte de dialogue propriétés.
4. Changer Origin z: à 170 et Dimension z: à 60. Ceci fournit une série de valeurs initiales pour
les valeurs de cotes des couches du modèle MODFLOW. Plus tard, on va interpoler les valeurs
des cotes des différentes couches.
5. Clique sur OK pour sortie de la boîte de dialogue Grid Frame.

Ali HAMMANI, 2021

Initiation à GMS 73

Création de la grille
Maintenant que les couvertures et le cadre du maillage sont créés, on est près pour créer le maillage.
1. Dans le menu principale de GMS, executer la commande Feature Objects | Map[F0E0?] 3D
Grid.
Noter que le maillage est dimensionnée en utilisant les données du cadre de maillage. Si le cadre
du maillage n'existe pas, le maillage va couvrir le modèle conceptuel. De plus, le nombre de maille
dans les directions x et y ne peut pas être changé du fait que le nombre de lignes et de colonnes et
la localisation des mailles limites vont être contrôlé par les données d'affinage entrer dans les puits.
1. Dans Z-Dimension changer Num cells à 2.
2. Clique sur OK.

Initialisation des données de MODFLOW

Maintenant le maillage est construit et les zones actives/inactives sont délimitées. L'étape suivante
consistera à convertir le modèle conceptuel à modèle numérique basé sur le maillage. Avant cela,
on doit tout d'abord initialiser les données de MODFLOW :
1. Clique sur 3D Grid Data dans Project Explorer.
2. Dans le menu principal exécuter la commande MODFLOW | New Simulation.
3. Cliquer sur OK.

Defining the Active/Inactive Zones

Maintenant le maillage est créer, l'étape suivante est de definer les zones actives et inactives du
modèle. Ceci est accompli automatiquement en utilisant les informations sur la couverture des termes
Puits/Sources local.

1. Cliquer sur Map Data dans Project Explorer.

2. Cllquer sur l'outil Select Polygons .

3. Sélectionner un des polygones.
4. Cliiquer sus le bouton Properties .
5. Confirmer que l'affectation des couches est 1 à 2 et cliquer sur OK.
6. Dans le menu principal executer la commande Feature Objects | Activate cells in Coverage(s).
Chacune des mailles à l'intérieur de chaque polygone à l'intérieur de la couverture Puits/Sources est
désignée comme active et chaque maille à l'extérieur est désignée comme inactive. Noter que les
mailles sur la limite sont activées pour que la limite à flux nul au nord et la limite et la limite à potentiel
impose au Sud et à l'Est.

Interpolation des cotes des couches

Maintenant on a besoin de définir les côtes des couches et les charges initiales. Puisqu'on utilise le
package LPF, les cotes des toits et des murs sont définies pour chaque couche en considérant le
type de la couche. Pour un modèle à deux couches, on a besoin de définir les cotes du toit de la
couche 1 (altitudes du sol), le mur de la couche 1 (substratum imperméable) et le mur de la couche
2. Il est supposé que toit la couche 2 correspond au mur de la couche 1.
Pour définir les cotes, on importera des séries de points dispersés définissant les cotes et on interpolera
les cotes directement dans les mailles de chaque couche. Dans certains cas, cela est fait en utilisant
une série de points dispersés. Dans le présent cas, on va utiliser deux séries de points dispersés :

Ali HAMMANI, 2021

Initiation à GMS 74

la première pour l'altitude de la surface du sol et la seconde pour les cotes du substratum de la
couche 1 et le toit de la couche 2.

Importing the Ground Surface Scatter Points

Les points disperses ont déjà été importés puisqu'ils ont été inclus dans le projet que nous avons lu
au départ. Ces points viennent d'un fichier texte. Les séries de points sont cachés qu'on peut faire
apparaître par :
1. Dans Project Explorer, cocher les cases terrain et elevs de 2D Scatter Data. Les points appa-
raissent dans le modèle conceptuel

Interpolation des charges initiales et des cotes du terrain

1. Cliquer par le bouton droit sur le série de point terrain et exécuter la commande Interpo-
lation To | MODFLOW Layers du menu flottant.
Cette boîte de dialogue informe GMS sur quelles données de MODFLOW les points seront interpolés.
1. Cliquer sur ground_eleg et Starting Heads, et cliquer sur le Bouton bap.
2. Cliquer sur ground_elev et Top Elevations Layer 1 puis cliquer sur le button Map.
3. Cliquer sur OK pour réarliser l'interpolation.

Interpolation des cotes des couches

Pour interpoler les cotes des couches:
1. Dans Project Explorer Cliquer sur elevs sous 2D scatter Data pour le rendre actif.
2. Cliquer par le bouton droit sur elevs et exécuter Interpolation To | MODFLOW Layers dans le
menu flottant.
GMS va automatiquement faire correspondre les valeurs Bottom 1 et Bottom 2 au mur de la couche
1 et au mur de la couche 2 en se basant sur les noms des données.
1. Cliquer OK.

Ajustement de l’affichage
On peut maintenant cacher les series de points dispersés.
1. Décocher 2D scatter data dans Project Explorer.
2. Décocher Grid Frame dans Project Explorer.

Affichage d’une section transversal du modèle conceptuel

Pour vérifier l'interpolation on va afficher une section transversale.

1. Cliquer sur 3o Grid Data dans Project Explorer.

2. Clique sur une maille près du centre du maillage.

3. Cliquer sue Side View .

Pour avoir un meilleure affichage de la section, on va amplifier z.

4. Cliquer sur le bouton Display Options .

5. Entrer la valeur 5 pour Z magnification factor.
6. Cliquer OK.

Ali HAMMANI, 2021

Initiation à GMS 75

7. Cliquer le bouton Frame .

Cliquer sup les flèches de Mini Grid pour se déplacer d'une colonne à l'autre et afficher toutes les
sections.
Noter que dans la partie droite de la section transversale, le mur de la couche 2 et au- dessus de
son toot. Ceci doit être corrigé avant l'exécution du modèle.

Correction des rotes des couches

GMS fournit un certains nombre d'outils pratiques pour corriger les problèmes de cotes. Ces outils
sont inclus dans le vérificateur du modèle.
1. Dans le menu de GMS exécuter la commande MODFLOW | check Simulation.
2. Cliquer sur le bouton Run Check.
3. Cliquer sur le bouton Fix Layer Errors.
Noter que plusieurs erreurs ont été trouvées. Il existe plusieurs méthode de corriger les erreurs des
couches. On choisira l'option de tronquer le substratum (Truncate to bedrock) option. Toutes les
mailles au dessous du substratum seront considérées inactives.
1. Cliquer sur l'option Truncate to bedrotk.
2. Cliquer sur le bouton Fix Affecteted Layers.
3. Clique sur OK.
4. Cliquer sur le bouton Done.
Noter que les erreurs des couches on été corrigées. Les couches peuvent aussi être affichées en vue
de plan en cliquant sur le bouton Plan View .

Conversion du modèle Conceptual

1. Cliquer par le bouton droit sur le modèle conceptual ExempleMC et exécuter la commande
Map To | MODFLOW / MODPATH dans le menu fottant.
2. S'assurer que l'option All applicable coverages cocher puis cliquer sur OK.
Noter que les mailles interceptées par les drains, les puits et les limites à potential impose ont été
identifies par le modèle. Les charges hydrauliques et les cotes ont été déterminées par interpolation
linéaire. Les conductances des drains dans les mailles ont été déterminées en calculant la longueur
du drain qui intercepte la maille. De plus, les recharges et les conductivités hydrauliques ont été
affectées aux différentes mailles.

Vérification de la simulation
1. Cliquer sur 3D Grid Data dans Project Explorer.
2. Exécuter la commande MODFLOW | Check Simulation du menu principal.
3. Cliquer sur le bouton Run Check.
4. Cliquer sur le bouton Done.

Enregistrement du projet
1. Cliquer sur la bouton Save .

Ali HAMMANI, 2021

Initiation à GMS 76

Exécuter MODFLOW
1. Exécuter la commande MODFLOW | Run MODFLOW du menu GMS.
2. Quand la solution est obtenue cliquer sur le bouton Close.

Affichage des lignes piézométriques

1. Exécuter la commande du menu Data | Contour Options.
2. Changer et tester différentes options.

Affichage du niveau de la nappe dans une vue latérale

1. Cliquer sur l'outil Select Cell .

2. Sélectionner une maille dans le modèle.

3. Cliquer sur le bouton Side View .

4. Cliquer sur le bouton Plan View .

Viewing the Flow Budget

1. Cliquer sur Map Data dans Project Explorer.

2. Cliquer sur l'outil Select Arcs .

3. Clique sur la polyligne du drain.
4. Noter que le flux total à travers la polyligne est affiché en bas . Dans ce qui suit on affichera
l'écoulement ver la rivière.
5. Cliquer sur l'une des polyligne à charge imposée.
6. Maintener la touche Shift et sélectionner tous les polylignes de la charge imposée.
7. Noter que le flux total est affiché pour l'esemble des polylignes sélectionnées.
8. Sélectionner le dossier 3D Grid Data dans l'explorateur du projet.
9. Sélectionner un group de maille en dessinant un rectangle autour des mailles.
10. Sélectionner les données | Flow Budget command.
11. Cette boîte de dialogue affiche montre le bilan global pour les mailles sélectionnées.
12. Cliquer sur Done pour sortir de la boîte de dialogue.
13. Clique n'importe où à l'extérieur du modèle pour désélectionner les mailles.

Ali HAMMANI, 2021

UtilisationdupackagePython
Chapitre 5

Flopypourlamodélisationdes
écoulementsouterrains
Modélisation en 2D un aquifère en
charge en régime permanent
Ce Jupyter Notebook parcourra la création d'un aquifère captif monocoucheen utilisant F lopy et
M ODF LOW − 2005 pour exécuter une simulation en régime permanent de l'écoulement des eaux
souterraines dans un domaine de 100 mx 100 m avec des conditions aux limites en terme de charge
hydrauliquee de 10m d'un côté et de 0 m de l'autre côté, comme indiqué ci-dessous :

Introduction
Flopy est un package python, développé pour créer et exécuter des modèles d'écoulement des eaux
souterraines par le modèle MODFLOW. Le schéma ci-dessous illustre comment Flopy communique avec
MODFLOW :
1. Les propriétés de l'aquifère sont collectées ou estimées à partir d'un système aquifère, y compris
les dimensions, les élévations, la conductivité hydraulique, le coefficient d'emmagasinement et
l'emplacement/les caractéristiques des puits, des rivières, des lacs ou d'autres caractéristiques
hydrologiques.
2. Un modèle est créé dans un script python à l'aide du package Flopy pour configurer les
différents packages que MODFLOW utilise pour exécuter ses modèles d'écoulement des eaux
souterraines, (DIS, BAS, LPF, RIV, WEL, ETC, … ), en tant qu'objets python.
3. Flopy est utilisé pour écrire les fichiers d'entrée MODFLOW à partir de ces objets ;
4. Une fonction Flopy envoie les fichiers à un programme exécutable MODFLOW spécifié qui sera
utilisé pour résoudre le modèle d'écoulement souterraine ;
5. Le programme exécutable MODFLOW génère des données binaires dans des fichiers head et
budget ;
6. Flopy lit les données binaires des charges hydrauliques et du fichier du bilan à partir des
fichiers de sortie MODFLOW ;
7. Les données sont visualisées à l'aide des capacités de représentation graphique de Flopy.

77
Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 78

Contenu du notebook
Création de packages MODFLOW dans Flopy :
1. Le package DIS
• Discrétisation du temps et de l'espace et configuration du fichier MODFLOW DIS
2. Le package BAS
• Spécification de l'activité de la cellule et configuration du fichier MODFLOW BAS.
3. Le package LPF
• Définition des propriétés du modèle et configuration du fichier MODFLOW LPF.
4. Le package OC
• Spécification des données que l'exécutable MODFLOW enregistrera comme sortie pendant
l'exécution du modèle et configuration du package MODFLOW OC.
5. Le package Paquet PCG. -> Affectation du solveur de gradient conjugué préconditionné au
modèle d'écoulement des eaux souterraines en configurant le package MODFLOW PCG.

Écriture de fichiers et exécution du modèle :

• Écrire les fichiers d'entrée MODFLOW.

• Exécutez l'exécutable de MODFLOW.

Post-processing
6. lecture de la sortie MODFLOW.
7. Représentation graphique des résultats.
• Nous allons d'abord importer les packages appropriés pour exécuter le modèle Modflow sous
Python et afficher les données
[3]: import flopy
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mp

# Spécifique à Jupyter -- Permet d'afficher les graphiques dans le in notebook

%matplotlib inline

# formater l'affichage des coefficient de la matrice

np.set_printoptions(formatter={'float': '{: 0.5f}'.format})

Créer un objet pour le modèle

Nous définissons un objet de classe de modèle Modflow en utilisant la fonction de Flopy suivante.
Les arguments indiquent à l'objet « modèle » le nom qu'on doit lui donner `modelname' et où se
trouve le fichier MODFLOW.exe pour exécuter le modèle sur votre ordinateur.
Vous devrez modifier la variable exe_name pour diriger Flopy vers le fichier exécutable de Modflow
sur votre ordinateur.
Notez que l'espace de travail du modèle par défaut est l'emplacement de votre notebook Python,
donc tous les fichiers MODFLOW y seront créés, mais vous pouvez spécifier un espace de travail
différent si vous souhaitez que les fichiers soient créés ailleurs.

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 79

[4]: modelname = "model1"

#m = flopy.modflow.Modflow(modelname, exe_name = 'mf2005')
m = flopy.modflow.Modflow(modelname, exe_name = 'Exe/MODFLOW-NWT_64')

Le package DIS de Modflow

Maintenant, nous attachons le package DIS à notre modèle. Le package DIS spécifie la discrétisation
du modèle dans l'espace et le temps. Il contient des informations concernant : - Discrétisation spatiale
(le nombre de lignes, de couches et de colonnes, ainsi que la taille des mailles); - Discrétisation
temporelle (nombre de périodes de contraintes, modèle permanent/transitoire dans une période de
pompage; durée des périodes de pompage et nombre de pas de temps dans les périodes de
pompage); - Unité des longueurs et du temps.

Discrétisation spatiale
• Nous définissons un domaine de modèle de 100x100m , et le discrétisons en 10 lignes et 10
colonnes.
[5]: # Fixer les variables de discrétisation
Lx = 100.
Ly = 100.
ztop = 0.
zbot = -50.
nlay = 1
nrow = 10
ncol = 10
dx = Lx/ncol
dy = Ly/nrow
dz = (ztop - zbot) / nlay

Discrétisation temporelle
Pour discrétiser le temps, il faut d'abord préciser combien de périodes de pompage seront simulées
dans ce modèle. Cela se fait via la variable `nper' de Flopy. Comme nous résolvons un modèle en
régime permanent, nous ne spécifierons qu'une seule période de contrainte.
[6]: # Specifier le nombre de périodes de pompage
nper = 1

Nous créons ensuite une variable appelée « steady » sous la forme d'une liste d'indicateurs booléens
« True/False », un pour chaque période de pompage, indiquant si le solveur de différences finies
doit résoudre un modèle en régime permanent ou transitoire. True=régime permanent, False=régime
transitoire. Notre liste ne contiendra qu'un seul booléen pour la seule période de pompage en
régime permanent.
[7]: # Specifier si la péiode de pompage est en régime transitoire ou en régime permanent
steady = [True]

print("steady-state data: \n", steady)

steady-state data:
[True]
Il existe plusieurs autres variables utilisées dans le package DIS qui permettent à l'utilisateur de
spécifier les propriétés des pas de temps, les données sur la période de pompage (longueur de la
période de pompage, nombre de pas de temps par période de pompage et multiplicateur de pas
de temps). Comme ce modèle est un modèle à période unique en régime permanent, ces options ne
s'appliquent pas et nous ne les spécifierons donc pas.

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 80

Créer un objet DIS

Les discrétisations spatiales et temporelles sont utilisées pour créer l'objet dis de Flopy qui sera
ensuite utilisé pour créer le package MODFLOW DIS.
[8]: # Creer un objet Flopy de discretization, les longueurs et le temps sont␣
,→respectivement en mètre (2) and est en jour (4)

dis = flopy.modflow.ModflowDis(model=m, nlay=nlay, nrow=nrow, ncol=ncol,

delr=dx, delc=dy, top=ztop, botm=zbot,
itmuni = 4, lenuni = 2,
nper=nper, steady=steady)

Si vous souhaitez voir la grille que vous avez créée, Flopy a des capacités de représentations
graphiques pour visualiser les attributs ``modelmap'' tels que la grille, qu'elle extrait de l'objet ``dis''.
Notez que modelmap utilise Matplotlib, vous devrer donc vous assurer que ce package est également
importé en haut de votre script !
[9]: # Utiliser Flopy pour représenter graphiquement la grille de différences finies du␣
,→modèle 'm'

modelmap = flopy.plot.PlotMapView(model=m, layer=0)

grid = modelmap.plot_grid()
plt.xlabel('$L_x$ (m)',fontsize = 14)
plt.ylabel('$L_y$ (m)',fontsize = 14)
plt.title('Grille de différences finies', fontsize = 15, fontweight = 'bold')
plt.show()

## Le package BAS de MODFLOW

Ensuite, nous attachons le package BAS à notre modèle. Le package BAS spécifie le type des mailles
(comment les cellules participeront à l'exécution du modèle) ainsi que les charges hydrauliques de
départ à utiliser dans les mailles au début de la simulation de l'écoulement des eaux souterraines.

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 81

Définir le type de la maille : IBOU N D

Nous attribuons l'activité de la cellule à l'aide d'une variable appelée IBOUND. Ibound est défini
comme un tableau d'entiers de dimensions (nlay, nrow, ncol), avec la valeur 1 pour chaque maille
et avec des indicateurs pour le type de mailles : - Active (la maille a des équations aux différences
finies actives calculant la charge hydraulique pour cette maille ainsi que les flux entre ses voisines); -
Inactive (la maille n'a pas d'équations aux différences finies actives qui lui sont associées); - Charge
constante (la maille a une équation aux différences finies active pour les flux entre les cellules mais
pas pour la charge hydraulique au niveau de la maille car elle est maintenue constante tout au
long de la simulation).
Pour ce modèle, la première et la dernière colonne sont définies comme étant des mailles à charges
constantes et le reste du domaine est défini comme actif.
[10]: # Créer la variable ibound comme une matrice dont tous les coefficients sont égaux␣
,→à 1 (de type int : entier)

ibound = np.ones((nlay, nrow, ncol), dtype=np.int32)

# Affecter les valeurs -1 aux mailles des limites gauche et droite du domaine où la␣
,→charge est imposée

ibound[:, :, 0] = -1
ibound[:, :, -1] = -1

print("Les valeurs de ibound : \n", ibound)

Les valeurs de ibound :

[[[-1 1 1 1 1 1 1 1 1 -1]
[-1 1 1 1 1 1 1 1 1 -1]
[-1 1 1 1 1 1 1 1 1 -1]
[-1 1 1 1 1 1 1 1 1 -1]
[-1 1 1 1 1 1 1 1 1 -1]
[-1 1 1 1 1 1 1 1 1 -1]
[-1 1 1 1 1 1 1 1 1 -1]
[-1 1 1 1 1 1 1 1 1 -1]
[-1 1 1 1 1 1 1 1 1 -1]
[-1 1 1 1 1 1 1 1 1 -1]]]

Définir les charges hydrauliques initiales

Les valeurs des charges hydrauliques initiales sont données maille par maille avec la variable strt
définie comme un tableau de valeurs de type réell ayant les dimensions (nlay, nrow, ncol). Pour ce
modèle, les limites gauche et droite sont définies pour être des charges constantes dans la variable
IBOU N D, de sorte que les charges initiales affectées à ces mailles seront définies dans tout le
modèle. Puisque nous nous intéressons à l'écoulement en régime permanent à partir d'une limite
gauche de h = 10m et d'une limite droite de h = 0m, nous affectons ces conditions dans la variable
strt.
[11]: # Créer la variable strt comme une matrice dont les coefficients sont égaux à 1␣
,→(type floats : réel)

strt = np.ones((nlay, nrow, ncol), dtype=np.float32)

# Affecter la valeur d'une charge hydraulique égale à 10 m pour les mailles de␣
,→gauche

strt[:, :, 0] = 10.
# Affecter la valeur d'une charge hydraulique égale à 0 m pour les mailles de droite
strt[:, :, -1] = 0.

print("Valeurs des charges initiales: \n", strt)

Valeurs des charges initiales:

[[[ 10.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 82

1.00000 1.00000 0.00000]

[ 10.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 0.00000]
[ 10.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 0.00000]
[ 10.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 0.00000]
[ 10.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 0.00000]
[ 10.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 0.00000]
[ 10.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 0.00000]
[ 10.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 0.00000]
[ 10.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 0.00000]
[ 10.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 0.00000]]]

Create BAS object

We now use strt and IBOUND to make the `bas' flopy object which will later be used to create the
MODFLOW BAS package.

Créer un objet BAS

Nous utilisons maintenant la variable strt et la variable IBOU N D pour créer l'objet `bas' de Flopy
qui sera plus tard utilisé pour créer le package MODFLOW BAS.
[12]: # Créer l'objet bas de Flopy bas
bas = flopy.modflow.ModflowBas(m, ibound=ibound, strt=strt)

Si vous souhaitez visualiser la variable ibound affectée à ce package bas, M odelM ap de _Flopy-
peut également être utilisée pour afficher ibound.
[13]: # Représenter graphiquement la grille et ibound
modelmap = flopy.plot.PlotMapView(model=m, layer=0)
grid = modelmap.plot_grid()
ib = modelmap.plot_ibound()
# Ajouter la légende
plt.xlabel('$L_x$ (m)',fontsize = 14)
plt.ylabel('$L_y$ (m)',fontsize = 14)
plt.title('Ibound', fontsize = 15, fontweight = 'bold')
plt.legend(handles=[mp.patches.Patch(color='blue',label='Charge␣
,→constante',ec='black'),

mp.patches.Patch(color='white',label='Maille active',ec='black'),
mp.patches.Patch(color='black',label='Maille␣
,→inactive',ec='black')],

bbox_to_anchor=(1.6,1.0))
plt.show(modelmap)

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 83

Le package LPF de MODFLOW

Le package LPF spécifie les propriétés de la couche associées à l'écoulement dans chaque maille.
Il contient des informations concernant la conductivité hydraulique horizontale et verticale, le type
d'aquifère de la couche (en charge/à surface libre), le coefficient d'emmagasinement spécifique et
le coefficient d'ammagasinement parmi d'autres indicateurs pour l'exécution du modèle pour lequel
nous utiliserons les paramètres par défaut de Flopy.

Définir la conductivité hydraulique

Ici, nous définissons ``hk'' comme une variable pour la conductivité hydraulique horizontale de 1 m/j
qui sera la même dans les directions x et y, et ``vka'' comme une variable pour la conductivité
verticale de 1 m/j. hk et vka sont des tableaux de valurs réelles de dimensions (nlay, nrow, ncol)
avec des valeurs de conductivité pour chaque maille. Vous pouvez également spécifier un rapport
d'anisotropie dans la direction y pour la fonction LPF de Flopy, si hk doit être anisotrope dans le
plan xy.
[14]: # Définir la conductivité hydraulique horizontale
hk = np.ones((nlay,nrow,ncol), dtype=np.float32)
vk = np.ones((nlay,nrow,ncol), dtype=np.float32)

print("Valeurs de la conductivité hydraulique horizontale : \n", hk,

"\n Valeurs de la conductivité hydraulique verticale: \n", vk)

Valeurs de la conductivité hydraulique horizontale :

[[[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]
[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]
[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 84

1.00000 1.00000 1.00000]

[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]
[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]
[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]
[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]
[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]
[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]
[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]]]
Valeurs de la conductivité hydraulique verticale:
[[[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]
[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]
[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]
[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]
[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]
[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]
[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]
[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]
[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]
[ 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000]]]

Définir le coeficient d’emmagasinement spécifique

Le coeficient d'emmagasinement spécifique est défini dans un format similaire à la conductivité hy-
draulique avec une valeur de 1x10−5 m−1 pour chaque cellule.
[15]: # Déefinir le coefficient d'emmagasinnement specifique
ss = np.ones((nlay,nrow,ncol), dtype=np.float64)
ss[:,:,:] = 1e-5

print("valeurs du coefficient d'emmagasinnement specifique: \n", ss)

valeurs du coefficient d'emmagasinnement specifique:

[[[ 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001
0.00001 0.00001 0.00001]
[ 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001
0.00001 0.00001 0.00001]
[ 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001
0.00001 0.00001 0.00001]
[ 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001
0.00001 0.00001 0.00001]
[ 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 85

0.00001 0.00001 0.00001]

[ 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001
0.00001 0.00001 0.00001]
[ 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001
0.00001 0.00001 0.00001]
[ 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001
0.00001 0.00001 0.00001]
[ 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001
0.00001 0.00001 0.00001]
[ 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001
0.00001 0.00001 0.00001]]]

Définir le type d’aquifère

Le type de couches peuvent être spécifié sur flopy avec la variable laytyp, une matrice d'indicateurs
d'entiers (nlay) pour dire si la couche du modèle est :
• =0 en charge : Les couches aquifères seront traitées comme étant en charge;
• 0 Convertible : Les couches aquifères seront traitées comme étant en charge si le
niveau d'eau est au-dessus du la limite supérieur de la couche aquifère et libre si le
niveau d'eau est en dessous du la limite supérieure de la couche aquifère.
Ce modèle monocouche sera traité comme étant en charge.
[16]: # Défine le type de la couche comme t=étant en charge
laytyp = np.zeros((nlay,), dtype=np.int32)

print("Valeurs des types de couches : \n", laytyp)

Valeurs des types de couches :

[0]

Définir l’objet LPF

Nous utilisons maintenant hk, vk, ss et laytyp pour créer l'objet `lpf ' de Flopy qui sera utilisé pour créer
le package MODFLOW LPF. ipakcb est un indicateur pour dire si nous allons enregistrer les données
de bilan de la maille (flux, etc.) (ipakcb > 0 = enregistrer les données). (Noter : la documentation
Flopy affirme que les données sur le bilan sont enregistrées en utilisant les paramètres par défaut de
ipakcb, mais il semble que ipakcb doit être explicitement spécifié pour être enregistré.)
[17]: lpf = flopy.modflow.ModflowLpf(model=m, hk=hk, vka=vk, ss=ss, laytyp=laytyp,␣
,→ipakcb=1)

Le Package OC de MODFLOW
Le package OC indique à MODFLOW quels attributs des solutions de la charge hydraulique et de
débit doivent être affichés et enregistrés pour chaque période de pompage pendant l'exécution
du modèle.

Créer des données de période de pompage pour le contrôle des sorties du modèle
Pour désigner les données à sauvegarder, nous créons la variable spd comme dictionnaire de don-
nées de la période de pompage à utiliser comme argument pour la fonction Flopy oc. Cela prend
la forme :
spd = {(période de stress, pas de temps): ['IMPRIMER TETE', 'IMPRIMER BUDGET', 'ENREGISTRER TETE'

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 86

Notez que la liste des chaînes de texte peut être modifiée si vous ne souhaitez enregistrer qu'une
partie des données. Ici, nous choisirons d'enregistrer et de sauvegarder les valeurs du bilan et de la
charge pour notre modèle en régime permanent avec une seule période de pomage qui a un seul
pas de temps, comme nous l'avons spécifié dans l'objet « dis ».
[18]: # Créer les données oc de la période de pompage.
spd = {(0, 0): ['print head', 'print budget', 'save head', 'save budget']}

print("données oc de la période de pompage: \n", spd)

données oc de la période de pompage:

{(0, 0): ['print head', 'print budget', 'save head', 'save budget']}

Définir l’objet OC
Nous utiliserons spd pour créer l'objet oc de _Flopy- qui sera ensuite utilisé pour créer le package
MODFLOW OC.
[19]: oc = flopy.modflow.ModflowOc(model=m, stress_period_data=spd, compact=True)

Le packace PCG de MODFLOW

L'attachement du package PCG à un modèle MODFLOW indique à MODFLOW d'utiliser le package
de gradient conjugué préconditionné pour résoudre l'équation aux différences finies. Étant donné
que ce modèle est en régime permanent dans un domaine très simple, le package P CG devrait avoir
peu de difficulté à résoudre l'équation d'écoulement parabolique des eaux souterraines avec ses
paramètres par défaut. Cependant, MODFLOW dispose de nombreux autres packages de solveurs
(Newton-Raphson (NWT), Direct Solver (DE4), Strongly Implicit Procedure (SIP), etc.) qui peuvent être
applicables à d'autres modèles. De plus, vous souhaiterez peut-être modifier le critère de résolution
en fonction de votre modèle.

Define PCG Object

Here we will create the flopy pcg object, with its default solver settings, which will be used to create
the MODFLOW PCG package.

Définir l’objet PCG

Ici, nous allons créer l'objet pcg de _Flopy$, avec ses paramètres de solveur par défaut, qui sera
utilisé pour créer le package MODFLOW PCG.
[31]: pcg = flopy.modflow.ModflowPcg(model=m)

Écrire les fichiers d’entrée MODFLOW

Maintenant que nous avons créé tous les objets de Flopy pour les packages que nous souhaitons
attacher à notre modèle (dis, bas, lpf , oc et pcg), nous pouvons utiliser Flopy pour écrire les fichiers
d'entrée de MODFLOW à lire par l'exécutable MODFLOW. Cela se fait simplement avec :
[32]: # Ecrire les fichiers d'entrée de MODFLOW
m.write_input()

Flopy devrait avoir créé 6 nouveaux fichiers dans votre répertoire de travail qui ressemblent à
quelque chose comme :
Le fichier auquel faire attention est le fichier ``my_model.nam''. Ce fichier détermine quels packages
seront joints à votre fichier MODFLOW. À l'intérieur du nom du fichier, il répertoriera tous les packages

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 87

à attacher à votre modèle et spécifiera les fichiers BINARY DATA qu'il écrira en tant que sortie.
(mon_modèle.hds, mon_modèle.cbc)
Un autre fichier auquel faire attention, après l'exécution, est le fichier LIST. Celui là contiendra des
informations supplémentaires sur la façon dont les packages sont attachés à votre modèle ainsi que
les bilans massiques de MODFLOW pendant l'exécution du modèle.

Exécuter le modèle MODFLOW

Nous pouvons maintenant exécuter le modèle MODFLOW avec les packages spécifiés dans le fichier
.nam en utilisant l'exécutable spécifié dans la définition de notre modèle de disquette `m' :
[34]: # Exécuter le modèle
success, mfoutput = m.run_model(pause=False, report=True)
if not success:
raise Exception("MODFLOW ne s'est pas terminé normalement.")

FloPy is using the following executable to run the model: Exe/MODFLOW-NWT_64.exe

MODFLOW-NWT-SWR1
U.S. GEOLOGICAL SURVEY MODULAR FINITE-DIFFERENCE GROUNDWATER-FLOW MODEL
WITH NEWTON FORMULATION
Version 1.1.4 4/01/2018
BASED ON MODFLOW-2005 Version 1.12.0 02/03/2017

SWR1 Version 1.04.0 09/15/2016

Using NAME file: model1.nam

Run start date and time (yyyy/mm/dd hh:mm:ss): 2021/12/23 12:42:32

Solving: Stress period: 1 Time step: 1 Groundwater-Flow Eqn.

Run end date and time (yyyy/mm/dd hh:mm:ss): 2021/12/23 12:42:32
Elapsed run time: 0.010 Seconds

Normal termination of simulation

Si l'exécution a réussi, vous devriez avoir vu une fenêtre rapportant les données d'exécution.

Lire les données de charge/débit

Nous allons maintenant utiliser le lecteur de fichiers binaires et les capacités de de représentation
graphique de Flopy pour lire et afficher les résultats du modèle.

Lire les résultats de la charge hydraulique simulée

Les drésultats de la charge hydraulique simulée du fichier model1.hds seront stockées en tant qu'objet
de tête de disquette.
[36]: # Extraire les données binaire du fichier des charges hydrauliques en tant qu'objet␣
,→headobj de Flopy

headobj = flopy.utils.binaryfile.HeadFile(modelname+'.hds')

print("Objet charges hydrauliques de Flopy: \n", headobj)

Objet charges hydrauliques de Flopy:

<flopy.utils.binaryfile.HeadFile object at 0x000002B99459DA00>

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 88

Cet objet de charge hydraulique a différents attributs qui peuvent être extraits de celui-ci, tels que
les temps utilisés dans l'exécution du modèle, les données sue les charges hydrauliques, les octets de
données, etc.
Pour ce modèle, nous allons extraire les données des charges hydrauliques pour l'unique période de
pompage (totim = 1, 0).
[24]: # Extraire les données des charges hydraulique de l'objet headobj
head = headobj.get_data(totim=1.0)

Lire les données de flux

Le lecteur de données binaires de Flopy sur les bilans sera utilisé pour lire les données des bilans du
fichier model1.cbc en tant qu'objet budget de Flopy.
[37]: # Extraire les données binaire du fichier des bilans comme un objet budget de Flopy
budgobj = flopy.utils.binaryfile.CellBudgetFile(modelname+'.cbc')

print("flopy budget object: \n", budgobj)

flopy budget object:

<flopy.utils.binaryfile.CellBudgetFile object at 0x000002B99459D790>
We can now extract the data from this by passing a string in the CellBudgetFile function's ``text''
argument to call which part of the data we want to get. Here we will access the flow data for the
right face and front face of each grid cell and save the data as variables `frf' and `fff' which are
output as arrays.
Nous pouvons maintenant en extraire les données en passant une chaîne de caractères dans
l'argument ``text'' de la fonction CellBudgetFile pour appeler la partie des données que nous voulons
obtenir. Ici, nous accéderons aux données de flux pour la face droite et la face avant de chaque
maille de la grille et enregistrerons les données en tant que variables « f rf » et « f f f » qui sont
sorties sous forme de tableaux.
[42]: frf = budgobj.get_data(text='flow right face', totim=1.0)
fff = budgobj.get_data(text='flow front face', totim=1.0)

print("Flux à travers la face droite de la maille (m^3/d) \n", frf,

"\n Flux à travers la face avant de la maille m^3/d \n", fff)

Flux à travers la face droite de la maille (m^3/d)

[array([[[ 55.55556, 55.55556, 55.55556, 55.55556, 55.55556,
55.55556, 55.55556, 55.55556, 55.55556, 0.00000],
[ 55.55556, 55.55556, 55.55556, 55.55556, 55.55556,
55.55555, 55.55556, 55.55556, 55.55556, 0.00000],
[ 55.55556, 55.55556, 55.55556, 55.55556, 55.55556,
55.55555, 55.55556, 55.55556, 55.55556, 0.00000],
[ 55.55556, 55.55556, 55.55556, 55.55556, 55.55556,
55.55556, 55.55556, 55.55556, 55.55556, 0.00000],
[ 55.55556, 55.55556, 55.55556, 55.55556, 55.55556,
55.55556, 55.55556, 55.55556, 55.55556, 0.00000],
[ 55.55556, 55.55556, 55.55556, 55.55556, 55.55556,
55.55556, 55.55556, 55.55556, 55.55556, 0.00000],
[ 55.55556, 55.55556, 55.55556, 55.55556, 55.55556,
55.55556, 55.55556, 55.55556, 55.55556, 0.00000],
[ 55.55556, 55.55556, 55.55556, 55.55556, 55.55556,
55.55556, 55.55556, 55.55556, 55.55556, 0.00000],
[ 55.55556, 55.55556, 55.55556, 55.55556, 55.55556,
55.55556, 55.55556, 55.55556, 55.55556, 0.00000],
[ 55.55556, 55.55556, 55.55556, 55.55556, 55.55556,

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 89

55.55556, 55.55556, 55.55556, 55.55556, 0.00000]]],

dtype=float32)]
Flux à travers la face avant de la maille m^3/d
[array([[[ 0.00000, -0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
-0.00000, -0.00000, 0.00000, 0.00000],
[ 0.00000, -0.00000, -0.00000, -0.00000, 0.00000, 0.00000,
0.00000, 0.00000, -0.00000, 0.00000],
[ 0.00000, -0.00000, 0.00000, 0.00000, -0.00000, -0.00000,
0.00000, 0.00000, 0.00000, 0.00000],
[ 0.00000, -0.00000, -0.00000, -0.00000, -0.00000, -0.00000,
0.00000, 0.00000, 0.00000, 0.00000],
[ 0.00000, 0.00000, 0.00000, 0.00000, -0.00000, 0.00000,
-0.00000, 0.00000, -0.00000, 0.00000],
[ 0.00000, 0.00000, -0.00000, -0.00000, -0.00000, -0.00000,
0.00000, -0.00000, -0.00000, 0.00000],
[ 0.00000, -0.00000, -0.00000, -0.00000, -0.00000, -0.00000,
0.00000, -0.00000, 0.00000, 0.00000],
[ 0.00000, -0.00000, -0.00000, 0.00000, 0.00000, 0.00000,
0.00000, -0.00000, -0.00000, 0.00000],
[ 0.00000, -0.00000, 0.00000, -0.00000, -0.00000, 0.00000,
0.00000, 0.00000, 0.00000, 0.00000],
[ 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
0.00000, 0.00000, 0.00000, 0.00000]]], dtype=float32)]
Ici, il est bon de vérifier si les résultats semblent raisonnables ; pour ce modèle en régime permanent, il
est logique que nous voyions une valeur constante dans le tableau de sortie ``Flow Right Face'' car
la charge la plus élevée se trouve sur le côté gauche du modèle induit un flux de gauche à droite.
Le tableau ``Flow Front Face'' contient de très petites valeurs positives et négatives, suggérant qu'il
existe des oscillations numériques potentielles dans la solution. Cependant, ces fluctuations sont si
faibles que la sortie peut toujours être appropriée en fonction de l'application/de la précision que
nous recherchons.
Il est bon de rechercher des nuances numériques dans la pratique de la modélisation, mais les don-
nées sont beaucoup plus faciles à visualiser en exploitant les capacités de représentation graphique
de Flopy.

Représentation graphique
[43]: # Représenter graphiquement les résultats
plt.figure(figsize=(10,10)) # Créer une figure de taille 10 x 10
modelmap = flopy.plot.PlotMapView(model=m, layer=0) # Utiliser modelmap pour␣
,→attache le graphique au modèle

grid = modelmap.plot_grid() # Représenter graphiquement la grille

contour_levels = np.linspace(head[0].min(),head[0].max(),11) # Configurer la␣
,→représentation des lignes izopièzes

head_contours = modelmap.contour_array(head, levels=contour_levels) # Créer les␣

,→lignes izopièzes

flows = modelmap.plot_discharge(frf[0], fff[0], head=head) # Créer les flèches de␣

,→flux

# Afficher les paramètres

plt.xlabel('$L_x$ (m)',fontsize = 14)
plt.ylabel('$L_y$ (m)',fontsize = 14)
plt.title('Modèle de régime permanent, Flux $(m^3/j)$ et charge (m) ', fontsize =␣
,→15, fontweight = 'bold')

plt.colorbar(head_contours,aspect=5).set_label('Charge hydraulique (m)', fontsize =␣

,→15)

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 90

plt.show(modelmap)

C:\Users\aliha\miniconda3\lib\site-packages\flopy\plot\map.py:819:
DeprecationWarning: plot_discharge() has been deprecated and will be replaced in
version 3.3.5. Use plot_vector() instead, which should follow after
postprocessing.get_specific_discharge()
warnings.warn(
C:\Users\aliha\miniconda3\lib\site-packages\flopy\plot\plotutil.py:1630:
DeprecationWarning: centered_specific_discharge() has been deprecated and will
be removed in version 3.3.5. Use postprocessing.get_specific_discharge()
instead.
warnings.warn(

3D Surface Plot
Note: This is not part of flopy's plotting, and uses a matplotlib function from a 3d projection toolkit.
Flopy takes cell indexing where the top left of the grid is the (0,0) cell index while the Length & Width
units start at 0 in the lower left of the grid. This seems a little wierd, but it makes it easy to match up
an array of cell by cell numbers to their respective locations on the grid when assigning properties
or observing output data. Flopy's plotting functions automatically flip the resultant head data to
display on its proper grid, however, for the 3d plots below, you'll see the function np.flipud() is used
to flip the data array to plot in the same direction.

Tracé de surface 3D
Remarque : Cela ne fait pas partie de fonctionalités de représentation graphique de Flopy et
utilise une fonction Matplotlib d'une boîte à outils de projection 3D. Flopy prend l'indexation des

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 91

mailles où le coin supérieur gauche de la grille est l'index de cellule (0, 0) tandis que les unités de
longueur et de largeur commencent à 0 dans le coin inférieur gauche de la grille. Cela semble
un peu étrange, mais cela permet de faire facilement correspondre un tableau de numéros de
maille par maille à leurs emplacements respectifs sur la grille lors de l'attribution de propriétés ou de
l'observation des données de sortie. Les fonctions de représentation graphique de Flopy retournent
automatiquement les données des charges hydrauliques résultantes pour les afficher sur sa grille
appropriée, cependant, pour les graphiques 3D ci-dessous, vous verrez que la fonction np.flipud()
est utilisée pour retourner le tableau de données afin de tracer dans la même direction.
[50]: # Importer toolkit des axes 3d de matplotlib
from mpl_toolkits.mplot3d import Axes3D

# Créer une figure 3D

fig_3d = plt.figure(figsize=(12,15))
ax = fig_3d.gca(projection='3d')

# Affecter X, Y, Z variables pour la représentation graphique 3d pour être le␣

,→domaine du modèle et la solution des charges hydrauliques

X = np.arange(0,Lx,dx)
Y = np.arange(0,Ly,dy)
X, Y = np.meshgrid(X, Y)
Z = np.flipud(head[0])

# Créer la surface et légende

surf = ax.plot_surface(X,Y,Z, cmap = plt.cm.coolwarm, linewidth=0,␣
,→antialiased=False, label='head')

fig_3d.colorbar(surf,shrink=0.5,aspect=5).set_label('Head␣
,→(m)',fontsize=10,fontweight='bold')

ax.set_xlabel('Lx (m)', fontsize=15, fontweight='bold')

ax.set_ylabel('Ly (m)', fontsize=15, fontweight='bold')
ax.set_title('Steady-State Model Head Profile', fontsize=15, fontweight='bold')
plt.show(surf)

C:\Users\Public\Documents\Wondershare\CreatorTemp/ipykernel_25112/3839266141.py:
6: MatplotlibDeprecationWarning: Calling gca() with keyword arguments was
deprecated in Matplotlib 3.4. Starting two minor releases later, gca() will take
no keyword arguments. The gca() function should only be used to get the current
axes, or if no axes exist, create new axes with default keyword arguments. To
create a new axes with non-default arguments, use plt.axes() or plt.subplot().
ax = fig_3d.gca(projection='3d')

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 92

Congratulations! You just used flopy to run a single-layer, steady-state model in MODFLOW-2005!
Activity: Try going back to where we defined the starting heads in the LPF package and changing
the left and right boundary conditions. How does the head surface change?
What happens if you add another row of constant head cells to the model by modifying the IBOUND
variable and set those to a different value of starting heads?
Félicitations! Vous venez d'utiliser Flopy pour exécuter un modèle monocouche en régime permanent
dans MODFLOW !
Activity : Essayez de revenir à l'endroit où nous avons défini les charges initiales dans le package
LPF et de modifier les conditions aux limites gauche et droite. Comment la surface de la tête change-
t-elle?
Que se passe-t-il si vous ajoutez une autre ligne de maille à charge constante au modèle en
modifiant la variable IBOU N D et en les définissant sur une valeur différente de la charge intiale ?
[ ]:

Flopy 2 : Puits de pompage en régime

permanent
Ce script parcourra la création d'un modèle d'un aquifère captif monocouche avec Flopy et utilisera
MODFLOW pour exécuter une simulation en régime permanent de l'écoulement des eaux souterraines
dans un domaine de 100 m x 100 m avec des conditions aux limites de charge hydraulique de

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 93

10 m surl'ensemble du pourtour du domaine et un puits de pompage au centre du domaine, comme

indiqué ci-dessous :
Ce notebook se concentre sur le package MODFLOW WEL et s'appuiera sur le modèle monocouche
en régime permanent créé dans le notebook ``Flopy 1'' afin que des aspects similaires soient créés
très rapidement. Si vous ne comprenez pas une partie du code fait, n'hésitez pas à revenir à ce
notebook pour revoir ce qui se passe.

Importer des packages et créer un objet pour le modèle

[1]: #import des packages
import flopy
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mp

%matplotlib inline

Remarque : Vous devrez modifier la variable exen ame pour diriger Flopy vers le fichier exécutable
modflow sur votre ordinateur !
[2]: # Créer l'objet du modèle
modelname = "model2"
#m = flopy.modflow.Modflow(modelname, exe_name = 'mf2005')
m = flopy.modflow.Modflow(modelname, exe_name = 'Exe/MODFLOW-NWT_64')

Discrétiser le modèle et joindre le package DIS

[5]: # Fixer les variables de discrétisation
Lx = 100.
Ly = 100.
ztop = 0.
zbot = -50.
nlay = 1
nrow = 25
ncol = 25
dx = Lx/ncol
dy = Ly/nrow
dz = (ztop - zbot) / nlay

# Specifier le nombre de périodes de pompage

nper = 1

# Specifier si la péiode de pompage est en régime transitoire ou en régime permanent

steady = [True]

# Creer un objet Flopy de discretization, les longueurs et le temps sont␣

,→respectivement en mètre (2) and est en jour (4)

dis = flopy.modflow.ModflowDis(model=m, nlay=nlay, nrow=nrow, ncol=ncol,

delr=dx, delc=dy, top=ztop, botm=zbot,
itmuni = 4, lenuni = 2,
nper=nper, steady=steady)

[8]: # Vérifie la grille

modelmap = flopy.plot.PlotMapView(model=m, layer=0)
grid = modelmap.plot_grid()

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 94

plt.xlabel('$L_x$ (m)',fontsize = 14)

plt.ylabel('$L_y$ (m)',fontsize = 14)
plt.title('Grille de différences finies', fontsize = 15, fontweight = 'bold')
plt.show()

Attribuer le type de maille/les charges initiales et attacher le package BAS

[10]: # Créer la variable ibound comme une matrice dont tous les coefficients sont égaux␣
,→à 1 (de type int : entier)

ibound = np.ones((nlay, nrow, ncol), dtype=np.int32)

# Affecter les valeurs -1 aux mailles des limites gauche et droite du domaine où la␣
,→charge est imposée

ibound[:, :, 0] = -1
ibound[:, :, -1] = -1
ibound[:,0,:] = -1
ibound[:,-1,:] = -1

# Affecter la valeur 10m aux éléments de la matrice strt

strt = np.ones((nlay, nrow, ncol), dtype=np.float32)
strt[:, :, :] = 10.

# Créer l'object Flopy bas

bas = flopy.modflow.ModflowBas(m, ibound=ibound, strt=strt)

[14]: # VERIFIER IBOUND

# Utiliser Flopy pour représenter la grille et ibound
modelmap = flopy.plot.PlotMapView(model=m, layer=0)
grid = modelmap.plot_grid()

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 95

ib = modelmap.plot_ibound()
# Ajouter les étiquettes et la légende
plt.xlabel('$L_x$ (m)',fontsize = 14)
plt.ylabel('$L_y$ (m)',fontsize = 14)
plt.title('Ibound', fontsize = 15, fontweight = 'bold')
plt.legend(handles=[mp.patches.Patch(color='blue',label='Charge␣
,→constante',ec='black'),

mp.patches.Patch(color='white',label='Maille active␣
,→Cell',ec='black'),

mp.patches.Patch(color='black',label='Maille␣
,→inactive',ec='black')],

bbox_to_anchor=(1.6,1.0))
plt.show(modelmap)

Définir les propriétés de la couche et attacher le package LPF

[15]: # Définir la conductivité hydraulique horizontale
hk = np.ones((nlay,nrow,ncol), dtype=np.float32)
vk = np.ones((nlay,nrow,ncol), dtype=np.float32)

# Définir le coefficient d'emmagasinement specifique

ss = np.ones((nlay,nrow,ncol), dtype=np.float64)
ss[:,:,:] = 1e-5

# Définir la couche comme captive

laytyp = np.zeros((nlay,), dtype=np.int32)

# Créer l'objet flopy 'layer property flow'

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 96

lpf = flopy.modflow.ModflowLpf(model=m, hk=hk, vka=vk, ss=ss, laytyp=laytyp,␣

,→ipakcb=1)

# V. Create Wells and Attach MODFLOW WEL Package

The MODFLOW WEL Package simulates a volumetric flux into or out of a specified model cell in units
L3
T . Active wells are specified per stress period and pump/inject at a constant rate throughout that
period.
In flopy, wells are specified in a dictionary of stress period data where the stress period key points
to a list with a list of data for each active well. This sounds a little confusing, but takes the form:
wel_spd = {stress period: [[lay, row, col, flux], [lay, row, col, flux]….]}
See https://modflowpy.github.io/flopydoc/mfwel.html for more documentation.
For this model, we will create a single well at the center of the domain:
[8]: #Create Single Well at center of domain with [lay, row, col, flux] list
pumping_rate = -100 #in m^3/d, negative for pumping/positive for injection
well_1 = [0,ncol/2,nrow/2,pumping_rate]

print("Well 1 [layer, row, column, flux]: \n", well_1)

Well 1 [layer, row, column, flux]:

[0, 12.5, 12.5, -100]
Note that the row and columns are not integer values but there is no such thing as row ``12.5''.
When implementing your well, flopy will therefore round the layer/rows/columns to the nearest integer
layer/row/column and put the well there. Hence, the well in this example will be placed in cell (0,13,13).
(We will check this by using flopy to plot the well locations after attaching the well package below.)
Additionally, flopy's cell indexing starts at (0,0) at the top left corner of the grid whereas the unit
axes (Lx, Ly) start at the bottom left. This is a little wierd, but makes it easy to match up an array of
cell by cell numbers to their respective locations on the grid when assigning properties or observing
output data. This is just something to keep in mind when assigning a new well's cell location.
We now create a dictionary with a list of all the wells active during our single stress period model.
Since we only have 1 well, this list will just contain the information for well 1. (Remember, python
indexing starts at zero, so the key 0 corresponds to stress period 1!)
[9]: #Create Dictionary With Stress Period Data
wel_spd = {0: [well_1]}

print("Well Stress Period Data: \n", wel_spd)

Well Stress Period Data:

{0: [[0, 12.5, 12.5, -100]]}
Now that we've created our stress period data, we use wel_spd to make the `wel' flopy object which
will later be used to attach the MODFLOW WEL package to our model.
[10]: #Create flopy wel object
wel = flopy.modflow.ModflowWel(model=m, stress_period_data=wel_spd)

To make sure our well is where we want it, we can visualize its location using the plot_bc() function
of flopy ModelMap and specifying the file type we wish to plot as the well file. (ftype=`WEL')
[11]: #CHECK WELL LOCATION
#use flopy to plot grid, ibound, and wells
modelmap = flopy.plot.PlotMapView(model=m, layer=0)

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 97

grid = modelmap.plot_grid()
ib = modelmap.plot_ibound()
wel = modelmap.plot_bc(ftype='WEL')
#add labels and legend
plt.xlabel('Lx (m)',fontsize = 14)
plt.ylabel('Ly (m)',fontsize = 14)
plt.title('Well and Boundary Conditions', fontsize = 15, fontweight = 'bold')
plt.legend(handles=[mp.patches.Patch(color='red',label='Well',ec='black'),
mp.patches.Patch(color='blue',label='Const Head',ec='black'),
mp.patches.Patch(color='white',label='Active Cell',ec='black'),
mp.patches.Patch(color='black',label='Inactive␣
,→Cell',ec='black')],

bbox_to_anchor=(1.5,1.0))
plt.show(modelmap)

### VI. Specify Data Output and Attach OC Package

[12]: #create oc stress period data.
spd = {(0, 0): ['print head', 'print budget', 'save head', 'save budget']}

#create flopy output control object

oc = flopy.modflow.ModflowOc(model=m, stress_period_data=spd, compact=True)

### VII. Assign PCG as Finite Difference Solver and Attach PCG Package
[13]: #assign groundwater flow solver
pcg = flopy.modflow.ModflowPcg(model=m)

### VIII. Create MODFLOW files and Run Model

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 98

[14]: #write MODFLOW input files

m.write_input()

[15]: # Run the model

success, mfoutput = m.run_model(pause=False, report=True)
if not success:
raise Exception('MODFLOW did not terminate normally.')

FloPy is using the following executable to run the model: .\MODFLOW-NWT_64.EXE

MODFLOW-NWT-SWR1
U.S. GEOLOGICAL SURVEY MODULAR FINITE-DIFFERENCE GROUNDWATER-FLOW MODEL
WITH NEWTON FORMULATION
Version 1.1.4 4/01/2018
BASED ON MODFLOW-2005 Version 1.12.0 02/03/2017

SWR1 Version 1.04.0 09/15/2016

Using NAME file: my_pumping_model.nam

Run start date and time (yyyy/mm/dd hh:mm:ss): 2021/10/17 9:03:24

Solving: Stress period: 1 Time step: 1 Groundwater-Flow Eqn.

Run end date and time (yyyy/mm/dd hh:mm:ss): 2021/10/17 9:03:24
Elapsed run time: 0.009 Seconds

Normal termination of simulation

### IX. Extract Head and Flow Data
[16]: #extract binary data from head file
headobj = flopy.utils.binaryfile.HeadFile(modelname+'.hds')
head = headobj.get_data(totim=1.0)

[17]: #extract binary data from budget file

budgobj = flopy.utils.binaryfile.CellBudgetFile(modelname+'.cbc')
frf = budgobj.get_data(text='flow right face', totim=1.0)
fff = budgobj.get_data(text='flow front face', totim=1.0)

## X. Plotting

Head Contour and Flow Plot

[18]: #plot results

plt.figure(figsize=(10,10)) #create 10 x 10 figure
modelmap = flopy.plot.PlotMapView(model=m, layer=0) #use modelmap to attach plot to␣
,→model

grid = modelmap.plot_grid() #plot model grid

contour_levels = np.linspace(head[0].min(),head[0].max(),15) #set contour levels␣
,→for contouring head

head_contours = modelmap.contour_array(head, levels=contour_levels) #create head␣

,→contours

flows = modelmap.plot_discharge(frf[0], fff[0], head=head) #create discharge arrows

#display parameters
plt.xlabel('Lx (m)',fontsize = 14)
plt.ylabel('Ly (m)',fontsize = 14)
plt.title('Steady-State Pumping, Flow(m^3/d) and Head(m) Results', fontsize = 15,␣
,→fontweight = 'bold')

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 99

plt.colorbar(head_contours, aspect=5)
plt.show(modelmap)

Head Contour Plot

The aquifer head may be more easily visualized with a filled contour plot using matplotlib's contourf()
function.
Note: Flopy takes cell indexing where the top left of the grid is the (0,0) cell index while the Length &
Width units start at 0 in the lower left of the grid. This seems a little wierd, but it makes it easy to match
up an array of cell by cell numbers to their respective locations on the grid when assigning properties
or observing output data. Flopy's plotting capabilities automatically flip the resultant head data to
display on its proper grid, however, for the rest of the plots below, you'll see the function np.flipud()
is used to flip the data array to plot in the same direction.
[23]: #create plot
plt.figure(figsize=(9,7))

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 100

cf = plt.contourf(np.flipud(head[0,:,:]),levels=contour_levels)
#display parameters
plt.xlabel('Lx (m)',fontsize = 14)
plt.ylabel('Ly (m)',fontsize = 14)
plt.title('Steady-State Pumping, Head(m) Results', fontsize = 15, fontweight =␣
,→'bold')

plt.colorbar(cf,aspect=10)
plt.show()

In this simulation, it's interesting to note the effect of the boundary conditions on the head profile
near the edges of the model, as well as the effect of the grid on the head profile near the well.
The contour lines are slightly squared at the points closest to the corners of the grid as there is a
greater distance to a constant head cell and hence there may be more draw down at these points.
Additionally, the head contours are distorted to be sort of diamond shaped at the point of pumping
as the coarse grid size compromises the ability to resolve greater head changes at this point.
Activity: Rerun the model with a refined grid spacing (increase nrow and ncol). How does this change
the shape of the head contours near the point of pumping?
Advanced Activity: Write a code to create an ibound variable that sets constant head cells in a
more ``circular'' boundary around the well. What happens to the head contours near the boundaries
now?

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 101

3D Head Surface Plot

[24]: #import 3d axes toolkit from matplotlib

from mpl_toolkits.mplot3d import Axes3D

#create 3d figure
fig_3d = plt.figure(figsize=(12,5))
ax = fig_3d.gca(projection='3d')

#set X, Y, Z variables for 3d plot to be our model domain and head solution
X = np.arange(0,Lx,dx)
Y = np.arange(0,Ly,dy)
X, Y = np.meshgrid(X, Y)
Z = np.flipud(head[0])

#create surface and labels

surf = ax.plot_surface(X,Y,Z, cmap = plt.cm.coolwarm, linewidth=0,␣
,→antialiased=False, label='head')

fig_3d.colorbar(surf,shrink=0.5,aspect=5).set_label('Head␣
,→(m)',fontsize=10,fontweight='bold')

ax.set_xlabel('Lx (m)', fontsize=15, fontweight='bold')

ax.set_ylabel('Ly (m)', fontsize=15, fontweight='bold')
ax.set_title('Steady-State Pumping, Head Surface', fontsize=15, fontweight='bold')
plt.show(surf)

Head Transect Plot

It may also be helpful to visualize the head profile along a single row. Below the head head data
is extracted halfway up Ly and plotted across Lx.
[19]: #plot head head transect
plt.figure(figsize = (10,4))
x = np.arange(0,Lx,dx)
plt.plot(x,np.flipud(head[0])[int(nrow/2)][:])
plt.title('Head Transect across X-Domain',fontweight = 'bold', fontsize = 14)
plt.xlabel('X Distance (m)',fontsize = 12)
plt.ylabel('Head (m)',fontsize = 12)
plt.show()

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 102

Activities: - Increase pumping rates and observe the changes in head surface/drawdown extent. -
Add another pumping well variable, well_2, to the well package and rerun the model with two wells. -
Change the location of well_1 and observe different interactions with the constant head boundary
conditions. - Change the IBOUND to accomodate active flux boundaries on the top and bottom of
the model edges and observe differences in head surface symmetry.

Tutoriel Fopy 3 : Modèle transitoire

avec rivière
Dans ce script on créera un aquifère captif monocouche avec Flopy et on utilisera MODFLOW pour
exécuter une simulation de l'écoulement en régime transitoire des eaux souterraines dans un domaine
de 1 000 m x 1 000 m avec des conditions aux limites à charge hydraulique imposée sur les limites
gauche et droite de l'aquifère et une rivière dans le centre dont le niveau varie de façon transitoire.
Ce tutorial se concentre sur les packages MODFLOW RIV & CHD et les différences entre les états
transitoires et le modèle à l'état permanent créé dans ``Flopy Tutorial 1''. Des aspects similaires au
Tutoriel 1 seront créés très rapidement. Si vous ne comprenez pas une partie du code, n'hésitez pas
à revenir à ce notebook pour revoir ce qui se passe.

Import des packages et créer un objet modèle

[62]: #import packages
import flopy
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mp
import pandas as pd

#jupyter specific--included to show plots in notebook

%matplotlib inline

Remarque : vous devrez modifier la variable exen ame pour diriger Flopy vers le fichier exécutable
modflow sur votre ordinateur.
[63]: #create model object
modelname = "my_river_model"
#m = flopy.modflow.Modflow(modelname, exe_name = 'Exe/mf2005')

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 103

m = flopy.modflow.Modflow(modelname, exe_name = 'Exe/MODFLOW-NWT_64')

Discrétiser le modèle et joindre le package DIS

[64]: # Attribuer les variables de discrétisation
Lx = 1000.
Ly = 1000.
ztop = 0.
zbot = -50.
nlay = 1
nrow = 25
ncol = 25
dx = Lx/ncol
dy = Ly/nrow
dz = (ztop - zbot) / nlay

C'est ici où le modèle transitoire commence à être différent de la simulation en régime permanent :
• nper = 3 veut dire que la simulation va se faire pour 3 périodes de pompage.
• la variable steady se transforme en une liste avec des entrées booléennes ``True/False'' pour
chacune de ces périodes de pompage. La première période est désigné comme « True » afin
de mettre en place une solution en régime permanent pour notre aquifère à utiliser comme
base pour les futures périodes de de pompage en régime transitoire. (Si nous commençons
simplement à résoudre les périodes de pompage transitoires avec des hauteurs de départ
arbitraires, cela ne constituerait pas un système d'eaux souterraines naturel très réaliste.) Da-
vantage de périodes de périodes sont ajoutées à une simulation s'il y a des changements
dans le « pompage » dans l'aquifère, c'est-à-dire pompage, les niveaux des rivières ou les
taux de recharge. (Remarque : certaines choses ne peuvent pas être modifiées par période
de contrainte, comme par ibound, la conductivité hydraulique ou d'autres propriétés liées à la
géologie de l'aquifère.)
• perlen désigne la durée de chaque période de pompage. Pour la période de pompage de
régime permanent, ce que nous attribuons ici n'a pas vraiment d'importance puisque l'équation
ne dépend pas du temps. Pour les périodes de pompage transitoires, nous attribuons des
entrées « perlen » qui correspondent à ce que nous avons défini comme étant nos unités de
longueur. Pour ce modèle, chaque période de pompage durera 5 jours.
• nstp est le nombre de pas de temps pour lesquelles le solveur de MODFLOW atteindra des
solutions au cours de chaque période de pompage. Il peut être ajusté en fonction de ce
que vous souhaitez résoudre dans la simulation ; si vous souhaitez voir des changements
de haute résolution dans les profils de charge hydraulique, un nstp plus élevé peut être utilisé.
Cependant, cela augmentera le temps nécessaire pour résoudre la simulation puisque le solveur
doit générer plus de solutions, donc si vous êtes simplement intéressé par le comportement
général de la charge hydraulique sur un grand nombre de périodes de contrainte, un ``nstp''
plut petit peut être utilisé.
[43]: # Specifier le nombre de périodes de pompage
nper = 3

# Specifier si la période de pompage est transitoire ou permanent

steady = [True, False, False]
perlen = [1.0,5.0,5.0]
nstp = [1,10,10]

Nous créons le package dis incluant les paramètres de discrétisation temporelle et spatiale ci-dessus :

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 104

[44]: # Creer un objet Flopy de discretization, les longueurs et le temps sont␣

,→respectivement en mètre (2) and est en jour (4)

dis = flopy.modflow.ModflowDis(model=m, nlay=nlay, nrow=nrow, ncol=ncol,

delr=dx, delc=dy, top=ztop, botm=zbot,
itmuni = 4, lenuni = 2,
nper=nper, steady=steady, perlen=perlen, nstp=nstp)

[45]: # Vérifier la grille à partir du package _DIS_

modelmap = flopy.plot.PlotMapView(model=m, layer=0)
grid = modelmap.plot_grid()
plt.xlabel('$L_x$ (m)',fontsize = 14)
plt.ylabel('$L_y$ (m)',fontsize = 14)
plt.title('Grille de différences finies', fontsize = 15, fontweight = 'bold')
plt.show()

Specifier le type de maille et attacher le package BAS

[46]: # Créer la variable ibound comme un tableau ints = 1
ibound = np.ones((nlay, nrow, ncol), dtype=np.int32)

# Créer strt comme un tableau floats = 1m

strt = np.ones((nlay, nrow, ncol), dtype=np.float32)

# Créer l'objet _Flopy_ _bas_

bas = flopy.modflow.ModflowBas(m, ibound=ibound, strt=strt)

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 105

Définir les propriétés de la couche et attacher le package LPF

• Pour plus de simplicité, nous attribuons des propriétés de couche uniformes en tant que valeur
unique à appliquer à chaque maille.
[47]: # définir les propriétés de la couche
laytyp = 1
hk = 1.0 #m/j
vka = 1.0 #m/j
sy = 0.1 #1/m
ss = 1.e-4 #1/m

# attacher le package _LPF_

lpf = flopy.modflow.ModflowLpf(model=m, hk=hk, vka=vka, sy=sy, ss=ss,␣
,→laytyp=laytyp, ipakcb=53)

## V. Assign Constant Head Boundaries and Attach CHD Package

• The CHD package simulates constant head cells, and is used to assign boundary conditions
for this model rather than the ibound variable because the heads set in this package may be
adjusted transiently. Each Constant Head cell takes the following assignments:
cdh_cell = [layer, row, column, shead, ehead]
where: layer, row, column indicate the cell index and shead is the head in the cell at the start of the
stress period and ehead is the head in the cell at the end of the stress period. If these are different,
a linear interpolation is used to simulate head throughout time steps in the stress period.
• Constant head cells are initiated during a stress period by putting them in the chd_spd variable
which is set up as a dictionary of lists containing constant head info for each active cell during
the stress period:
chd_spd = {stress period: [chd_cell(s)]}
(More info in flopy documentation: http://modflowpy.github.io/flopydoc/mfchd.html)
• Here, we define a boundary condition that will carry through all the stress periods by default,
but it may be adjusted per stress period depending on what is most realistic for your model.

V. Attribuer des limites à charges hydrauliques constantes et attacher le pack

• Le package -CHD- simule des mailles à charge constante et est utilisé pour attribuer des
conditions aux limites pour ce modèle plutôt que la variable ibound car les charges définies
dans ce package peuvent être ajustées de manière transitoire. Chaque maille à charge
constante prend les affectations suivantes :
cdhc ell = [layer, row, column, shead, ehead]
Où: layer, row, column indiquent l'indece de la maill et shead est la charge hydraulique dans la
maille au début de la période de pompage et ehead est la charge dans la maille à la fin de la
période de pompage. Si ceux-ci sont différents, une interpolation linéaire est utilisée pour simuler la
charge à travers les pas de temps de la période de contrainte.
• Les mailles à charge constante sont initiées pendant une période de pompage en les mettant
dans la variable chds pd qui est configurée comme un dictionnaire de listes contenant des
informations sur les charges constantes pour chaque maille active pendant la période de
pompage :
chds pd = priode de stress : [chdc ell(s)]
• Ici, nous définissons une condition aux limites qui portera sur toutes les périodes de pompage
par défaut, mais elle peut être ajustée par période de pompage en fonction de ce qui est le
plus réaliste pour votre modèle.

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 106

[48]: # Attribuer les charges au début et à la fin de la période de pompage

strt_head=2
end_head=2
# Créer une liste pour contenir dans les mailles à charges constantes pendant la␣
,→période de pompage

bound_sp1 = []

# Attribuer une condition aux limites à charge constante aux mailles de gauche et␣
,→de droite du domaine

for lay in range(nlay):

for row in range(nrow):
bound_sp1.append([lay,row,0,strt_head,end_head])
bound_sp1.append([lay,row,ncol-1,strt_head,end_head])

print('Période de pompage 1 mailles CHD : \n',bound_sp1)

Période de pompage 1 mailles CHD :

[[0, 0, 0, 2, 2], [0, 0, 24, 2, 2], [0, 1, 0, 2, 2], [0, 1, 24, 2, 2], [0, 2,
0, 2, 2], [0, 2, 24, 2, 2], [0, 3, 0, 2, 2], [0, 3, 24, 2, 2], [0, 4, 0, 2, 2],
[0, 4, 24, 2, 2], [0, 5, 0, 2, 2], [0, 5, 24, 2, 2], [0, 6, 0, 2, 2], [0, 6, 24,
2, 2], [0, 7, 0, 2, 2], [0, 7, 24, 2, 2], [0, 8, 0, 2, 2], [0, 8, 24, 2, 2], [0,
9, 0, 2, 2], [0, 9, 24, 2, 2], [0, 10, 0, 2, 2], [0, 10, 24, 2, 2], [0, 11, 0,
2, 2], [0, 11, 24, 2, 2], [0, 12, 0, 2, 2], [0, 12, 24, 2, 2], [0, 13, 0, 2, 2],
[0, 13, 24, 2, 2], [0, 14, 0, 2, 2], [0, 14, 24, 2, 2], [0, 15, 0, 2, 2], [0,
15, 24, 2, 2], [0, 16, 0, 2, 2], [0, 16, 24, 2, 2], [0, 17, 0, 2, 2], [0, 17,
24, 2, 2], [0, 18, 0, 2, 2], [0, 18, 24, 2, 2], [0, 19, 0, 2, 2], [0, 19, 24, 2,
2], [0, 20, 0, 2, 2], [0, 20, 24, 2, 2], [0, 21, 0, 2, 2], [0, 21, 24, 2, 2],
[0, 22, 0, 2, 2], [0, 22, 24, 2, 2], [0, 23, 0, 2, 2], [0, 23, 24, 2, 2], [0,
24, 0, 2, 2], [0, 24, 24, 2, 2]]

[49]: # Créer un dictionnaire avec les donnée de la période de pompage

chd_spd={0: bound_sp1}

print('Données CHD Sde la période de pompage : \n', chd_spd)

Données CHD Sde la période de pompage :

{0: [[0, 0, 0, 2, 2], [0, 0, 24, 2, 2], [0, 1, 0, 2, 2], [0, 1, 24, 2, 2], [0,
2, 0, 2, 2], [0, 2, 24, 2, 2], [0, 3, 0, 2, 2], [0, 3, 24, 2, 2], [0, 4, 0, 2,
2], [0, 4, 24, 2, 2], [0, 5, 0, 2, 2], [0, 5, 24, 2, 2], [0, 6, 0, 2, 2], [0, 6,
24, 2, 2], [0, 7, 0, 2, 2], [0, 7, 24, 2, 2], [0, 8, 0, 2, 2], [0, 8, 24, 2, 2],
[0, 9, 0, 2, 2], [0, 9, 24, 2, 2], [0, 10, 0, 2, 2], [0, 10, 24, 2, 2], [0, 11,
0, 2, 2], [0, 11, 24, 2, 2], [0, 12, 0, 2, 2], [0, 12, 24, 2, 2], [0, 13, 0, 2,
2], [0, 13, 24, 2, 2], [0, 14, 0, 2, 2], [0, 14, 24, 2, 2], [0, 15, 0, 2, 2],
[0, 15, 24, 2, 2], [0, 16, 0, 2, 2], [0, 16, 24, 2, 2], [0, 17, 0, 2, 2], [0,
17, 24, 2, 2], [0, 18, 0, 2, 2], [0, 18, 24, 2, 2], [0, 19, 0, 2, 2], [0, 19,
24, 2, 2], [0, 20, 0, 2, 2], [0, 20, 24, 2, 2], [0, 21, 0, 2, 2], [0, 21, 24, 2,
2], [0, 22, 0, 2, 2], [0, 22, 24, 2, 2], [0, 23, 0, 2, 2], [0, 23, 24, 2, 2],
[0, 24, 0, 2, 2], [0, 24, 24, 2, 2]]}

[50]: # Créer un objet CHD de Flopy, et l'attacher au modèle

chd = flopy.modflow.ModflowChd(model=m, stress_period_data=chd_spd)

Définir l’emplacement et les propriétés de la rivière et joindre le package RIV

The MODFLOW simulation of a river may be visualized as the following:
A stage in the river is higher than the aquifer head in the cell directly below it, will induce ﬂow into

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 107

the aquifer depending on the head difference between the river and the aquifer and conductance
of the sediment bottom. If the reverse is true of the head and stage, water will ﬂow from the aquifer
to the river.

MODFLOW calculates flow between river and aquifer as the following: Qf lux = C ∗ ∆H
Where:
∆H = Hriver − Haquif er
k(dx ∗ dy)
C=
t
And:
L2
C = Conductance ( )
T
L
k = Conductivity of River Bed Sediment ( )
T
dx = Cell width (L)
dy = Cell Length (L)
t = Thickness of River Bed Sediment (L)

Specified Parameters: Location, Conductance, River Bottom, and River Stage are all specified for
each river cell in a MODFLOW model. In flopy, conductance is specified as cond , the river bottom
is specified as rbot and river stage is specified as r_stage . Layer, Row, and Column of the river cell
is also specified.
• The RIV Package takes stress period data very similarly to the CHD package with river cell
information defined in a list:
riv_cell = [layer, row, column, stage, conductance, bottom]
• And stress period data assigned in a dictionary of the form:
riv_spd = {stress period: [riv_cell(s)]}
(More info in flopy documentation: http://modflowpy.github.io/flopydoc/mfriv.html)
• Here, we define a column of river cells halfway across the domain and set the head in the river
to fluctuate between 1 and 5 during the stress periods.
La simulation MODFLOW d'une rivière peut être illustrée comme suit :
Sil le niveau de la rivière est plus élevé que la charge hydraulique de l'aquifère dans la maille située
juste en dessous, induira un écoulement vers l'aquifère en fonction de la différence de charges entre
la rivière et l'aquifère et la conductance du fond sédimentaire. Si c'est l'inverse, l'eau s'écoulera de
l'aquifère vers la rivière.

MODFLOW calcule le débit entre la rivière et l’aquifère comme suit :

Qf lux = C ∗ ∆H
Où:
∆H = Hriver − Haquif er
k(dx ∗ dy)
C=
t
Et:
L2
C = Conductance ( )
T
L
k = Conductivité des sédiments du lit de la rivière ( )
T

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 108

dx = Largeur de cellule (L)

dy = Longueur de cellule (L)
t = Épaisseur des sédiments du lit de la rivière (L)

Paramètres spécifiés :
La localisation, la conductance, le fond de la rivière et le niveau de la rivière sont tous spécifiés
pour chaque maille de la rivière dans le modèle MODFLOW. Dans Flopy, la conductance est spécifiée
comme cond, le fond de la rivière est spécifié comme rbot et le niveau de la rivière est spécifié comme
rs tage. La couche, la ligne et la colonne de la maille de la rivière sont également spécifiées.
• Le package RIV prend les données de période de pompage de manière très similaire au
package CHD avec des informations sur les maille de la rivière définies dans une liste :
rivc ell = [layer, row, column, stage, conductance, bottom]
• Et des données de période de pompage assignées dans un dictionnaire de la forme :
rivs pd = {période de pompage : [rivc ell(s)]}
• Ici, nous définissons une colonne de la maille de rivière à mi-hauteur du domaine et définissons
la charge dans la rivière pour qu'elle fluctue entre 1 et 5 pendant les périodes de pompage.
[51]: # Définir les rivières

# Mailles rivières de la période de pompage 1

riv_sp1 = [] # Créer une liste pour stocker toutes les mailles rivières pour la␣
,→période de pompage 1

k_rivbott = 1 # Conductivité hydraulique du lit de fond de la rivière e m/j

sed_thick = 1 # Epaisseur du lit de la rivière en m
cond = k_rivbott*(dy)*(dx)/(sed_thick) # Conductance du lit de fond de la rivière␣
,→m^2/j

r_stage = 1 # Niveau d'eau dans la rivière (période de pompage 1)

r_bott = 0 # Fond de la rivière
# Afffecter les données de la rivière aux mailles au centre la colonne
for i in range(nrow):
riv_sp1.append([0, i, ncol/2, r_stage, cond, r_bott])

# Mailles rivières de la période de pompage 2

riv_sp2 = [] #Créer une liste pour stocker toutes les mailles rivières pour la␣
,→période de pompage 2

r_stage = 5 #Niveau d'eau dans la rivière (période de pompage 2)

for i in range(nrow):
riv_sp2.append([0, i, ncol/2, r_stage, cond, r_bott])

# Mailles rivières de la période de pompage 3

riv_sp3 = [] #Créer une liste pour stocker toutes les mailles rivières pour la␣
,→période de pompage 3

r_stage = 1 #Niveau d'eau dans la rivière (période de pompage 3)

for i in range(nrow):
riv_sp3.append([0, i, ncol/2, r_stage, cond, r_bott])

# Créer le dictionnaire des donnée pour la période de pompage

riv_spd = {0: riv_sp1, 1: riv_sp2, 2: riv_sp3}

[52]: # Attacher le package rivière

riv = flopy.modflow.ModflowRiv(model=m,stress_period_data = riv_spd)

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 109

Check CHD and RIV Boundary Conditions

Vérifier les conditions aux limites CHD et RIV

[53]: # VERIFIER LES CONDITIONS AUX LIMITES
# Utiliser _Flopy_ pour représenter graphiquement la grille, le type de maille␣
,→ibound, les rivières et

# les conditions aux limites générales

modelmap = flopy.plot.PlotMapView(model=m, layer=0)
grid = modelmap.plot_grid()
ib = modelmap.plot_ibound()
riv_plot = modelmap.plot_bc(ftype='RIV')
chd_plot = modelmap.plot_bc(ftype='CHD')
# Ajouter les étiquettes et la légende
plt.xlabel('Lx (m)',fontsize = 14)
plt.ylabel('Ly (m)',fontsize = 14)
plt.title('Conditions aux limites', fontsize = 15, fontweight = 'bold')
plt.legend(handles=[mp.patches.Patch(color='green',label='Rivière',ec='black'),
mp.patches.Patch(color='navy',label='Limites à charges␣
,→constantes',ec='black'),

mp.patches.Patch(color='white',label='Maille active',ec='black'),
mp.patches.Patch(color='black',label='Maille␣
,→inactive',ec='black')],

bbox_to_anchor=(2,1))
plt.show(modelmap)

Spécifier la sortie des données et joindre le package OC

[54]: # Créer les données OC pour la période de pompage
oc_spd = {}
for kper in range(nper):
for kstp in range(nstp[kper]):
oc_spd[(kper, kstp)] = ['Enregistrer les charges',

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 110

'Enregistrer les rabattements',

'Enregistrer le bilant',
'Afficher les charges',
'Afficher le bilant']
oc = flopy.modflow.ModflowOc(model=m, stress_period_data=oc_spd,
compact=True)
print('Controle des sorties \n', oc_spd)

Controle des sorties

{(0, 0): ['Enregistrer les charges', 'Enregistrer les rabattements',
'Enregistrer le bilant', 'Afficher les charges', 'Afficher le bilant'], (1, 0):
['Enregistrer les charges', 'Enregistrer les rabattements', 'Enregistrer le
bilant', 'Afficher les charges', 'Afficher le bilant'], (1, 1): ['Enregistrer
les charges', 'Enregistrer les rabattements', 'Enregistrer le bilant', 'Afficher
les charges', 'Afficher le bilant'], (1, 2): ['Enregistrer les charges',
'Enregistrer les rabattements', 'Enregistrer le bilant', 'Afficher les charges',
'Afficher le bilant'], (1, 3): ['Enregistrer les charges', 'Enregistrer les
rabattements', 'Enregistrer le bilant', 'Afficher les charges', 'Afficher le
bilant'], (1, 4): ['Enregistrer les charges', 'Enregistrer les rabattements',
'Enregistrer le bilant', 'Afficher les charges', 'Afficher le bilant'], (1, 5):
['Enregistrer les charges', 'Enregistrer les rabattements', 'Enregistrer le
bilant', 'Afficher les charges', 'Afficher le bilant'], (1, 6): ['Enregistrer
les charges', 'Enregistrer les rabattements', 'Enregistrer le bilant', 'Afficher
les charges', 'Afficher le bilant'], (1, 7): ['Enregistrer les charges',
'Enregistrer les rabattements', 'Enregistrer le bilant', 'Afficher les charges',
'Afficher le bilant'], (1, 8): ['Enregistrer les charges', 'Enregistrer les
rabattements', 'Enregistrer le bilant', 'Afficher les charges', 'Afficher le
bilant'], (1, 9): ['Enregistrer les charges', 'Enregistrer les rabattements',
'Enregistrer le bilant', 'Afficher les charges', 'Afficher le bilant'], (2, 0):
['Enregistrer les charges', 'Enregistrer les rabattements', 'Enregistrer le
bilant', 'Afficher les charges', 'Afficher le bilant'], (2, 1): ['Enregistrer
les charges', 'Enregistrer les rabattements', 'Enregistrer le bilant', 'Afficher
les charges', 'Afficher le bilant'], (2, 2): ['Enregistrer les charges',
'Enregistrer les rabattements', 'Enregistrer le bilant', 'Afficher les charges',
'Afficher le bilant'], (2, 3): ['Enregistrer les charges', 'Enregistrer les
rabattements', 'Enregistrer le bilant', 'Afficher les charges', 'Afficher le
bilant'], (2, 4): ['Enregistrer les charges', 'Enregistrer les rabattements',
'Enregistrer le bilant', 'Afficher les charges', 'Afficher le bilant'], (2, 5):
['Enregistrer les charges', 'Enregistrer les rabattements', 'Enregistrer le
bilant', 'Afficher les charges', 'Afficher le bilant'], (2, 6): ['Enregistrer
les charges', 'Enregistrer les rabattements', 'Enregistrer le bilant', 'Afficher
les charges', 'Afficher le bilant'], (2, 7): ['Enregistrer les charges',
'Enregistrer les rabattements', 'Enregistrer le bilant', 'Afficher les charges',
'Afficher le bilant'], (2, 8): ['Enregistrer les charges', 'Enregistrer les
rabattements', 'Enregistrer le bilant', 'Afficher les charges', 'Afficher le
bilant'], (2, 9): ['Enregistrer les charges', 'Enregistrer les rabattements',
'Enregistrer le bilant', 'Afficher les charges', 'Afficher le bilant']}

Affecter PCG en tant que solveur de différences finies et attacher le package

[55]: # Affecter le solveur du problème d'écoulement souterrain
pcg = flopy.modflow.ModflowPcg(model=m)

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 111

Créer des fichiers -MODFLOW- et exécuter le modèle

[56]: #écrire les fichiers d'entrée de MODFLOW
m.write_input()

[57]: # Exécuter le modèle

success, mfoutput = m.run_model(pause=False, report=True)
if not success:
raise Exception('MODFLOW did not terminate normally.')

FloPy is using the following executable to run the model: .\mf2005.EXE

MODFLOW-2005
U.S. GEOLOGICAL SURVEY MODULAR FINITE-DIFFERENCE GROUND-WATER FLOW MODEL
Version 1.12.00 2/3/2017

Using NAME file: my_river_model.nam

Run start date and time (yyyy/mm/dd hh:mm:ss): 2019/06/06 15:05:00

Solving: Stress period: 1 Time step: 1 Ground-Water Flow Eqn.

Solving: Stress period: 2 Time step: 1 Ground-Water Flow Eqn.
Solving: Stress period: 2 Time step: 2 Ground-Water Flow Eqn.
Solving: Stress period: 2 Time step: 3 Ground-Water Flow Eqn.
Solving: Stress period: 2 Time step: 4 Ground-Water Flow Eqn.
Solving: Stress period: 2 Time step: 5 Ground-Water Flow Eqn.
Solving: Stress period: 2 Time step: 6 Ground-Water Flow Eqn.
Solving: Stress period: 2 Time step: 7 Ground-Water Flow Eqn.
Solving: Stress period: 2 Time step: 8 Ground-Water Flow Eqn.
Solving: Stress period: 2 Time step: 9 Ground-Water Flow Eqn.
Solving: Stress period: 2 Time step: 10 Ground-Water Flow Eqn.
Solving: Stress period: 3 Time step: 1 Ground-Water Flow Eqn.
Solving: Stress period: 3 Time step: 2 Ground-Water Flow Eqn.
Solving: Stress period: 3 Time step: 3 Ground-Water Flow Eqn.
Solving: Stress period: 3 Time step: 4 Ground-Water Flow Eqn.
Solving: Stress period: 3 Time step: 5 Ground-Water Flow Eqn.
Solving: Stress period: 3 Time step: 6 Ground-Water Flow Eqn.
Solving: Stress period: 3 Time step: 7 Ground-Water Flow Eqn.
Solving: Stress period: 3 Time step: 8 Ground-Water Flow Eqn.
Solving: Stress period: 3 Time step: 9 Ground-Water Flow Eqn.
Solving: Stress period: 3 Time step: 10 Ground-Water Flow Eqn.
Run end date and time (yyyy/mm/dd hh:mm:ss): 2019/06/06 15:05:00
Elapsed run time: 0.087 Seconds

Normal termination of simulation

Extraire les données de la charge hydraulique et de débit

[58]: # extraire les données binaires du fichier des charges
times = [perlen[0],perlen[0]+perlen[1],perlen[0]+perlen[1]+perlen[2]] # temps␣
,→d'extraction à la fin de chaque période de stress

head = {} # créer un dictionnaire pour stocker les données de la charge hydraulique␣

,→à la fin de chaque période de pompage

frf = {} # créer un dictionnaire pour stocker les flux à travers la face droite de␣
,→la cellule à la fin de chaque période de stress

fff = {} #créer un dictionnaire pour stocker les flux à travers la face avant de la␣
,→cellule à la fin de chaque période de stress

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 112

headobj = flopy.utils.binaryfile.HeadFile(modelname+'.hds') # obtenir les données␣

,→de la charge en tant qu'objet python

budgobj = flopy.utils.binaryfile.CellBudgetFile(modelname+'.cbc') #obtenir des␣

,→données de flux en tant qu'objet python

# obtenir des données à partir d'objets python

for stress_per, time in enumerate(times): # itérer à travers le temps à la fin de␣
,→chaque période de pompage

head['sp%s'%(stress_per)] = headobj.get_data(totim=time) # ajouter tête à tête␣

,→liste pour chaque période de pompage

frf['sp%s'%(stress_per)] = budgobj.get_data(text='FLOW RIGHT FACE',totim=time)␣

,→# ajouter le flux de la face droite à la liste frf pour chaque période de stress

fff['sp%s'%(stress_per)] = budgobj.get_data(text='FLOW FRONT FACE',totim=time)␣

,→# ajouter le flux de la face avant à la liste fff pour chaque période de pompage

Représentation graphique
Graphiques des lignes piézométriques et de flux

[59]: # tracer les résultats pour toutes les périodes de période

for i in range(len(times)):
plt.figure(figsize=(9,9)) #Créer une figure 10 x 10
modelmap = flopy.plot.PlotMapView(model=m, layer=0) # utiliser modelmap pour␣
,→attacher le tracé au modèle

grid = modelmap.plot_grid() # Représenter la grille du modèle

riv_plot = modelmap.plot_bc(ftype='RIV') # Représenter les mailles rivières
chd_plot = modelmap.plot_bc(ftype='CHD') # Représenter les mailles ghb
contour_levels = np.linspace(head['sp%s'%i][0].min(),head['sp%s'%i][0].
,→max(),15) # définir les niveaux de contour pour la charge de contour

head_contours = modelmap.contour_array(head['sp%s'%i][0],␣
,→levels=contour_levels) # Créer les lignes piézométriques

flows = modelmap.plot_discharge(frf['sp%s'%i][0], fff['sp%s'%i][0],␣

,→head=head['sp%s'%i]) # Créer les flèches du flux

# Afficher les paramètres

plt.xlabel('Lx (m)',fontsize = 14)
plt.ylabel('Ly (m)',fontsize = 14)
plt.title('Résultat Charge/Débit : période de pompage %s'%(i+1), fontsize = 13,␣
,→fontweight = 'bold')

plt.colorbar(head_contours, aspect=5)
plt.legend(handles=[mp.patches.Patch(color='green',label='Rivière',ec='black'),
mp.patches.Patch(color='navy',label='Limite à charge␣
,→imposée',ec='black'),

mp.patches.Patch(color='white',label='Mailles␣
,→actives',ec='black')],

bbox_to_anchor=(1.8,1.0))
plt.show(modelmap)

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 113

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 114

Tracés 3D de la surface des charges hydrauliques

Remarque : Flopy prend l'indexation des mailles où le coin supérieur gauche de la grille est l'index de
la maille (0,0) tandis que les unités de longueur et de largeur commencent à 0 dans le coin inférieur
gauche de la grille. Cela semble un peu étrange, mais cela permet de faire facilement correspondre
un tableau de numéros de cellule par cellule à leurs emplacements respectifs sur la grille lors de
l'attribution de propriétés ou de l'observation des données de sortie. Les capacités de traçage
de Flopy retournent automatiquement les données de tête résultantes pour les afficher sur sa grille
appropriée, cependant, pour les graphiques 3D ci-dessous, vous verrez que la fonction np.f lipud()
est utilisée pour retourner le tableau de données pour tracer dans la même direction.
[60]: #importer la boîte à outils des axes 3D de matplotlib
from mpl_toolkits.mplot3d import Axes3D

# définir les variables X, Y, Z pour le tracé 3D comme notre domaine de modèle et␣
,→notre solution principale

X = np.arange(0,Lx,dx)
Y = np.arange(0,Ly,dy)
X, Y = np.meshgrid(X, Y)

for i in range(len(times)):
# créer une figure 3D
fig_3d = plt.figure(figsize=(12,5))
ax = fig_3d.gca(projection='3d')
Z = np.flipud(head['sp%s'%i][0])

# créer une légende et des étiquettes

surf = ax.plot_surface(X,Y,Z, cmap = plt.cm.coolwarm, linewidth=0,␣
,→antialiased=False, label='head')

fig_3d.colorbar(surf,shrink=0.5,aspect=5).set_label('Charge hydraulique␣
,→(m)',fontsize=10,fontweight='bold')

ax.set_xlabel('$L_x$ (m)', fontsize=15, fontweight='bold')

ax.set_ylabel('$L_y$ (m)', fontsize=15, fontweight='bold')

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 115

ax.set_title('Surface des charges: Période de pompage %s'%(i+1), fontsize=15,␣

,→fontweight='bold')
plt.show(surf)

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 116

• Notice the difference between the plots for Stress-Periods 1 and 3: They both have the same
boundary conditions, however, plot 1 displays the steady-state solution, whereas plot 3 repre-
sents conditions that have not entirely relaxed to the steady-state drawdown. Increasing the
length (perlen) of this last stress period should yield a head surface that looks more similar to
stress period 1.
• Notez la différence entre les représentations graphiques pour les périodes de pompage 1
et 3 : ils ont tous les deux les mêmes conditions aux limites, cependant, le graphe 1 affiche
la solution à l'état permanent, tandis que le graphe 3 représente les conditions qui ne se
sont pas entièrement stavilisées jusqu'au rabattement à l'état permanent. L'augmentation de
la longueur (perlen) de cette dernière période de contrainte devrait donner une surface de
charges hydrauliques qui ressemble plus à la période de pompage 1.

Tracé des séries temporelle de la charge dans une seule maille

get time series for cell
[61]: # tracer une série chronologique dans la maille à gauche de la rivière
# obtenir des séries temporelles pour la cellule
cell_id = (0, int(nrow/2) - 1, int(ncol/2) - 1)
time_series = headobj.get_ts(cell_id)

# Créer le tracé
plt.subplot(1, 1, 1)
plt.title('Charge hydraulique à la maille ({0},{1},{2})'.format(cell_id[0] + 1,
cell_id[1] + 1,
cell_id[2] + 1),fontweight='bold')
plt.xlabel('Temps (days)',fontweight='bold')
plt.ylabel('Charge (m)',fontweight='bold')
plt.plot(time_series[:, 0], time_series[:, 1], 'bo-')
plt.show()
print('Données sur la série temporelle de la charge hydraulique: \n', time_series)

Ali HAMMANI, 2021

Utilisation du package Python Flopy pour la modélisation des écoulement souterrains 117

Time Series Head Data:

[[ 1. 1.0880829]
[ 1.5 3.2954612]
[ 2. 3.9382052]
[ 2.5 4.2017574]
[ 3. 4.343931 ]
[ 3.5 4.4343276]
[ 4. 4.497604 ]
[ 4.5 4.5445633]
[ 5. 4.5807242]
[ 5.5 4.6092424]
[ 6. 4.632088 ]
[ 6.5 2.4432013]
[ 7. 1.8155264]
[ 7.5 1.56431 ]
[ 8. 1.4322644]
[ 8.5 1.3502005]
[ 9. 1.2937887]
[ 9.5 1.252489 ]
[10. 1.2209972]
[10.5 1.196332 ]
[11. 1.1766671]]

Ali HAMMANI, 2021

Utilisationconjuguéedeseaux
Chapitre 6

desurfaceetdeseauxsouter-
raines

Introduction
Dans les régions arides et semi-arides l'irrigation permet l'exploitation des terres agricoles pendant
les périodes de l'année où l'évapotranspiration excède les précipitations. Le développement de
l'irrigation s'accompagne souvent de problèmes d'engorgement des sols par remontée de la nappe
et de problèmes de salinité rendant les sols non productifs. Par ailleurs, la surexploitation des
eaux souterraines dans certaines régions a entraîné un rabattement excessif des nappes souter-
raines (Souss-Massa) et l'intrusion marine (région de Bir Jdid). Cela est généralement le résultat du
développement séparé de l'exploitation des ressources en eau de surface et des ressources en
eau souterraine. L'utilisation conjuguée des eaux de surface et des eaux souterraines combine les
avantages peut constituer une manière efficace pour la gestion efficace des ressources en eau.
Cependant, les eaux souterraines et les eaux de surface doivent être gérées d'une manière opti-
male dans l'objectif de minimiser les effets indésirables de l'utilisation séparée de chacune des deux
ressources.
L'utilisation conjuguée est définie comme étant une gestion coordonnée, dans le temps et dans
l'espace, de l'eau de surface et de l'eau souterraine pour tirer profit de leurs propriétés complémen-
taires. Les Le terme `` conjuguée ou conjointe '' est utilisé lorsque les deux sources qui sont utilisées
séparément, sont utilisées pour un meilleur effet sous une gestion intégrés des lâchers d'eau de sur-
face à partir des barrages, des dérivations de rivière et d'exploitation des eaux souterraines. Sous
des programmes d'utilisation conjuguée, les eaux de surface et les eaux souterraines sont utilisées
pendant les années à déficit pluviométrique pour des apports supplémentaires d'eau d'irrigation.
L'utilisation conjuguée peut être également envisagée pour augmenter le débit des canaux et avoir
une main d'eau suffisante. Lorsque les eaux souterraines sont de mauvaise qualité l'utilisation con-
juguée peut être une manière de l'améliorer par son mélange avec les eaux de surface. Pendant
les périodes où la pluie dépasse nettement l'évapotranspiration les eaux de surface sont utilisée au
maximum et l'excédant est utilisé pour la recharge des nappes souterraines.
L'utilisation conjuguée des eaux de surface et des eaux souterraines ne consiste pas uniquement
à l'utilisation des eaux souterraines en plus des eaux de surface mais consiste en une exploitation
optimale des deux ressources en considérant la dynamique des terres, des eaux et des systèmes
écologiques.
Un autre aspect de l'utilisation conjuguée consiste en l'utilisation mélangée de l'eau de surface de
bonne qualité et de l'eau souterraine de mauvaise qualité pour réduire la salinité du mélange à une
limite acceptable.

118
Utilisation conjuguée des eaux de surface et des eaux souterraines 119

La gestion des ressources en eau de sur-

face et souterraines
Gestion des ressources en eau
La demande en eau est de plus en plus croissante dans les différents secteurs (agriculture, eau
potable, industrie, tourisme…(figure ??) alors que l'offre est en réduction constante notamment dans
les pays aride est semi aride. Si l'utilisation des eaux de surface peut directement être maîtrisée par
les pouvoirs publics, les ressources en eaux souterraines restent mal contrôlées du moment ou la mise
en place des ouvrages de captage est individuelle.

L'eau souterraine est une ressource extrêmement importante pour beaucoup de pays en développe-
ment. En Inde, par exemple, les puits tubulaires donnaient à la fin des années 80 plus de la moitié de
l'eau utilisée sur la superficie irriguée nette. Du fait que la gestion des couches aquifères doit tenir
compte de l'interaction complexe entre la société et l'environnement physique, elle pose des prob-
lèmes redoutables au stade de la conception des politiques. En présence de couches aquifères
surexploitées, la gestion ou la réglementation doit répondre à deux catégories de choix collectifs.
D'une part, pour ce qui est de gérer l'eau, les décisions se fondent sur :
• le taux annuel approprié de pompage;
• la distribution géographique du pompage;
• la possibilité de complément par les eaux superficielles;
• la possibilité de la recharge artificielle de la couche aquifère.
D'autre part, l'autre catégorie de décisions, à savoir coordonner le pompage, détermine:

Ali HAMMANI, 2021

Utilisation conjuguée des eaux de surface et des eaux souterraines 120

• quelles doivent être les institutions et les politiques qui répartissent le taux d'extraction entre
utilisateurs individuels potentiels et classes d'utilisateurs, et influencent les modalités de pom-
page;
• comment les règles qui limitent le pompage sont appliquées et leur respect surveillé.
L'exploitation des eaux souterraines est traditionnellement extensive et s'attache encore souvent à
découvrir et à utiliser des ressources nouvelles. Cette exploitation doit devenir moins agressive pour
les milieux naturels, plus économe, et à l'échelle du système aquifère, permettre la multiplication des
utilisations en veillant à préserver la ressource, que ce soit dans son volume ou dans sa qualité. Le
recours aux ressources nouvelles doit être concerté et doit tenir compte du droit des générations
futures : la gestion doit intégrer le long terme.

Les problèmes liés à la gestion des ressources en eau

Les problèmes de gestion de l'eau commencent dès que le volume pompé dépasse les possibilités
de recharge naturelle. On observe alors en premier lieu un épuisement de plus en plus rapide des
réserves d'eau souterraine et un accroissement en conséquence des coûts de pompage; l'intrusion
d'eau de moindre qualité dans les couches que l'on exploite; des intrusions d'eau saline sous l'effet
d'un pompage à trop fort débit à proximité des côtes; enfin, des dépôts minéralisés alternant avec
des eaux de meilleure qualité.
Dans une perspective générale, l'exploitation des couches aquifères peut provoquer, séparément
ou simultanément, des dilemmes sociaux de deux types. D'une part, la surexploitation entre dans la
catégorie des problèmes liés aux ressources, que l'on baptise souvent problèmes de ``bien commun''.
Une ressource de bien commun se définit par deux caractéristiques :
• le principe de soustraction (ce qui signifie qu'une unité de ressource utilisée par un usager n'est
plus disponible pour un autre);
• le coût élevé d'exclusion d'usagers potentiels de l'exploitation de cette ressource.
Les problèmes liés aux ressources de bien commun ont leur origine dans l'insuffisance du cadre
économique et institutionnel dans lequel sont exploitées lesdites ressources. Les ressources de bien
commun sont classiquement utilisées dans une perspective de ``libre accès'', dans lequel elles sont
exploitées selon des règles de capture. Quand personne n'est propriétaire de la ressource, les
usagers ne sont nullement incités à la conserver en prévision de l'avenir, et l'intérêt égoïste de
l'individu porte celui-ci à la surexploitation. La question fondamentale posée par la gestion des
ressources de bien commun est celle de la définition des institutions économiques qui gouvernent
leur usage.

Les avantages des eaux souterraines

Les eaux souterraines peuvent produire une ressource en eau supplémentaire ; aussi bien comme
moyen de stockage que de distribution et de traitement, qui peut être avantageusement combinée
aux ressources en eau de surface et à leurs aménagements. Six différents ressources gérables peuvent
être considérées dans un bassin hydrogéologique : la production sûre, le volume des réserves (volume
potentiellement exploitable), la capacité de stockage à long terme, le système de transfert, la qualité
de l'eau et la ressource d'énergie.
Pour relever les atouts de l'utilisation conjuguée de l'eau de surface et de l'eau souterraine on peut
évoquer l'originalité des eaux souterraines qui sont caractérisées par :
• La possibilité d'une extension dans l'espace facilitant le captage des eaux au droit des points
d'utilisation, ce qui n'est pas le cas de l'eau de surface.
• Les aquifères sont des réservoirs naturels de long terme et sont à l'abri de l'évapotranspiration
(excepté pour les nappes phréatiques peu profondes).
• La réponse des eaux souterraines, A la différence des eaux superficielles, est moins sujette aux
aléas climatiques, puisque l'eau souterraine, stockée au niveau de la nappe phréatique ou

Ali HAMMANI, 2021

Utilisation conjuguée des eaux de surface et des eaux souterraines 121

une nappe plus profonde, n'est pas directement soumise l'influence de facteurs climatiques, en
l'occurrence la température qui stimule l'évaporation et la demande climatique.
• Les ouvrages de captage de l'eau souterraine occupent des espaces très réduits, et l'effort
de leur réalisation est moins important que celui du captage des eaux superficielles. Par
conséquent, les investissements, pour l'exploitation et l'accès A 1'eau souterraine, sont réduits
et peuvent être à la portée des petites collectivités et même des individus qu'ils ne le sont
pour les eaux de surface. Toutefois. ceci n'est pas le cas, lorsqu'un programme de recharge
artificielle de la nappe est envisagé.
Par ailleurs, le même auteur signale qu'une sur-exploitation de la nappe, qui, dans certaines situa-
tions, peut conduire par exemple à une intrusion d'eau salée, compromettrait ainsi la durabilité de
l'utilisation conjuguée des eaux de surface et des eaux souterraines.
D'autre part, si l'investissement de base pour l'exploitation de l'eau souterraine, est réduit par rap-
port à celui de l'eau de surface, le coût d'exploitation du m3 d'eau souterraine est généralement
supérieur du fait qu'il est supporté individuellement par les agriculteurs. Par contre, l'investissement
pour l'ouvrage de captage de l'eau de surface fait profiter un grand nombre d'agriculteurs qui
se trouvent à l'aval de l'ouvrage. Par conséquent, cette disparité des prix peut être une con-
trainte financière la réussite de programmes d'utilisation conjuguée, et à une gestion optimale des
ressources en eau. Cette contrainte peut être réduite par un effort conscient d'équilibrer les prix des
eaux souterraines et des eaux de surface, comme c'est le cas dans une partie du périmètre irrigué
du Tadla où le coût de mètre cube d'eau d'irrigation est le même lorsque l'agriculteur pompe les
eaux souterraines ou lorsqu'il recours à l'eau de surface.

Avantages du stockage souterrain de l’eau

• Les réservoirs (barrages) de surface sont mis en place d'une manière définitive. Avec le temps,
il y aura plus de site convenable pour la mise en place de nouveau réservoir. Le stockage
souterrain de l'eau par des recharges artificielle serait plus sûr ;
• La qualité physique et chimique de l'eau souterraine est plu uniforme que l'eau de surface. Les
eaux souterraines sont moins sujettes à la pollution que les eaux de surface. Même si elles sont
polluées, les polluants sont éliminés ou dilués lors du mouvement de l'eau. Les eaux souterraines
ne sont pas chargées en sédiments et ne subissent pas des fluctuations de température comme
les eaux de surface ;
• Le stockage de l'eau souterraine peut se faire sans aucune perte de terres agricoles ;
• L'eau souterraine peut être utilisée à n'importe quel endroit et à n'importe quel moment sans
risque de fuite et de perte par évaporation durant la période de stockage ;
• Il y a moins de risques écologiques le stockage des eaux souterraines que pour le stockage
des eaux de surface ;
• Le coût de stockage des eaux souterraines est plus faible que celui des eaux de surface.

Le problème générés par une utilisation non raisonnée des ressources en eau
Salinité des eaux et du sol
La salinisation se produit sous l'effet combiné d'un mauvais drainage et d'un taux élevé d'évaporation,
qui entraîne une concentration les sels dans les terres irriguées; le phénomène se produit principale-
ment dans les régions arides ou semi-arides. Même une eau d'irrigation de bonne qualité contient
une certaine quantité de sels en solution, et peut en abandonner plusieurs tonnes par hectare
chaque année. A moins que ces sels ne soient entraînés par l'eau au-dessous de la couche dans
laquelle plongent les racines, le sol devient salin. Plusieurs facteurs influencent la salinité, notamment
la profondeur de la nappe phréatique, les caractéristiques capillaires du sol et les pratiques de
conduite de l'irrigation, notamment la quantité d'eau apportée en excès de l'évapotranspiration
effective des plantes pour lessiver les sels.

Ali HAMMANI, 2021

Utilisation conjuguée des eaux de surface et des eaux souterraines 122

Dans le périmètre irrigué du Tadla, il se trouve actuellement affecté par la salinité des eaux souter-
raines, principalement à l'aval hydraulique des périmètres des Béni-Amir et des Beni-Moussa ainsi
qu'au droit des zones à problème de remontée de la nappe. La qualité chimique des eaux souter-
raines du périmètre présente une très grande irrégularité du taux de salure qui varie de moins de
0,5 g/l dans les Béni- Moussa Est à plus de 4g/l dans les Beni-Moussa de l'Ouest et les Beni-Amir.

La surexploitation des ressources en eau souterraines

La surexploitation des nappes est un problème très préoccupant tant par son ampleur que par les
difficultés qu'a l'administration à faire respecter la réglementation et par ses capacités à enrayer le
phénomène.
La qualité des eaux souterraines se dégrade au fur et à mesure de leur exploitation et surtout de
leur surexploitation. La désalinisation deviendra de plus en plus une nécessité.
Les conséquences de la surexploitation mettent un certain temps à être perçues car le stock
disponible d'eau dans les aquifères représente en général plusieurs fois les apports. Aussi, il est
difficile de sensibiliser les exploitants aux risques qu'ils encourent et qu'ils font encourir à l'ensemble
des autres exploitants d'un même aquifère. A terme, les surexploitations vont conduire à un abaisse-
ment des niveaux piézométriques et inexorablement à un tarissement des puits ou des forages ou
précédemment à une salinisation des eaux due à des intrusions marines pour les aquifères côtiers ou
au déplacement du front salin à proximité des Chotts. Certes, les années excédentaires arrivent à
rétablir les équilibres pour certaines nappes, mais il faut craindre les phénomènes de dégradation
irréversibles.
Les rythmes d'implantation des moyens d'exhaure et d'accroissement des volumes pompés conduisent
à prévoir un développement rapide de la surexploitation et de ses conséquences.
Au Maroc, l'exemple le plus marquant de la sur-exploitation des eaux souterraines et ce lui de la
nappe de Souss où une intensification agricole a été opérée pendant très longtemps. La nappe
de Saïss (Région de Fès-Meknes) commence à connaître le même sort. Dans la région de Bir Jdid
(entre Casablanca et Azemmour), un exploitation non planifiée des eaux souterraines a entraîné
un rabattement excessif de la nappe induisant une forte intrusion marine et l'abandon de l'activité
agricole (figure ??).

Dans le périmètre irriguée de Tadla, après une période de remontée excessive de la nappe phréa-
tique pendant les années suivant la mise en eau du périmètre, l'eau souterraine a constitué une
ressource supplémentaire utilisée notamment pendant les années sèche par la mise en place de
puits et forages (figure ??). Le cas de Tadla est un cas où l'on peut mettre en œuvre l'utilisation
conjuguée des eaux de surface et des eaux souterraines. La nappe peut être rechargé pendant
les périodes pluvieuses.

Ali HAMMANI, 2021

Utilisation conjuguée des eaux de surface et des eaux souterraines 123

Gestion d’une utilisation conjuguée des

eaux de surface et des eaux souter-
raines
Eléments de l’utilisation conjuguée
La planification d'une utilisation conjuguée doit prendre en considération les spécificités locales
de chaque zone et certains éléments doivent être listés pour chaque zone. Les éléments généraux
d'une utilisation conjuguée consiste en une étude des potentialités des ressources en eau (de sur-
face et souterraines) leur distribution, les dispositifs de drainage, l'étendu et la capacité des réseaux
de canaux, les changement dynamiques dans le régime et la qualité des ressources en eaux, les
conditions agro-climatiques, les besoins en eau dans l'espace et dans le temps, les impacts environ-
nementaux et écologiques ainsi que les possibles changements socio-économiques.
L'utilisation conjuguée nécessite la mise en place d'équipement :
• La distribution de l'eau ;
• La recharge artificielle ;
• Le pompage de l'eau souterraine.
La mise en place de ces équipements nécessite une connaissance détaillée des paramètres de
l'aquifère, des estimations des débits de pompage et recharge, des informations sur le niveau des
nappes et de la qualité des eaux souterraines, etc. L'utilisation optimale des ressources en eau
disponible nécessite en plus un modèle pour simuler la réponse du bassin à la variation de la
recharge naturelle et artificielle, au pompage et à d'autres stratégies alternatives.
D'autres paramètres doivent également connus pour la planification de l'utilisation conjuguée. Il
s'agit notamment de schéma d'occupation des sols et des assolements, la capacité des canaux,
la capacité des stations de pompage, l'espacement entre les puits, l'étendu des revêtement des
canaux, les besoins en drainage, le niveau minimal auquel la nappe doit être maintenu, etc. Ceci
nécessite généralement des simulations pour tester différents scénarios de conception de l'utilisation
conjuguée.
Trois composantes importantes peuvent être définie pour une planification de l'utilisation conjuguée :
• Les apports en eau de surface et son transport ;
• Les apports en eau souterraine et rabattement de la nappe ;

Ali HAMMANI, 2021

Utilisation conjuguée des eaux de surface et des eaux souterraines 124

• La nature des terres agricoles.

Les précipitations sont les sources principales des apports d'eau. Elles peuvent être stockées dans
les réservoirs, dans les petits ou grands barrages selon l'intensité de la pluie et l'extension du bassin
versant. Les eaux de surface peuvent également être disponibles par des transferts entre bassins
versants. Elles peuvent être transférées au point de leur utilisation par un réseau de canaux revêtus
ou non revêtus.
Le stockage de l'eau souterraine diffère quantitativement et qualitativement en fonction du temps
et de l'espace selon la disponibilité des aquifères, le taux de recharge et la qualité inhérente de
l'eau souterraine. L'exploitation des eaux souterraines à partir de nappes peu ou très profondes
contenant de l'eau douce ou salée diffère très largement. La conception d'ouvrage de captage
des eaux souterraines peut varier d'un puits traditionnel individuel à de forages profonds et peu
profonds individuel ou collectif. Les puits individuels pour l'irrigation peuvent avoir leur propre système
de distribution de l'eau pompée. Ce genre de puits a plusieurs avantages comme par exemple un
système séparé de fourniture en énergie. Les puits ou forages collectifs qui ont une plus grande
capacité sont plus économiques comparés aux puits de faible capacité. Cependant, ils présentent
un problème de gestion.
La nature des terres agricoles et différentes situations qui peuvent lui être associées influence aussi
la planification de l'utilisation conjuguée. Les facteurs liés aux terres agricoles comprennent égale-
ment le regroupement au non des exploitations agricoles et le degré de modernisation de ces
exploitations.

Etapes de planification de l’utilisation conjuguée

Les étapes d'une planification conjuguée des eaux de surface et des eaux souterraines peuvent
être résumées comme suit :
1. Evaluer les situations hydrologiques et quantification des différentes composantes du bilan
hydrologique dans le périmètre irrigué ;
2. Collecter, traiter et analyser toutes les informations ou données disponibles sur l'occupation
des sols, les assolements, les besoins en eau, les conditions hydrologique et météorologique.
3. Evaluer la disponibilité des eaux souterraines et des eaux de surface ainsi que la demande en
eau. Puis adapter la demande à l'offre ;
4. Identifier les zones critiques de point de vue engorgement des sols et salinité ;
5. Développer un modèle numérique pouvant simuler les conditions hydrogéologiques dans la
zone d'étude et qui peut générer des scénarios ;
6. Elaborer des programmes et des stratégies de développement des ressources en eau souter-
raine conjointement avec les eaux de surface pour arriver à un développement optimal des
ressources en eau ;
7. Faire des analyses des coûts/bénéfices du programme de développement de l'utilisation con-
juguée.

Les avantages de l’utilisation conjuguée des eaux de surface et des eaux sout
L'utilisation conjuguée des eaux de surface et des eaux souterraines pour l'irrigation combine les
avantages l'exploitation des eaux souterraines et ceux de l'exploitation des eaux de surface et sert
de mesure corrective pour une gestion et une utilisation efficace des ressources en eau. Les avan-
tages majeurs de l'utilisation de l'utilisation conjuguée des eaux de surface et des eaux souterraine
peuvent être résumés comme suit :
• L'utilisation des eaux souterraines contribue dans la réduction de la demande de point pour
l'irrigation, à la réduction de la taille des canaux et par conséquent des coûts ;
• Les apports supplémentaires d'eau souterraine garantie une programmation propre des irriga-
tions, permet la diversification des systèmes de culture et mis en culture précoce des terres ;

Ali HAMMANI, 2021

Utilisation conjuguée des eaux de surface et des eaux souterraines 125

• L'eau souterraine accroît et garanti la disponibilité des ressources en eau dans les queues de
réseau ;
• L'exploitation des eaux souterraines rabat les niveaux des nappes phréatiques et réduit les
risques d'engorgement des sols et par conséquent d'éviter le gaspillage d'eau dans le cas où
il faut éliminer les excès d'eau par drainage ;
• L'accès à l'eau souterraine permet de réduire le volume des bassins de stockage et donc
réduction des coûts des projets de goutte à goutte ;
• Avec l'utilisation des eaux conjuguée, les écoulements de surface sont minimisés entraînant une
réduction des débits de crues ;
• Quand l'utilisation conjuguée est intégrée à un projet de recharge de nappe souterraine, les
besoins en revêtement des canaux sont réduits du fait que les infiltrations à partir des canaux
recharge la nappe ;
• L'utilisation conjuguée permet de réduire les pertes en d'eau par évaporation directe à partir
des retenues des eaux de surface.

Contraintes de la mise en œuvre d’une utilisation conjuguée

• Possibilité de la dégradation de la qualité des eaux souterraines due au lessivage vers la
nappe des sels du sol. Cette dégradation peut aussi résulter de transferts latéraux ou verticaux
des sels à partir de nappes adjacentes ou soujacente vers les zones de pompage où la qualité
de l'eau est meilleure ;
• Augmentation de la consommation d'énergie pour le pompage dans des eaux dans les puits. Il
y a risque de perturbation des apports d'eau souterraine suite à un manque d'énergie pendant
les périodes critiques. La fluctuation du niveau de la nappe peut aussi causer une réduction
des rendements des pompes ;
• La mise en œuvre, la supervision et le contrôle de l'utilisation conjuguée et de la recharge
artificiel sont très complexes ;
• Les difficultés administratives dans la répartition acceptable et équitable de l'eau sachant
qu'il faut prévoir des motivations et des incitations pour l'utilisation de l'eau souterraine lorsque
l'eau de surface est disponible.

Contrôle de l’utilisation de l’eau

Il est certain, qu'on ne peut parler de la possibilité d'une coordination efficace de l'utilisation de
l'eau de surface et de l'eau souterraine dans le temps et dans l'espace pour l'irrigation sans évoquer
le problème de contrô1e de l'utilisation de cette eau. Ce contrô1e a, en effet, reçu une attention
particulière de la part des spécialistes du fait de ses implications socio-politiques et économiques.
Ainsi, deux approches principales relatives à l'utilisation conjuguée pour le contrô1e de l'exploitation
des ressources en eaux souterraines ont été présentées; celle d'un contrô1e public ou puits collectifs
et celle d'un contrôle individuel ou puits privés.

Puits collectifs
Cette approche concerne l'instauration de puits collectifs qui vont être exploités à 1'échelle ré-
gionale pour une utilisation conjuguée avec l'eau de surface. Ces puits, gérés par la collectivité,
vont permettre d'approvisionner les exploitations agricoles en eau, par l'intermédiaire d'un réseau
d'irrigation.
L'expérience Pakistanaise avec les puits publics, mérite une attention particulière. En effet, à la suite
de problèmes d'engorgement et de salinisation qu'à connu le périmètre Indus (Pakistan), plusieurs
chercheurs et experts de la Banque Mondiale ont recommandé l'installation de puits publics, dans
le but de réaliser un drainage vertical des sols par le rabattement du niveau de la nappe en
permettant par la suite d'amé1iorer la productivité des sols.

Ali HAMMANI, 2021

Utilisation conjuguée des eaux de surface et des eaux souterraines 126

Ainsi, même si le projet de démonstration du contrô1e de la salinité et de l'engorgement des sols a

montré aux agriculteurs que le niveau de la nappe peut être rabattu par une distribution uniforme
des puits, ceci n'a pas empêché des milliers d'entre eux d'investir dans des puits privés, en se limitant
uniquement aux zones où la qualité de la nappe est bonne. De ce fait, 1'échec de l'expérience
Pakistanaise a amené les experts à nuancer, de plus en plus, la recommandation d'une telle alter-
native. On estime que cet échec est dû d'une part à la participation effective quasi absente des
agriculteurs (les premiers concernés) et d'autre part à la sous estimation des difficultés de la mise
en oeuvre et de la gestion d'un projet de puits collectifs. De plus, le comportement des agriculteurs
est compréhensible du fait qu'ils ne cherchent qu'à avoir un contrô1e personnel et donc sûr sur la
ressource en eau.
Au Maroc une expérience similaire des puits collectifs a été menée dans le Tadla au périmètre des
Beni Moussa. Un système de drainage vertical par des batteries de puits a été installé dans une
zone reconnue par le problème de remontée de nappe. L'eau pompée, est rejetée dans les canaux
d'irrigation en complément de l'eau de surface. Cette expérience, a donné de bons résultats d'une
part, par le rabattement du niveau de la nappe qu'elle assure et d'autre part par le complément
d'eau d'irrigation qu'elle fourni, qui arrive a satisfaire 45% des besoins d'irrigation d'un secteur ayant
une surface d'environ 1000 ha.

Puits individuels
Les agriculteurs ont une préférence pour un contrô1e individuel ou privé sur les approvisionnements
en eau. En effet, signale que depuis que l'accès A la ressource en eau est déterminant pour la
stabilité des revenus des agriculteurs dans les régions où l'eau est une ressource rare, ces derniers
préfèrent, naturellement, avoir un contrô1e sur cette ressource.
Par conséquent, le comportement pragmatique des agriculteurs pour la recherche d'un contrô1e
sur les ressources en eau explique l'accroissement et le développement des investissements dans les
puits privés. En effet, grâce à un puits situé au niveau de l'exploitation agricole, l'agriculteur a une
eau à sa demande, qu'il peut exploiter le jour comme la nuit. De plus, ce dernier est conscient que
le puits va le mettre à 1'abri des contraintes l'approvisionnement en eau de surface. Cependant,
l'installation non contrôlé des puits privés peut entraîner une inégalité à l'accès à l'eau souterrain
puisque seul les agriculteurs qui ont le plus de moyens peuvent l'exploiter.

Utilisation efficient de l’eau d’irrigation

Au niveau des régions où la rareté de l'eau se fait de plus en plus sentir, et que, par suite, les limites
sur les dotations des ressources en eau sont atteintes, l'accent est mis sur le développement d'une
utilisation rationnelle des ressources en eau en particulier sur une utilisation conjuguée efficiente de
l'eau de surface et de l'eau souterraine
La majorité des auteurs reconnaissent que la difficulté majeure à l'achèvement d'une utilisa-
tion efficiente des ressources en eau, est l'interdépendance entre les usagers sur les sources
d'approvisionnement en eau. En effet, les agriculteurs qui dépendent par exemple du même aquifère
ont des coûts 1iés du fait que si quelques uns pompent suffisamment jusqu'à rabattre le niveau de la
nappe au voisinage des puits voisins. Cela va entraîner pour les voisins une augmentation des coûts
de 1'énergie pour le pompage. Si 1'extraction est très importante, jusqu'à induire un tarissement des
puits (et donc une mise hors service), 1'investissement exigé va être encore plus important. Ainsi,
un sur-pompage par quelques agriculteurs engendre une perte pour leurs voisins. Réciproquement,
si une situation de recharge positive conduit à un engorgement et à une salinisation des sols, un
sur-pompage par quelques agriculteurs peut entraîner une diminution des coûts pour les voisins par
l'amé1ioration de la productivité des sols qu'un faible niveau de la nappe va amener. Dans ce cas,
un sur-pompage dans l'aquifère conduit à un bénéfice pour les agriculteurs voisins.
Il apparaît donc que pour une utilisation efficiente de l'eau d'irrigation à 1'échelle régionale, il faut
que les utilisateurs se reconnaissent les effets mutuels sur l'augmentation ou la diminution des coûts
d'exploitation de la ressource en eau.

Ali HAMMANI, 2021

Utilisation conjuguée des eaux de surface et des eaux souterraines 127

Méthodes institutionnelles pour une utilisation conjuguée efficiente

n précédente, longtemps reconnue comme la source potentielle de 1'échec d'une allocation effi-
ciente des ressources en eau de surface et souterraines, plusieurs économistes ont, alors, proposé
diverses solutions d'ordre institutionnel qui ont eu des supports pratiques au niveau de certaines
régions à travers le monde, elles portent sur:
3.3.1.1. La taxation de l'utilisation de l'eau
La solution relative à la taxation de l'utilisation de l'eau souterraine, de manière à ce que les
agriculteurs soient pénalisés ou rémunérés, pour ajuster leurs coûts individuels afin de refléter les
vrais coûts sociaux qui traduisent les impacts négatifs qu'une exploitation locale peut avoir sur les
autres utilisateurs.
3.3.1.2. Centralisation du contrôle
Ce contrô1e peut se faire par le biais d'associations ou de coopératives, qui s'occuperaient de la
gestion des ressources en eau pour le bien de la collectivité.
3.3.1. 3. Affectation de droits de propriétés sur les ressources en eau
Dans ce cas, le système d'allocation des ressources en eau sera basé sur le principe de droits d'eau
reconnus, qui peuvent être sujets à des transactions dans le cadre d'un `` marché de l'eau ''.
Cependant, les alternatives de la taxation du pompage et de la gestion centralisée de l'aquifère
souffrent des problèmes de contrô1e et de suivi tant que les possesseurs de puits sont libres de
choisir quand pomper, et plusieurs peuvent le faire la nuit par exemple. Le même auteur rajoute
que ce n'est pas avec des employés, pour assurer la police de l'eau, que cette situation peut être
redressée. En revanche. La troisième alternative concernant l'affectation de droits sur les ressources,
peut être rendue efficace par l'intermédiaire d'une réforme de la structure 1égale et institutionnelle
relative à l'utilisation de l'eau.
Finalement, une participation des agriculteurs dans le processus de planification et de gestion des
ressources en eau, par 1 'encouragement de leurs organisations sous la forme d'associations ou de
coopératives des usagers de 1'eau d'irrigation, va éliminer ou du moins réduire, l'ampleur des pra-
tiques inefficientes évoquées ci-dessus. En effet, l'agriculteur sera de cette façon, responsabilisé des
conséquences que présente une utilisation non rationnelle de l'eau, l'incitant par la même occasion,
à en faire un usage efficient et équitable tout en encourageant son économie et sa préservation

Méthodes analytiques de l’utilisation conjuguée efficiente

La gestion d'une utilisation conjuguée efficiente est à la fois complexe et techniquement difficile du
fait qu'une gestion efficiente est contrainte par la disponibilité de gestionnaires habiles, de bons
techniciens et par les restrictions institutionnelles qui vont dé1imiter la marge des choix de gestion
possibles. Néanmoins, la disponibilité d'outils scientifiques, pour assister les gestionnaires dans les
processus de prise de décision, est d'un intérêt certain.
La question de base pour la gestion de la nappe concerne: quand et combien peut être extrait
du réservoir souterrain à long terme, sans pour autant, causer des effets indésirables notamment,
l'augmentation des coûts de pompage, ou l'intrusion d'une eau de mauvaise qualité, occasionnés
par des niveaux faibles de la nappe.
Plusieurs modèles de simulation et d'optimisation ont été développés, pour répondre à ses interroga-
tions, si importantes pour les gestionnaires de la ressource en eau souterraine. Ainsi, Gorelick (1988)
a présenté une revue des différentes méthodes de modé1isation de la gestion de l'eau souterraine
(nappe). Il y a deux types de modèles, de gestion de la nappe selon la catégorie de décisions:
3.3.2.1. Modèles de gestion hydraulique
Ce genre de modèle est basé sur un couplage des modèles d'écoulement de surface et des modèles
d'écoulement souterrain (figure ??). Le modèle doit pouvoir simuler les plans d'action qui prennent
en considération les interactions entre la nappe, la zone saturée et la surface du sol.
Les modèles de gestion de la nappe incorporent un modè1e de simulation d'un système particulier
d'aquifère comme contraintes dans le modèle de gestion. Ces derniers, visent la gestion en particulier

Ali HAMMANI, 2021

Utilisation conjuguée des eaux de surface et des eaux souterraines 128

des stress de la nappe tels que le pompage et la recharge, qui sont traités avec les charges
hydrauliques comme des variables de décision du modèle de gestion. La composante modèle de
simulation, est basée sur 1'équation de 1'écoulement souterrain dans les zones saturées.

3.3.2.2. Modèles d'allocation et d'évaluation de politique

Ce type de modèles est valable pour des problèmes complexes lorsque la gestion hydraulique
n'est pas le seul souci du planificateur. Ces modé1es analysent les problèmes d'allocation d'eau et
d'investissement incluant des objectifs de gestion économique et ont été appliques à des problèmes
de grande échelle qui étudient comment une économie agricole va répondre aux politiques institu-
tionnelles et à l'optimisation de l'utilisation conjuguée de l'eau de surface et de l'eau souterraine.

Considération de la qualité de l’eau

pour l’utilisation conjuguée
Une eau convient ou non à 1'irrigation selon la quantité et le type de sels qu'elle contient. Avec
une eau de qualité médiocre, on peut atteindre divers problèmes pédologiques et agronomiques.
Il faut alors mettre en oeuvre des méthodes d'aménagements spéciales afin de maintenir une pleine
productivité agricole.
Il est ainsi évident que pour une meilleure utilisation conjuguée de l'eau de surface et de l'eau
souterraine, la qualité de l'eau revêt la même importance que la quantité avec laquelle cette eau
est disponible. D'autre part, lorsque l'une ou les deux ressources sont salines (ou de mauvaise
qualité), se pose alors le problème de la manière l'utilisation de cette eau pour l'irrigation, sans
porter préjudice aux cultures, et par la suite, affecter négativement les rendements.
Il y a deux alternatives pour l'utilisation de l'eau saline lors d'un programme l'utilisation conjuguée
l'eau de surface et l'eau souterraine pour l'irrigation : (i) l'utilisation séparée de l'eau de surface et
de l'eau souterraine et (ii) le Mélange de l'eau de surface et de l'eau souterraine.

Ali HAMMANI, 2021

Utilisation conjuguée des eaux de surface et des eaux souterraines 129

Utilisation séparée de l’eau de surface et de l’eau souterraine

L'alternative la plus utilisée consiste à irriguer séparément avec une eau de surf ace douce et une
eau souterraine saline ou inversement. Dans ce cas, l'eau douce est utilisée pour permettre la germi-
nation des cultures alors que l'eau saline est utilisée pour les dernières irrigations. En effet, plusieurs
plantes sont très sensibles à la salinité pendant la phase de germination et aux stades précoces de
croissance. Une fois elles sont établies et ont approfondi leur système racinaire, elles peuvent alors
résister à l'effet de l'irrigation avec une eau saline sans pour autant que la concentration des sels
ne dépasse les limites de sensibilité des cultures à la salinité.

Mélange de l’eau de surface et de l’eau souterraine

La seconde alternative pour l'utilisation conjuguée est de mé1anger l'eau saline avec l'eau douce
dans des proportions de telle manière à ce que le mé1ange résultant ne soit pas excessivement salin
restant donc au voisinage des seuils de tolérance des cultures. Dans ce cas, il est possible l'irriguer
continuellement avec l'eau l'irrigation saline mé1angée sans causer une réduction des rendements des
cultures. Les limites du mé1ange de l'eau de surface et de l'eau souterraine vont donc dépendre de
la tolérance des cultures pratiquées à la salinité du mé1ange. Les résultats de ces derniers, donnent
les valeurs du seuil de tolérance des cultures à la salinité définies comme étant la limite à laquelle
aucune perte de rendement ne va se produire à cause de l'utilisation d'une eau salée (tableau …).

Plante Seuil de tolérance Pente de la courbe % Classement

dS/m par dS/m
Orge 8 5 T
Haricot 1 19 S
Maïs 1.7 12 MS
Coton 7.7 5.2 T
Arachide 3.2 29 MS
Riz 3 12 S
Seigle 11.4 10.8 T
Sorgho 6.8 16 MT
Soja 5.0 20.0 MT
Betterave sucre 7 5.9 T
Canne à sucre 1.5 7.9 MS
Blé tendre 6 7.1 MT
Blé dur 5.9 3.8 T
Bersim 1.5 7.5 MS
Tomate 2.5 9.9 MS
Oignon 1.2 16 S
T : tolérant ; MT : moyennement tolérant ; S : sensible ; MS : Moyennement sensible

Le problème de l'utilisation de l'eau saline se pose particulièrement pour l'eau souterraine, du fait
que l'accumulation des sels au niveau de la nappe phréatique, et de la zone racinaire, est une
conséquence inévitable de 1'évapotranspiration des cultures.
On peut utiliser un ratio de conductivité électrique, ECR, qui donne la proportion l'eau souterraine
saline qui peut être mé1angée avec l'eau douce du canal de la façon suivante :

Ali HAMMANI, 2021

Utilisation conjuguée des eaux de surface et des eaux souterraines 130

ECc − ECsw
ECR = (6.1)
ECgw − ECsw

Avec:
• ECc : Seuil de tolérance de la culture A la salinité de l'eau l'irrigation ;
• ECsw : Conductivité électrique de l'eau de surface ;
• ECgw : Conductivité de l'eau souterraine.
Ce ratio ECR exprime le rapport du gain de qualité permis par l'utilisation d'une eau de surface
douce, et de la perte causée par une eau souterraine saline, compte tenu du seuil de tolérance
de la culture à la salinité. De plus, la détermination de ce ratio va permettre d'estimer la quantité
maximale de l'eau souterraine saline (MAXgw ) en terme de volume, qui va être mé1angée avec l'eau
de surface sans causer une réduction des rendements des cultures:

M axgw = IW · ECR (6.2)

ET − RF
IW = IE
(6.3)
IE + LF IE−LF

Avec:
• IW : Quantité totale de l'eau de surface et de l'eau souterraine utilisée pour l'irrigation ;
• ET : Evapotranspiration de la culture durant la période considérée ;
• RF : Précipitation totale au cours de la même période ;
• LF : Dose de fraction de lessivage en % ;
• IE : Efficience de l'irrigation basée sur le % l'eau d'irrigation infiltrée qui est stockée dans la
zone racinaire par pratique normale de l'irrigation.
La quantité IW représente les besoins nets l'irrigation compte tenu de la demande climatique (ET R),
des précipitations, de l'efficience de l'irrigation à la parcelle (pour considérer les doses réelles l'eau
qui arrivent au niveau de la zone racinaire) et de la dose de lessivage qui va dépendre de
la conductivité électrique de l'eau dans la zone racinaire. Cette quantité IW est d'autant plus
grande que la salinité de l'eau au niveau de cette dernière est élevée nécessitant par la suite des
grandes doses pour le lessivage des sels.

Ali HAMMANI, 2021

Vous aimerez peut-être aussi

Notes de Cours Hydrogeologie Automn 2014 ULAVAL PDF
Pas encore d'évaluation
Notes de Cours Hydrogeologie Automn 2014 ULAVAL PDF
284 pages
Ledoux
Pas encore d'évaluation
Ledoux
133 pages
09-Ecoulement Vers Les Ouvrages de Captage
Pas encore d'évaluation
09-Ecoulement Vers Les Ouvrages de Captage
35 pages
Modélisation Numérique en Hydrogéologie
Pas encore d'évaluation
Modélisation Numérique en Hydrogéologie
56 pages
Detay Foragedeau Chapitre10 Biblioannexesindex
Pas encore d'évaluation
Detay Foragedeau Chapitre10 Biblioannexesindex
38 pages
GRM1412 Chapitre1 Hydrodynamique Ecoulement A Surface Libre
Pas encore d'évaluation
GRM1412 Chapitre1 Hydrodynamique Ecoulement A Surface Libre
16 pages
E S Débits Forages
100% (1)
E S Débits Forages
93 pages
Pour L'obtention Du Diplôme de Licence Sciences Et Techniques
Pas encore d'évaluation
Pour L'obtention Du Diplôme de Licence Sciences Et Techniques
55 pages
Cours Modelisation Master 13-10-2021
Pas encore d'évaluation
Cours Modelisation Master 13-10-2021
168 pages
G3 GR Cours Forages VF 111950
Pas encore d'évaluation
G3 GR Cours Forages VF 111950
55 pages
Srat Region de Zinder
Pas encore d'évaluation
Srat Region de Zinder
227 pages
Protection de Oued Fodda
Pas encore d'évaluation
Protection de Oued Fodda
154 pages
Méthodes d'hydrologie appliquées
Pas encore d'évaluation
Méthodes d'hydrologie appliquées
1 page
Rapport de La Modélisation CHEIKH PDF
Pas encore d'évaluation
Rapport de La Modélisation CHEIKH PDF
27 pages
Cours Hydrologie - DR AMARNI Chapitre I
Pas encore d'évaluation
Cours Hydrologie - DR AMARNI Chapitre I
15 pages
Sad Bennis
Pas encore d'évaluation
Sad Bennis
21 pages
Rapport de Stage Agence Bassin Hydraulique 1
Pas encore d'évaluation
Rapport de Stage Agence Bassin Hydraulique 1
19 pages
Cornette Nicolas
Pas encore d'évaluation
Cornette Nicolas
198 pages
Mémoire Hydraulique - Barrage en Terre
100% (1)
Mémoire Hydraulique - Barrage en Terre
77 pages
Analyse Du Fonctionnement Hydro-Sédimentaire Du Bassin Versant Tleta Du Rif Occidental Du Maroc À L'Aide Du Modèle Swat
Pas encore d'évaluation
Analyse Du Fonctionnement Hydro-Sédimentaire Du Bassin Versant Tleta Du Rif Occidental Du Maroc À L'Aide Du Modèle Swat
9 pages
Mobilisation Des Eaux
100% (1)
Mobilisation Des Eaux
2 pages
Cours D'hydraulique A Surface Libre - Chapitre 5 - Ecoulement en Canaux CC1
100% (2)
Cours D'hydraulique A Surface Libre - Chapitre 5 - Ecoulement en Canaux CC1
28 pages
Bonnes pratiques pour réaliser des forages
Pas encore d'évaluation
Bonnes pratiques pour réaliser des forages
48 pages
Essais de Pompage
Pas encore d'évaluation
Essais de Pompage
35 pages
"Eau Cachée du Sahara : Défis et Enjeux"
100% (1)
"Eau Cachée du Sahara : Défis et Enjeux"
17 pages
Creation de Sig Pour Les Reseaux D'eau Du
Pas encore d'évaluation
Creation de Sig Pour Les Reseaux D'eau Du
25 pages
Écoulement Graduellement Varié: Concepts et Méthodes
Pas encore d'évaluation
Écoulement Graduellement Varié: Concepts et Méthodes
26 pages
Cours Hydrodynamique Souterraine 1 M1 Hydro 2023 2024
Pas encore d'évaluation
Cours Hydrodynamique Souterraine 1 M1 Hydro 2023 2024
29 pages
Catalogue BRGM: Rapports Hydrogéologie 1963-1990
Pas encore d'évaluation
Catalogue BRGM: Rapports Hydrogéologie 1963-1990
152 pages
Recharge Artificielle des Nappes Souterraines
Pas encore d'évaluation
Recharge Artificielle des Nappes Souterraines
49 pages
Chap I Introduction
Pas encore d'évaluation
Chap I Introduction
25 pages
Rapport sur l'étude de perméabilité du sol
Pas encore d'évaluation
Rapport sur l'étude de perméabilité du sol
15 pages
TP 6 - SWMM
Pas encore d'évaluation
TP 6 - SWMM
21 pages
Réhabilitation des réseaux d'assainissement
Pas encore d'évaluation
Réhabilitation des réseaux d'assainissement
115 pages
Hydrochimie des aquifères du Cap Bon
100% (1)
Hydrochimie des aquifères du Cap Bon
290 pages
Polycopie de Cours-Hydrologie I
Pas encore d'évaluation
Polycopie de Cours-Hydrologie I
59 pages
Eaux Souterraines La Recharge Maitris e Des Nappes 1655488001
Pas encore d'évaluation
Eaux Souterraines La Recharge Maitris e Des Nappes 1655488001
12 pages
Méthodes de Captage Des Eaux Souterraines
Pas encore d'évaluation
Méthodes de Captage Des Eaux Souterraines
45 pages
Cours Hydrologie Générale - CH 05
Pas encore d'évaluation
Cours Hydrologie Générale - CH 05
9 pages
Hydrogéologie2 - Licence1 PDF
Pas encore d'évaluation
Hydrogéologie2 - Licence1 PDF
78 pages
Auscultation et surveillance des barrages
Pas encore d'évaluation
Auscultation et surveillance des barrages
144 pages
Écoulement en Terrains Fissurés: Analyse
Pas encore d'évaluation
Écoulement en Terrains Fissurés: Analyse
9 pages
Études pour Forage au Dispensaire de Dévégo
Pas encore d'évaluation
Études pour Forage au Dispensaire de Dévégo
23 pages
Cours Hydraulique Generale
Pas encore d'évaluation
Cours Hydraulique Generale
43 pages
Rechemisage d'un forage à La Défense
Pas encore d'évaluation
Rechemisage d'un forage à La Défense
4 pages
Fondamentaux de l'hydraulique et modélisation
Pas encore d'évaluation
Fondamentaux de l'hydraulique et modélisation
43 pages
Mémoire Epanet Messaoudi Guermazi
Pas encore d'évaluation
Mémoire Epanet Messaoudi Guermazi
19 pages
Information: Feuille
Pas encore d'évaluation
Information: Feuille
13 pages
Érosion hydrique des sols en Tunisie
Pas encore d'évaluation
Érosion hydrique des sols en Tunisie
36 pages
Les Inondations en Algerie
100% (1)
Les Inondations en Algerie
1 page
Modélisation Modflow
Pas encore d'évaluation
Modélisation Modflow
172 pages
Essai de Perméabilité À Charge Constante
Pas encore d'évaluation
Essai de Perméabilité À Charge Constante
6 pages
Ahbari
Pas encore d'évaluation
Ahbari
50 pages
Étude de l'Écoulement Souterrain en TP
Pas encore d'évaluation
Étude de l'Écoulement Souterrain en TP
6 pages
Introduction à l'hydrogéologie
Pas encore d'évaluation
Introduction à l'hydrogéologie
30 pages
These
Pas encore d'évaluation
These
131 pages
Modélisation Hydrogéologique: Concepts clés
Pas encore d'évaluation
Modélisation Hydrogéologique: Concepts clés
46 pages
Cours-Hec-Ras PDF
100% (2)
Cours-Hec-Ras PDF
43 pages
Cours Hec Ras PDF
100% (10)
Cours Hec Ras PDF
43 pages
Cours de Modélisation en Hydrogéologie1
Pas encore d'évaluation
Cours de Modélisation en Hydrogéologie1
7 pages
Ecls 1
Pas encore d'évaluation
Ecls 1
166 pages
Bourahouat Nisrine
Pas encore d'évaluation
Bourahouat Nisrine
33 pages
Présentation - PFE ZOUINE Modf
Pas encore d'évaluation
Présentation - PFE ZOUINE Modf
39 pages
Les Étapes de L'étude D'un Projet VRD Sur Covadis
100% (1)
Les Étapes de L'étude D'un Projet VRD Sur Covadis
73 pages
Stage Diagnostic Irrigation et Eau 2024
Pas encore d'évaluation
Stage Diagnostic Irrigation et Eau 2024
1 page
Calibration et validation du modèle SWAT
100% (1)
Calibration et validation du modèle SWAT
13 pages
Limnologie
Pas encore d'évaluation
Limnologie
60 pages
Recomcfbr2017 Dommagesincrementauxdifferentiels
Pas encore d'évaluation
Recomcfbr2017 Dommagesincrementauxdifferentiels
70 pages
Irrigation et Cultures au Nigeria
Pas encore d'évaluation
Irrigation et Cultures au Nigeria
3 pages
Modélisation Hydrologique de la Sisaony
Pas encore d'évaluation
Modélisation Hydrologique de la Sisaony
126 pages
2 Hydrologie
Pas encore d'évaluation
2 Hydrologie
44 pages
Systèmes d'assainissement urbain VRD
Pas encore d'évaluation
Systèmes d'assainissement urbain VRD
16 pages
Programme Prévisionnel Tanger
Pas encore d'évaluation
Programme Prévisionnel Tanger
1 page
Le Mali Hydrographie
Pas encore d'évaluation
Le Mali Hydrographie
1 page
Rapport Groupe 4 PDF
100% (1)
Rapport Groupe 4 PDF
31 pages
Évaluation sur l'eau et ses usages
Pas encore d'évaluation
Évaluation sur l'eau et ses usages
2 pages
Stations d'épuration et pollution
Pas encore d'évaluation
Stations d'épuration et pollution
16 pages
Cours Assainissement Agricole VF
Pas encore d'évaluation
Cours Assainissement Agricole VF
29 pages
UE GEO 121 Cours Prof Dzana Chap1
Pas encore d'évaluation
UE GEO 121 Cours Prof Dzana Chap1
18 pages
Temps de Concentration Dun Bassin Versant
Pas encore d'évaluation
Temps de Concentration Dun Bassin Versant
9 pages
Sodapdf
Pas encore d'évaluation
Sodapdf
81 pages
Vocabulaire essentiel de l'assainissement
Pas encore d'évaluation
Vocabulaire essentiel de l'assainissement
1 page
MÉTHODE-MAC-MATH Intro
Pas encore d'évaluation
MÉTHODE-MAC-MATH Intro
3 pages
Calculer la capacité d'un adoucisseur d'eau
Pas encore d'évaluation
Calculer la capacité d'un adoucisseur d'eau
9 pages
Approvisionnement en Eau Potable RDC
Pas encore d'évaluation
Approvisionnement en Eau Potable RDC
7 pages
2023-2024 Examen de Déchets, Assainissement Et Valorisation
Pas encore d'évaluation
2023-2024 Examen de Déchets, Assainissement Et Valorisation
1 page
Rapport d'expertise sur écoulements d'eau
Pas encore d'évaluation
Rapport d'expertise sur écoulements d'eau
18 pages
Chapitre 9 FC V24 Les Ouvrages Divers
Pas encore d'évaluation
Chapitre 9 FC V24 Les Ouvrages Divers
122 pages
Geologues 167 Eaux Souterraines Ouagadougou
Pas encore d'évaluation
Geologues 167 Eaux Souterraines Ouagadougou
6 pages
Energie Osmotique
Pas encore d'évaluation
Energie Osmotique
1 page
Ressouces en Eau 20241
Pas encore d'évaluation
Ressouces en Eau 20241
2 pages
Dissertation Enjeux Eau
Pas encore d'évaluation
Dissertation Enjeux Eau
1 page
Preliminaires Memoire Tonga Enstp
Pas encore d'évaluation
Preliminaires Memoire Tonga Enstp
14 pages
Calcul du Bilan Hydrique et ETR
Pas encore d'évaluation
Calcul du Bilan Hydrique et ETR
2 pages
Poster de Ressource D'eau Potable
Pas encore d'évaluation
Poster de Ressource D'eau Potable
1 page

Cours GMRES

Transféré par

Cours GMRES

Transféré par

Institut Agronomique et Vétérinaire Hassan II

Département du Génie Rural

3ème année Génie Rural

1 Chapitre 2: Modélisation numérique des écoulements souterrains. . . . . . . . . . . . . . . . . . . . 1

2 Chapitre 2: Modélisation numérique des écoulements souterrains. . . . . . . . . . . . . . . . . . . . 18

Les différents éléments du bilan

Equation générale du bilan

Equation du bilan hydrologique globale d’un bassin versant

Equation du bilan hydrogéologique dans une nappe phréatique

Ali HAMMANI, 2021

• Ip et Ii sont respectivement les infiltrations efficaces des eaux de précipitation et d'irrigation,

Evaluation des différents termes du bi-

7.3.1.1. Evaluation des apports par les précipitations

iii) Calcul de la lame d’eau par la méthode de Theissen

Ali HAMMANI, 2021

en pour-cent de la surface totale du bassin. Ce pourcentage est utilisé comme coefficient de

b. Calcul de la lame d’eau efficace

Ali HAMMANI, 2021

Ali HAMMANI, 2021

7.3.1.2. Evaluation des apports par irrigation

7.3.1.3. Evaluation des apports latéraux

Ali HAMMANI, 2021

Dans cette formule :

7.3.1.4 Les apports par drainance

Où K ′ et B ′ sont, respectivement, la conductivité hydraulique et l'épaisseur de la couche semi-

7.3.2.1. Evaluation des prélèvements par pompages et par drainage artificiel

7.3.2.2. Le drainage naturel et la drainance vers le bas

7.3.2.3. Evaporation et évapotranspiration

∆ = 170 + 8T ± 15cm (1.12)

Ali HAMMANI, 2021

T étant la température moyenne annuelle en o C.

Où : S est le coefficient d'emmagasinement de l'aquifère. Si la remontée n'est pas uniforme et

Exemple de bilan hydrogéologique :

Ali HAMMANI, 2021

[3]: import numpy as np

hmin = h1.min();hmax = h1.max()

xmin, ymin, xmax, ymax=275500, 275000, 335500, 325000

Ali HAMMANI, 2021

# Exécuter la méthode de krigeage ordinare, variogram Gauss

h1, ss1 = OK.execute('grid', grid_x, grid_y)

OK = OrdinaryKriging(x2, y2, h2, variogram_model='linear', verbose=True,␣

h2, ss2 = OK.execute('grid', grid_x, grid_y)

# Résultat de l'interpolation spatiale

contour3 = plt.contour(xintrp, yintrp, dh,20)

Ali HAMMANI, 2021

plt.title('Interpolation spatiale de la charqe hydraulique (%d points)' %␣

Using 'linear' Variogram Model

Calculating statistics on variogram model fit…

Adjusting data for anisotropy…

Using 'linear' Variogram Model

Calculating statistics on variogram model fit…

La variation de la réserve est : -24.400 Mm3

Ali HAMMANI, 2021

[4]: import plotly

go.Surface(z=h2, x=xintrp, y=yintrp, opacity=1.0, colorscale='Blues',␣

#go.Surface(z=dh, x=xintrp, y=yintrp, opacity=1.0, colorscale='Blues',␣

Ali HAMMANI, 2021

Bilan de la nappe de Béni Moussa

interactive(children=(FloatSlider(value=0.1, description='iir', max=0.5, step=0.005),␣

Ali HAMMANI, 2021

[5]: <function __main__.WaterBalance(iir, iip, ip, a, S)>

Ali HAMMANI, 2021

Différents types de modèles

Les modèles physiques

Les modèles analogiques

Il y a deux types d'analogie: analogie entre l'écoulement souterrain et le transfert de chaleur, et

Variable Loi régissant le Equation générale

Transfert de chaleur Température T Loi de Fourrier ∇2 T = 0

Electricité Potentiel V Loi d’Ohm ∇2 V = 0

L'analogie éléctrique est la plus utilisée.

Choix d’une technique numérique et

Ali HAMMANI, 2021

Ali HAMMANI, 2021

Méthode de différences finie en

Ali HAMMANI, 2021

• Φ est le potentielle de charge:

[5]: <function main.WaterBalance(iir, iip, ip, a, S)>