Hydrogéologie L3

François Métivier

► To cite this version:

François Métivier. Hydrogéologie L3. Licence. France. 2025, pp.119. cel-01877908v12

HAL Id: cel-01877908

https://hal.science/cel-01877908v12

Submitted on 13 Jan 2025 (v12), last revised 4 Jan 2026 (v13)

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Hydrogéologie
- L3 -

Version du 13 janvier 2025

François Métivier

Université de Paris, Institut de physique du globe de Paris

<figure><img src="image_1.png" /></figure>

<figure><img src="image_2.png" /></figure>

Université Paris Cité

<figure><img src="image_3.png" /></figure>
# Table des matières

## I Hydrogéologie 1D 6

### 1 Le sol et ses propriétés 7

1.1 Porosité totale 7

1.2 Contenu en eau et saturation 7

1.3 Porosité efficace et capacité de rétention 8

1.4 La nappe 10

1.5 Exercices 11

### 2 Hydrostatique 12

2.1 Pression dans un fluide 12

2.2 Équation de l'hydrostatique 13

2.3 Charge d'un récipient au repos 14

2.4 Piézométrie des nappes 15

2.5 Tension de surface et eau capillaire 17

2.6 Exercices 18

### 3 Écoulement dans un tube, application aux nappes captives 20

3.1 Introduction 20

3.2 Conservation de la masse 20

3.3 Loi de Darcy 22

3.4 Écoulement dans le tube 22

3.5 Application à une nappe captive 23

3.6 Exercices 23

### 4 Écoulement dans une nappe libre 26

4.1 Hypothèse de Dupuit 26

4.2 Etude du problème à une dimension 26

4.3 Terme source et équation de la diffusivité 28

4.4 Exercices 29

### 5 Polluants 31

5.1 Mélange instantané 31

5.2 Mélange classique 32

5.3 Advection et diffusion moléculaire 33

5.4 Équation de dispersion et nombre de Péclet 33

5.5 Solutions en régime permanent 33

5.6 Exercices 36

## II Écoulements dans le plan 37

### 6 Conservation de la masse 38

6.1 Coordonnées cartésiennes 38

6.2 Coordonnées cylindriques 39

6.3 Exercices 41
<table><tbody><tr><td>7</td><td>Equation de Laplace</td><td>43</td></tr><tr><td>7.1</td><td>Expression en coordonnées cartésiennes</td><td>43</td></tr><tr><td>7.2</td><td>Equation de Laplace en coordonnées Polaires et cylindriques</td><td>43</td></tr><tr><td>7.3</td><td>Signification</td><td>43</td></tr><tr><td>7.4</td><td>Conditions aux limites et recherche de solutions</td><td>44</td></tr><tr><td>7.5</td><td>Additivité des solutions</td><td>45</td></tr><tr><td>7.6</td><td>Exercices</td><td>45</td></tr><tr><td>8</td><td>Équipotentielles et lignes de courant</td><td>47</td></tr><tr><td>8.1</td><td>Définitions</td><td>47</td></tr><tr><td>8.2</td><td>Fonction courant</td><td>47</td></tr><tr><td>8.3</td><td>Exemple</td><td>48</td></tr><tr><td>8.4</td><td>Exercices</td><td>49</td></tr><tr><td>9</td><td>Écoulements autour d'un puits</td><td>52</td></tr><tr><td>9.1</td><td>Introduction</td><td>52</td></tr><tr><td>9.2</td><td>Écoulements plan en coordonnées polaire</td><td>52</td></tr><tr><td>9.3</td><td>Écoulement radial et uniforme dans un anneau</td><td>53</td></tr><tr><td>9.4</td><td>Loi de Darcy et débit d'un puits</td><td>55</td></tr><tr><td>9.5</td><td>Rabattement et importance des puits</td><td>56</td></tr><tr><td>9.6</td><td>Exercices</td><td>58</td></tr><tr><td>10</td><td>Écoulement dans une nappe captive</td><td>59</td></tr><tr><td>10.1</td><td>Description du problème et conditions aux limites</td><td>59</td></tr><tr><td>10.2</td><td>Résolution</td><td>60</td></tr><tr><td>10.3</td><td>Débit de la nappe</td><td>61</td></tr><tr><td>11</td><td>Écoulement sous une masse d'eau</td><td>63</td></tr><tr><td>11.1</td><td>Problème(s) à résoudre</td><td>63</td></tr><tr><td>11.2</td><td>Lac alimentant sa nappe</td><td>64</td></tr><tr><td>11.3</td><td>Ajout d'un niveau imperméable</td><td>69</td></tr><tr><td>11.4</td><td>Ajout de charges non nulles sur les bords du domaine</td><td>70</td></tr><tr><td>11.5</td><td>Le problème du Barrage</td><td>72</td></tr><tr><td>11.6</td><td>Exercices</td><td>73</td></tr><tr><td>12</td><td>Problèmes</td><td>75</td></tr><tr><td>12.1</td><td>Stockage de Fontsante</td><td>75</td></tr><tr><td>12.2</td><td>Nappe de Beauce</td><td>76</td></tr><tr><td>12.3</td><td>Plateau de Saclay</td><td>76</td></tr><tr><td>12.4</td><td>Lacs du Médoc</td><td>80</td></tr><tr><td>12.5</td><td>Lacs de Biscarosse</td><td>81</td></tr><tr><td>12.6</td><td>Nappe de Champigny</td><td>81</td></tr><tr><td>III</td><td>Solution partielle des exercices</td><td>87</td></tr></tbody></table>

IV Annexes 103

<table><tbody><tr><td>13</td><td>Précision et incertitudes de mesures</td><td>104</td></tr><tr><td>13.1</td><td>Objectif</td><td>104</td></tr><tr><td>13.2</td><td>Précision des mesures</td><td>104</td></tr><tr><td>13.3</td><td>Variables aléatoires et Incertitude</td><td>106</td></tr><tr><td>13.4</td><td>Érreur systématique</td><td>107</td></tr></tbody></table>

<table><tbody><tr><td>14</td><td>TP 1 : Porosité, capillarité, piézométrie</td><td>108</td></tr><tr><td>14.1</td><td>Dispositif</td><td>108</td></tr><tr><td>14.2</td><td>Capillarité</td><td>108</td></tr><tr><td>14.3</td><td>Porosité</td><td>109</td></tr></tbody></table>
<table><tr><td>15 TP 2 : Expérience de Darcy, Nappe Captive</td><td>111</td></tr><tr><td>15.1 Dispositif</td><td>111</td></tr><tr><td>15.2 Fonctionnement d'une nappe captive</td><td>112</td></tr><tr><td>15.3 Pompage dans une nappe</td><td>113</td></tr><tr><td>16 TP 3 : Nappe libre, transfert de polluant</td><td>115</td></tr><tr><td>16.1 Dispositif</td><td>115</td></tr><tr><td>16.2 Etude d'une nappe libre</td><td>115</td></tr><tr><td>16.3 Propagation d'un polluant</td><td>116</td></tr><tr><td>17 Changelog des versions sur Hal</td><td>118</td></tr></table>
# Cadre du cours

Ce cours<fnref n="1" /> ne décrit ni l'ensemble des écoulements observés en milieu souterrain, ni les méthodes d'étude et de suivi des nappes. Dans ce qui suit je vous propose d'aborder l'étude des écoulements dans un milieu poreux simple. Toutes les équations et développements du cours reposent sur les hypothèses suivantes.

1. Le milieu est **isotrope** : la valeur d'une propriété *P* du milieu, sa perméabilité par exemple, est la même dans les trois directions principales de l'espace. Quand le milieu est anisotrope ses propriétés prennent la forme d'un tenseur.

2. Le milieu est **homogène** : la valeur d'une propriété ne dépend pas de la position.

3. L'eau est **incompressible**, sa densité est constante. Elle ne dépend ni de la température ni de la pression. À l'exception de l'eau de mer, la densité n'est en outre pas modifiée par un mélange avec des éléments ou fluides miscibles.

4. Nous ne considérerons quasiment pas la zone dite non saturée du sol dans lequel les vides sont occupés par un mélange de liquides et de gaz

5. Nous considérerons les milieux dans lesquels se produisent des écoulements souterrains comme des **milieux granulaires** dans lesquels les vides sont tous connectés entre eux.

6. Enfin nous considérerons presque toujours des écoulements stationnaires (ou permanents). Nous n'étudierons pas les variations dans le temps des écoulements sauf dans le cas de modèles simples de pollution.

Malgré toutes ces simplifications nous verrons que les connaissances acquises vous permettront d'aborder au premier ordre et de façon raisonnée de nombreux problèmes d'hydrogéologie. Des références plus avancées ou historiques sont données en bibliographie.

Pour un cours d'hydrogéologie complet de niveau Licence on se reportera à [7] ou [10] ou, mieux encore, au livre de Freeze et Cherry[9]<fnref n="2" />. Pour des développements plus avancés en français [12] reste la référence.

La première partie traite de problèmes d'écoulements unidirectionnels en régime permanent c'est à dire ce qu'il y a de plus simple en hydrogéologie. Elle est adossée à des travaux pratiques afin d'aider les étudiants à percevoir et comprendre ce que sont les écoulements dans le milieu souterrain à partir de l'étude des configurations les plus simples possibles. L'essentiel du cours de L3 repose sur cette première partie. La deuxième partie s'intéresse aux écoulements dans un plan et propose notamment une étude des puits et des solutions de l'équation de Laplace dans différentes configurations. L'équation de Laplace est centrale dans de nombreux domaines de la géophysique, en particulier en hydrogéologie, dans des milieux homogènes. Comprendre cette équation et la façon dont on la résout, de façon analytique et numérique, permet d'acquérir une meilleure compréhension des écoulements souterrains. Enfin la troisième partie propose des solutions partielles de quelques exercices.

<footnote n="1">1. Ce cours est gratuit et distribué sous une Licence Creative Commons Attribution-ShareAlike 3.0 Unported. Il ne vient avec aucune garantie de quelque sorte et ni moi-même, ni mon université, ni mon institut de recherche ne saurions être tenus responsables des effets indésirables de toute mauvaise utilisation ou mauvaise compréhension de ces notes et de leur contenu. Je serais heureux de recevoir vos commentaires et suggestions.</footnote>

<footnote n="2">2. Une version gratuite de ce classique est en ligne au http://hydrogeologistswithoutborders.org/wordpress/original-groundwater-by-freeze-and-cherry-1979-now-available-online/.</footnote>
Première partie

Hydrogéologie 1D
# 1 Le sol et ses propriétés

## 1.1 Porosité totale

Aucun sol ou roche n'est totalement compact. il existe toujours des fissures ou des trous, que nous nommerons indifféremment des **pores**. Ceux-ci peuvent former un réseau dans lequel cohabitent différents fluides (air, eau ou hydrocarbures par exemple).

Soit un volume de sol ou de roche $V_s$. La **porosité totale** d'une roche est la proportion du volume $V_s$ des vides (ce volume peut être occupé par de l'air ou tout autre fluide).

$$ \omega = \frac{V_{vides}}{V_s} \qquad (1.1) $$

La table 1.1 donne des exemples de valeurs de la porosité totale de différents matériaux.

<table><caption>TABLE 1.1 – Porosité totale de quelques types de roches ou de sols</caption><thead><tr><td>Matériau</td><td>Porosité Totale %</td></tr></thead><tbody><tr><td>Quartzite</td><td>0.6</td></tr><tr><td>Granite et gneiss non altérés</td><td>0.02 - 1</td></tr><tr><td>Schistes, ardoises</td><td>0.5 - 7</td></tr><tr><td>Calcaires, Dolomies primaires</td><td>0.5 - 12</td></tr><tr><td>Dolomies secondaires</td><td>10 - 30</td></tr><tr><td>Craie</td><td>8 - 37</td></tr><tr><td>Grès</td><td>3.5 - 38</td></tr><tr><td>Tufs volcaniques</td><td>30 - 40</td></tr><tr><td>Sables</td><td>15 - 48</td></tr><tr><td>Argiles</td><td>44 - 53</td></tr><tr><td>Argiles gonflantes</td><td>jusqu'à 90</td></tr><tr><td>Sols de culture labourés</td><td>45 - 65</td></tr></tbody></table>

## 1.2 Contenu en eau et saturation

Soit un sol de volume $V_s$ et $V_w$ le volume d'eau contenue dans ses pores.

Le contenu volumétrique ou teneur en eau est le rapport entre volume de l'eau et volume du sol soit

Teneur en eau

$$ \theta = \frac{V_w}{V_s}, \quad \theta \in [0, \omega]. \qquad (1.2) $$

$\theta$ peut en théorie varier de 0 la valeur de la porosité $\omega$.

Le degré de saturation est définit comme
Saturation

$$S = \frac{\theta}{\omega}, S \in [0, 1] \qquad (1.3)$$

S ou θ sont utilisés indifféremment pour décrire le contenu en eau du sol.

## 1.3 Porosité efficace et capacité de rétention

### 1.3.1 Définitions

Considérons un sol, ou un matériau granulaire, de volume $V_s$, saturé en eau. On suppose que les pores sont connectés et définissent un réseau permettant la circulation de fluides. Les molécules d'eau situées à proximité des grains interagissent avec eux et forment sur le pourtour des pores une fine couche d'eau dite de rétention<fnref n="1" /> qui ne bouge pas sous l'effet de la gravité. Le volume correspondant à cette couche est appelé le **Volume (d'eau) lié(e)** et noté $V_l$. Cette couche à une épaisseur de l'ordre de grandeur du micron.

<figure><img src="image_4.png" /><figcaption>FIGURE 1.1 – Structure de la couche d'eau entre deux particules. Eau adsorbée et eau pelliculaire forment ce que l'on appelle l'eau de rétention (d'après [12]).</figcaption></figure>

Une circulation des fluides sous l'effet de la gravité ou d'un gradient de pression, n'a lieu que dans la partie des vides où les forces d'interactions moléculaire sont faibles et négligeable, au-delà de la couche de rétention. Nous appelerons le volume dans lequel une telle circulation peut avoir lieu le volume (d'eau) libre (ou gravitaire) et noté $V_g$. Par conservation on a

$$\omega = \frac{V_l + V_g}{V_s} \qquad (1.4)$$

Soit un volume $V_s$ de sol contenant les volumes $V_g$ et $V_l$ d'eau libre et d'eau liée respectivement. La porosité efficace (specific yield in anglais) est définie comme

$$\omega_{eff} = \frac{V_g}{V_s}. \qquad (1.5)$$

La capacité de rétention (retention capacity) du sol est définie comme

$$cr = \frac{V_l}{V_s}. \qquad (1.6)$$

<footnote n="1">1. On distinguera l'eau adsorbée et l'eau pelliculaire cf figure 1.1</footnote>
Il vient alors

$$ \omega = \omega_{eff} + cr. \qquad (1.7) $$

### 1.3.2 Comparaison, notion de réserve utile

La figure 1.2 permet de comparer ces trois grandeurs pour différentes tailles de grains.

<figure><img src="image_5.png" /><figcaption>FIGURE 1.2 – Comparaison entre porosité totale (Porosity), porosité efficace (Specific Yield) et Capacité de rétention (Specific retention) pour des sols de granulométries différentes [11, 12]</figcaption></figure>

La valeur de la porosité efficace dépend de la texture de la roche et de la distribution des tailles des grains minéraux qui la composent. En particulier, la taille des pores est une fonction croissante de la taille des grains ou des minéraux. On peut donc avoir un milieu de grande porosité dont les pores sont petits et dans lesquels l'eau n'est pas libre de s'écouler.

Pour aller un peu plus loin isolons un volume cubique d'un sol composé de grains sphériques uniformes de diamètre d arrangés de façon cubique. isolons un cube de coté d contenant un grains. La porosité totale ω sera donnée par

$$ \omega = 1 - \frac{\pi}{6}. \qquad (1.8) $$

On notera qu'elle est constante et ne dépend pas du diamètre des grains. La capacité de rétention correspond au volume d'épaisseur e entourant le grains. Au premier ordre elle vaut

$$ cr \approx \frac{\pi}{6} \left[ \left( \frac{d+2e}{d} \right)^3 - 1 \right] \qquad (1.9) $$

Dans les faits, le rayon d+e étant plus grand que le coté d du cube, le volume d'eau retenu correspond à la formule précédente moins les 6 calottes dépassant sur chaque face (cf exercice 1.5.3). La porosité efficace se déduit de ces deux grandeurs par

$$ \omega_{eff} = \omega - cr \approx 1 - \frac{\pi}{6} \left( \frac{d+2e}{d} \right)^3 \qquad (1.10) $$

Quand d devient grand, la porosité efficace tend vers la porosité totale et la capacité de rétention diminue. Au contraire quand d diminue et se rapproche de e c'est la porosité efficace qui diminue.

Dans les bilans hydrologiques on recourt souvent à un réservoir que l'on appelle la réserve utile ou encore réserve facilement utilisable du sol. Ce réservoir emmagasine l'eau et diminue jusqu'à sa saturaton la part des eaux qui s'infiltrent d'une part, il participe à l'évapotranspiration d'autre part. D'un point de vue physique cette réserve correspond peu ou prou à la capacité de rétention. Elle dépend donc de la nature du sol, en particulier du diamètre des grains qui le composent, et de sa profondeur.

pour un cube de sol de largeur L aura

$$ RFU \approx \left( \frac{L}{d} \right)^3 \left( d^3 - \frac{\pi}{6} (d+2e)^3 \right) \qquad (1.11) $$
<figure><img src="image_6.png" /><figcaption>FIGURE 1.3 – Nappes libres et captives</figcaption></figure>

Pour sol sableux fin ($d \sim 10^{-4}$ m, $e \sim 5 \cdot 10^{-7}$ m) et un volume unitaire de $1\text{m}^3$ la RFU sera de l'ordre de 0.017 m. Pour un sol limoneux ($d \sim 10^{-4}$ m, $e \sim 5 \cdot 10^{-7}$ m) elle sera de 0.17 m. On est bien dans l'ordre de grandeur de ce qui est décrit dans la littérature.

Dans la suite de ce cours, par souci de simplification, nous ne différencierons pas la porosité efficace de la porosité totale.

## 1.4 La nappe

**Formation aquifère** : Corps de roches perméables contenant une zone suffisamment conductrice pour permettre l'écoulement et le captage de quantités appréciables d'eau

**Formation aquiclude** : Formation réputée imperméable.

**Nappe** : Ensemble des eaux comprises dans la zone saturée d'un ou plusieurs terrains aquifères. La notion de nappe, propre au français, est source d'une ambiguïté constante par rapport à l'anglais qui ne connaît que "groundwater" et "aquifer". En effet le vocabulaire anglais ne reconnaît a priori qu'une nappe = groundwater qui peut s'écouler au travers de plusieurs aquifères alors qu'en France on a tendance à considérer que chaque aquifère saturé correspond à une nappe. C'est la raison pour laquelle par exemple on entendra parler de la nappe du Lutétien, de la nappe du Cuisien (les noms correspondent aux horizons rocheux saturés en eau) etc...Dans ce cours nous considérerons toujours une nappe unique.

**Nappe phréatique** : Étymologiquement nappe dans laquelle on puise. Traditionnellement c'est donc la plus proche de la surface et c'est ainsi qu'est généralement entendu ce terme (les puits profonds sont relativement récents). Il est intéressant de noter qu'en anglais nappe phréatique se traduit par "water table" et donc correspond à ce qu'on appelle la surface libre ou surface piézométrique (cf chapitre 2). Plusieurs formations aquifères
superposées auront la même surface libre dès lors qu'elles sont connectées. En France, souvent, on aura tendance
à les considérer comme des nappes différentes avec chacune leur surface piézométrique.

Nappe libre ou captive : voir Figure 1.3.

1.5 Exercices

1.5.1 Étude d'un sol

Un échantillon cylindrique de 10 cm de long et de 5cm de diamètre a été prélevé dans un sol glaciaire à l'aide d'une tarière. La masse volumique moyenne des minéraux composant le sol est ρ<sub>m</sub> = 2650 kg/m<sup>3</sup> Lors du prélèvement l'échantillon pesait 331.8 g. Après séchage à 105°C l'échantillon ne pèse plus que 302.4 g. Calculez

1. la masse volumique du sol,

2. sa porosité,

3. sa teneur en eau,

4. son pourcentage de saturation.

1.5.2 Porosité

Un sol saturé en eau est composé de grains de diamètre moyen *d* que l'on comptera en microns. L'arrangement des grains est cubique.

1. Quel est le volume $V_p$ et la surface $S_p$ du pore se trouvant à l'intérieur d'un cube de 8 grains de côté $d$ ?

2. La couche d'eau adsorbée contre les grains à une épaisseur d'environ 1µm, quelle est l'ordre de grandeur du volume d'eau liée $V_c$ occupée par cette couche ?

3. Quel est l'ordre de grandeur du volume d'eau libre $V_f$.

4. Quels sont les ordres de grandeur de la porosité totale $\omega_t$, de la porosité efficace $\omega_e$ et de la capacité de rétention cr du pore ?

5. Représentez graphiquement ces grandeurs en fonction du diamètre des grains pour des diamètres variant entre 10µm et 1cm.

1.5.3 porosité efficace et capacité de rétention

<figure><img src="image_7.png" /></figure>

FIGURE 1.4 – Vue en plan d'un cube contenant une sphère. À gauche, le cube et les deux volumes d'eau (libre = porosité efficace, liée = capacité de rétention. À droite schéma montrant les calottes qui dépassent du cube

Un sol est constitué de cubes de coté *d* qui chacun contiennent une sphère de diamètre *d* (empilement dit
cubique). On souhaite calculer la porosité efficace et la capacité de rétention du sol. On fait l'hypothèse que
l'eau emmagasinée par le sol non saturé occupe une couche d'épaisseur *e* autour des grains.

1. Calculer la porosité du volume de sol ainsi défini

2. En considérent un cube unitaire, calculer le volume de la calotte de la sphère de diamètre *d* + *e* qui dépasse de chaque face du cube;

3. En déduire la capacité de rétention et la porosité efficace du cume cubique.
2
Hydrostatique

2.1 Pression dans un fluide

Quand une force est exercée sur une surface le rapport entre cette force et la surface sur laquelle elle s'exerce définit ce que l'on appelle une contrainte (stress en anglais). Cette contrainte peut avoir trois composantes, une composante normale et deux composantes tangentielles à la surface sur laquelle s'applique la force. La composante normale est appelée la pression et notée p.

$$
p dS = \mathbf{F} \cdot \mathbf{n} \tag{2.1}
$$

où **n** est la normale à l'élément de surface sur lequel s'exerce la force **F** et dS la surface de l'élément.

Dans un fluide au repos toutes les contraintes sont normales aux surfaces sur lesquelles elles s'appliquent (dès qu'une contrainte cisaillante s'exerce le fluide s'écoule). Elles sont dirigées vers l'intérieur du fluide.

Considérons l'élément prismatique représenté par schéma de la figure 2.1. On peut écrire un bilan des forces
dans chacune des directions (x et z<fnref n="1" />). Dans la direction x on a

$$
\| \mathbf{F}_3 \| - \| \mathbf{F}_1 \| \sin \beta = 0, \qquad (2.2)
$$

soit

$$
p_3 dz dy - p_1 dy ds \sin \beta = 0. \tag{2.3}
$$

Dans la direction z on a

<figure><img src="image_8.png" /><figcaption>FIGURE 2.1 -</figcaption></figure>

$$
p_2 dx dy - p_1 dy ds \cos \beta - \rho g dV = 0. \tag{2.4}
$$

Des simplifications sont possibles grace aux relations géométriques suivantes :

— l'élément de volume dV s'écrit dV = dxdydz/2,

<footnote n="1">1. La direction des *y* est laissée à votre sagacité.</footnote>
— $ds \cos \beta = dx$,

— $ds \sin \beta = dz$.

On obtient alors

$$p_1 = p_3, \tag{2.5}$$

et

$$p_2 = p_1 + \rho g dz/2. \tag{2.6}$$

Quand $dz \to 0$ $p_2 \to p_1$ et

$$p_1 = p_2 = p_3. \tag{2.7}$$

Conclusion : en un point d'un fluide au repos la pression est la même dans toutes les directions. On vérifie bien ainsi que la pression ne dépend pas de l'orientation d'une surface d'application.

## 2.2 Équation de l'hydrostatique

<figure><img src="image_9.png" /><figcaption>FIGURE 2.2 -</figcaption></figure>

Nous allons effectuer un bilan de forces sur l'élément de volume dont une face est représenté sur le schéma de la figure 2.2. On a dans la direction des $x$ :

$$pdydz - \left( p + \frac{\partial p}{\partial x} \right) dydz = 0; \tag{2.8}$$

dans la direction des $y$ :

$$pdx dz - \left( p + \frac{\partial p}{\partial y} \right) dx dz = 0; \tag{2.9}$$

dans la direction des $z$ :

$$p dxdy - \left( p + \frac{\partial p}{\partial z} \right) dxdy - \rho g dxdydz = 0. \tag{2.10}$$

Après simplifications il vient

$$\frac{\partial p}{\partial x} = \frac{\partial p}{\partial y} = 0, \tag{2.11}$$

et

$$\frac{\partial p}{\partial z} = -\rho g. \tag{2.12}$$

Ces deux équations peuvent être réécrites sous une forme compacte

$$\mathbf{grad}(p + \rho gz) = \mathbf{0} \tag{2.13}$$
<figure><img src="image_10.png" /><figcaption>FIGURE 2.3 -</figcaption></figure>

Enfin, en intégrant (2.13), il vient

$$p + \rho g z = P_0 \qquad (2.14)$$

ou $P_0$ est une constante d'intégration dont la dimension est celle d'une pression. En divisant les deux termes de 2.14 il vient

$$\frac{p}{\rho g} + z = \frac{P_0}{\rho g} \qquad (2.15)$$

L'équation a maintenant la dimension d'une hauteur et le terme de gauche est appelé la **charge** et couramment noté $\phi$

$$\phi = \frac{p}{\rho g} + z \qquad (2.16)$$

L'équation de l'hydrostatique peut alors se réécrire

$$\mathrm{grad}(\phi) = \mathbf{0} \qquad (2.17)$$

## 2.3 Charge d'un récipient au repos

Soit un fluide au repos (figure 2.3), la pression en un point d'altitude $z$ d'un réservoir dont la surface libre se trouve à une altitude $h$ est donnée par (cf (2.13))

$$p = \rho g(h - z) + p_a, \qquad (2.18)$$

$\rho g(h - z)$ représente le poids de la colonne d'eau se trouvant au dessus du point d'altitude $z$, $p_a$ la pression atmosphérique. En négligeant cette dernière et en réarrangeant on arrive à

$$p + \rho g z = \rho g h. \qquad (2.19)$$

Comme $h$ est constante pour un réservoir au repos, (2.19) vérifie donc bien (2.13) et la charge s'écrit

$$\phi = \frac{p}{\rho g} + z = h \qquad (2.20)$$

En tout point d'un réservoir au repos, la **charge** $\phi$ est une constante qui vaut $h$, l'altitude de la surface libre. La charge à la dimension d'une hauteur que l'on nomme aussi hauteur piézométrique ([piezometric] head en anglais). C'est cette hauteur que l'on cherche à mesurer dans les milieux naturels et que l'on utilise pour étudier les écoulements.
## 2.4 Piézométrie des nappes

En hydrogéologie, un piézomètre *stricto sensu* est un trou dans le sol. On peut y mesurer le niveau d'eau si tant est que le trou est suffisamment profond pour atteindre la zone saturée de l'aquifère. Une mare, un puits traditionnel de ferme ou de maison de campagne constituent un piézomètre (fig. 2.4). Le développement des suivis de nappes et les pompages collectifs dans les nappes ont conduit au développement de forages moins larges que les puits traditionnels et dans lesquels on ne se contente pas de creuser un trou. Un piézomètre est aujourd'hui un trou dans laquel on a installé un tube, généralement en PVC, d'un diamètre de l'ordre de la dizaine de cm, tube qui est crépiné (creusé de trous) sur toute ou partie de sa longueur afin de laisser entrer l'eau qui s'écoule dans le sol (fig.2.4). Le trou se rempli donc d'eau, la hauteur du niveau d'eau correspond à ce que l'on appelle la surface piézométrique. Dans le cas où la nappe est libre cette surface correspond à la surface de la nappe.

<figure><img src="image_11.png" /><figcaption>FIGURE 2.4 – Trois puits et une mare à Madagascar</figcaption></figure>

La mesure de la hauteur de la surface piézométrique se fait en deux temps. on mesure tout d'abord la profondeur de la surface de l'eau par rapport au niveau de la surface topographique. Cette mesure se fait soit par des sondes dites manuelles que l'on plonge dans le puits et qui émettent un bruit ou un signal lumineux au contact avec l'eau (fig.2.4), soit par des sondes de pressions qui enregistrent la pression d'eau dans le forage et permettent d'obtenir des enregistrements sur la durée. Cette mesure relative doit ensuite être replacée dans un repère cartésien régional afin de pouvoir comparer les mesures de différents piézomètres. Il faut pour cela connaitre les coordonnées du puits soit au moyen d'un GPS, soit par nivellement topographique classique. Le repère utilisé en France est le niveau géographique de référence (NGF). La figure 2.4 (photo en haut à gauche) montre le géoréférencement d'un puits à l'aide d'un GPS cinématique dans le sud de Madagascar.

Si on dispose d'un réseau de puits ou de piézomètres, il devient alors possible de tracer une carte piézométrique. La figure 2.4 montre une telle carte au niveau du plateau de Saclay. Les piézomètres sont installés dans les sables de Fontainebleau. Il est important de prendre conscience que cette surface piézométrique est loin d'être plane. Dans le cas des sables de Fontainebleau, la surface de la nappe présente une topographie qui mime celle de la topographie réelle. La figure 2.4 montre cette surface en bloc diagramme et permet de s'en convaincre. Les zones blanches correspondent aux vallées de la Bièvre et de l'Yvette que la nappe alimente.
<figure><img src="image_12.png" /><figcaption>FIGURE 2.5 – Piézométrie. À gauche et de haut en bas : schéma d'un piézomètre, courbe piézométrique d'un piézomètre de la nappe de Beauce. À droite et de haut en bas : piézomètre en place, mesure piézométrique, sonde de mesure. Sources : wikimédia pour les images, base de données ADES (https://ades.cauefrance.fr/) pour la courbe piézométrique.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIGURE 2.6 – Carte piézométrique de la nappe des sables de Fontainebleau au niveau du plateau de Saclay</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 2.7 – Bloc diagramme de la surface piézométrique des sables de Fontainebleau au niveau du plateau de Saclay</figcaption></figure>

## 2.5 Tension de surface et eau capillaire

Si on verse de l'eau sur une surface plane elle s'étale. Si le volume versé est petit on constate que l'étalement n'est pas infini. le volume s'étale jusqu'à former une goutte plus ou moins sphérique (fig. 2.5).

De la même façon, si on plonge un tube de faible diamètre dans un récipient contenant de l'eau, on constate que l'eau monte dans le tube jusqu'à une certaine hauteur supérieure au niveau d'eau dans le récipient (fig 2.5). On constate en outre que l'interface entre l'eau et l'air dans le tube est courbe.

Ces phénomènes sont liés à l'existence d'une surface de contact entre deux fluides, dans le cas qui nous préoccupe l'eau et l'air. La grandeur qui caractérise ces effets s'appelle le coefficient de tension de surface $\gamma$. Elle a la dimension d'une énergie par unité de surface ou d'une force par unité de longueur (Nm$^{-1}$). Le coefficient de tension de surface d'un interface eau/air est $\gamma = 70 \cdot 10^{-3} Nm^{-1}$ à 20°C

Deux relations nous intéressent dans le cadre de ce cours (la démonstration de ces deux lois fait l'objet d'un exercice). La première donne la relation existant entre la différence de pression $\Delta p = p_{in} - p_{out}$ entre l'intérieur et l'extérieur d'une goutte de fluide de rayon $R$ et le coefficient de tension de surface. C'est la loi de Laplace.
<figure><img src="image_15.png" /><figcaption>FIGURE 2.8 – Étalement de gouttes sur une surface. le volume des gouttes augmente de droite à gauche</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 2.9 – Montée d'eau colorée dans un tube capillaire de diamètre 1mm.</figcaption></figure>

Loi de Laplace

$$ \Delta p = \frac{2\gamma}{R} \qquad (2.21) $$

La seconde relation donne la relation existant entre la hauteur $h$ de montée de l'eau dans un tube capillaire de rayon $r$ et le coefficient de tension de surface. Cette relation, la loi de Jurin, dépend en outre de l'angle du ménisque que fait l'eau contre la paroi du tube. Dans le cas d'un interface eau/air/verre l'angle est de l'ordre de $\theta \approx 60^{\circ}$.

Loi de Jurin

$$ h = \frac{2\gamma \cos \theta}{\rho g r} \qquad (2.22) $$

Dans le milieu souterrain les pores ont des diamètres variables mais qui sont couramment de l'ordre de la centaine de microns voire moins. Ces pores jouent donc le rôle de capillaires et, en vertu de la loi de Jurin, l'eau monte dans les pores. La hauteur de montée dans un capillaire est de l'ordre de 1cm pour un diamètre de 1mm. la loi de jurin nous montre donc que pour un diamètre de pores de 100microns, la remontée sera d'environ 1m au dessus du niveau piézométrique qui, en raison du diamètre des piézomètres ($\sim 0.1m$), sera très peu sensible aux effets de tension de surface. La zone saturée d'une nappe contient donc ce que l'on appelle une frange capillaire qui varie de quelques cm à plusieurs mètres en fonction de la taille des pores. Nous ne nous en préoccuperons pas dans la suite de ce cours mais il est utile de savoir que cette frange existe.

## 2.6 Exercices

### 2.6.1 Piézométrie

1. Une réservoir contient deux fluides immiscibles de l'eau et un hydrocarbure de densité 1.2. La surface libre du fluide de l'eau est située à une hauteur $h_1$, l'interface entre les deux fluides à une hauteur $h_2$. Deux piézomètres sont en contact l'un avec l'eau, l'autre avec l'hydrocarbure.
<figure><img src="image_17.png" /><figcaption>FIGURE 2.10 – Biseau salé</figcaption></figure>

(a) Quelles sont les charges (hauteurs de fluide) dans les deux piézomètres ?

(b) Quelle est la pression à la base du récipient ?

2. Répondre aux même questions si on considère un récipient remplit de trois liquides non miscibles de densités 0.7,1 et 1.6 dont les interfaces sont respectivement situées aux altitudes $h_1, h_2$ et $h_3$ ?

## 2.6.2 Tension de surface

1. Démontrer la loi de Laplace. Pour ce faire on peut calculer le travail des forces de pression et de tension de surface si on fait varier une goutte d'eau de rayon $r$ d'une quantité infinitésimale $dr$.

2. Démontrer la loi de Jurin. Pour ce faire on peut calculer la pression sous la surface de l'eau dans le capillaire de deux façons différentes.

3. Utiliser la loi de Jurin pour estimer la hauteur de montée capillaire en fonction du diamètre des pores.

4. Représenter graphiquement le résultat pour des diamètres de pores allant de 1 micron à 1cm.

## 2.6.3 Nappe côtière, biseau salé

On considère une nappe côtière de surface libre (figure 2.10). On cherche à creuser un puits afin d'alimenter une station balnéaire. Le puits serait creusé à une distance $L$ du rivage. La hauteur piézométrique à cet endroit est $h$. On cherche la profondeur maximale à laquelle on peut forer sans atteindre la limite eau douce – eau salée. On suppose que l'interface entre les deux types d'eau est abrupte et stationnaire.

1. Écrire en chaque point P de l'interface la pression correspondant à l'eau de mer et à l'eau salée.

2. Calculer la profondeur Z de l'interface sous un point où la surface piézométrique en surface est à hauteur $h$ sachant que la masse volumique de l'eau douce est 1000 kg/m3 et celle de l'eau salée 1025 kg/m3 (cela correspond à 32 grammes de sel par litre d'eau). Cette relation est connue sous le nom de « Principe de Ghyben-Herzberg ». Noter que l'on a supposé ici que le biseau est une droite.

3. Application : $h = 2$ m, $L = 200$ m, calculer la profondeur de la limite eau douce – eau salée. Quel est l'angle du biseau ?

4. On aménage un puits dont la base est à l'altitude –20m. Calculer le rabattement maximal admissible.
# Ecoulement dans un tube, application aux nappes captives

## 3.1 Introduction

L'objectif de ce chapitre est en premier lieu de comprendre à quoi peut ressembler l'écoulement dans une expérience comme celle que l'on voit dans la figure 3.1. Ceci nous amènera à définir ce qu'est la loi dite de Darcy et à appréhender le fonctionnement d'une nappe captive.

<figure><img src="image_18.png" /><figcaption>FIGURE 3.1 – Schéma d'une expérience de filtration horizontale. Dans l'expérience originiale de Darcy la filtra- tion était verticale.</figcaption></figure>

## 3.2 Conservation de la masse

La figure 3.2 schématise l'écoulement dans le tube horizontal de la figure 3.1. Cet écoulement est caractérisé par sa vitesse $U$ horizontale. Le flux de masse traversant la surface $dS$ dans le sens du vecteur normal à la surface $dS(x)$ s'écrit

$$F(x) = \rho U(x) dS(x) = \rho U(x) i(-dS i) = -\rho U(x) dS \quad (3.1)$$

Il est négatif. Le vecteur étant dirigé vers l'extérieur à la surface cela veut dire que le flux circule dans le sens opposé à la normale $dS$. Il entre dans le volume

Le flux traversant la surface $dS(x + dx)$ dans le sens indiqué par le vecteur surface s'écrit

$$F(x + dx) = \rho U(x + dx) dS(x + dx) = \rho U(x + dx) i(dS i) = \rho U(x + dx) dS \quad (3.2)$$
<figure><img src="image_19.png" /><figcaption>FIGURE 3.2 -</figcaption></figure>

Il est positif. Le vecteur étant dirigé vers l'extérieur à la surface cela veut dire que le flux circule dans le sens de la normale $dS$. Il sort du volume. L'origine de cette «inversion» des signes (entrant négatif, sortant positif) provient du fait que le vecteur normal à un élément de surface englobant un volume fermé est dirigé vers l'extérieur du volume.

La somme des flux donne

$$ \Sigma = F(x + dx) + F(x) = \rho(U(x + dx) - U(x))dS \quad (3.3) $$

Elle correspond à un flux sortant si elle est positive (la masse dans l'élément de volume diminue), entrant si elle est négative (la masse dans l'élément de volume augmente). Durant un intervalle de temps $dt$ la variation de masse $dM$ dans l'élément de volume sera donc

$$ (F(x + dx) + F(x))dt = -dM \quad (3.4) $$

ou encore

$$ \frac{dM}{dt} = -(F(x + dx) + F(x)) \quad (3.5) $$

Si la somme des flux est positive (il sort plus de masse qu'il n'en rentre) la variation de masse est négative ($dM/dt < 0$). Si la somme des flux est négative (il rentre plus de masse qu'il n'en sort) la variation de masse est positive ($dM/dt > 0$).

Si on souhaite écrire ce bilan de façon vectoriel on aura

$$ \frac{dM}{dt} = -\rho(\mathbf{U}(x)d\mathbf{S}(x) + \mathbf{U}(x+dx)d\mathbf{S}(x+dx)) \quad (3.6) $$

ou encore

$$ \frac{dM}{dt} = - \sum_{S} \rho \mathbf{U} d\mathbf{S} \quad (3.7) $$

Reprenons notre bilan (3.5) et remplaçons les flux par leurs expressions (3.1) et (3.2) on a

$$ \frac{dM}{dt} = \rho(U(x) - U(x+dx))dS \quad (3.8) $$

En effectuant un développement limité au premier ordre de $U(x + dx)$ on a

$$ \frac{dM}{dt} = \frac{d\rho}{dt}xdS = -\rho\frac{dU}{dx}xdS \quad (3.9) $$

soit

$$ \frac{d\rho}{dt} = -\rho \frac{dU}{dx} \quad (3.10) $$

Dans le cadre de ce cours<fnref n="1" /> $\rho$ est constante. on aboutit donc à l'équation de conservation

$$ \frac{dU}{dx} = 0 \quad (3.11) $$

dont la solution est

$$ U = Cte. \quad (3.12) $$

<footnote n="1">1. L'approximation $\rho = Cte$ reste valable dans la plupart des problèmes d'hydrogéologie.</footnote>
## 3.3 Loi de Darcy

L'expérience de Darcy [6] est publiée en 1856. Son principe est représenté figure 3.1. Deux réservoirs remplis d'eau dont la hauteur est maintenue constante sont reliés par un milieu poreux. La charge $\phi$ est constante dans chacun des deux réservoirs. Elle vaut $\phi = h_1$ dans le premier réservoir et $\phi = h_2$ dans le second. Cette propriété reste valable à la limite de chaque milieu poreux si bien qu'on peut écrire

$$ \phi(x_1) = h_1, \phi(x_2) = h_2 \qquad (3.13) $$

La loi de Darcy obtenue à partir de cette expérience s'écrit

$$ Q = KA \frac{\phi(x_1) - \phi(x_2)}{x_2 - x_1} \qquad (3.14) $$

ou A est la section de la colonne sableuse, et K une constante dépendant du milieu que l'on nomme la conductivité hydraulique ou encore perméabilité des hydrogéologues (dimension $LT^{-1}$)

On peut définir une vitesse à partir de (3.14) :

$$ U = \frac{Q}{A} = K \left( \frac{\phi(x_1) - \phi(x_2)}{x_2 - x_1} \right). \qquad (3.15) $$

C'est la vitesse de Darcy ou vitesse de filtration. C'est une vitesse fictive. Elle correspond à la vitesse moyenne qu'aurait l'écoulement si celui-ci s'effectuait au travers de toute la section (milieu solide inclus donc). On peut approximer la vitesse réelle du fluide en tenant compte de la porosité du milieu $\omega$ par l'expression

$$ v = \frac{U}{\omega} = \frac{K}{\omega} \left( \frac{\phi(x_1) - \phi(x_2)}{x_2 - x_1} \right). \qquad (3.16) $$

Dans les équations (3.14), (3.15) et (3.16) qui précèdent on constate que $(\phi(x_1) - \phi(x_2))/(x_2 - x_1)$ est sans dimension et représente une pente (figure 3.1a) ou un gradient (figure 3.1b). Si on réduit progressivement la distance $x_2 - x_1$ dans le cas de l'expérience horizontale on arrivera à

$$ \lim_{x_2 \to x_1} \frac{\phi(x_1) - \phi(x_2)}{x_2 - x_1} = -\frac{d\phi}{dx}. \qquad (3.17) $$

À partir de ce que nous venons de voir, l'équation (3.15) donnant la valeur algébrique de la vitesse de Darcy dans l'expérience éponyme s'écrit

$$ U = -K \frac{d\phi}{dx} \qquad (3.18) $$

Afin de généraliser nous remarquerons que la vitesse va des charges les plus fortes vers les charges les moins fortes. Le vecteur vitesse est donc opposé au vecteur gradient de charge soit in fine

Loi de Darcy

$$ \mathbf{U} = -K \frac{d\phi}{dx} \qquad (3.19) $$

## 3.4 Écoulement dans le tube

Si maintenant nous couplons loi de Darcy (3.20) et conservation de la masse (3.11) nous obtenons une équation différentielle du second ordre

Equation de Laplace 1D

$$ \frac{d^2\phi}{dx^2} = 0. \qquad (3.20) $$

Considérons notre expérience de Darcy (figure 3.1a). L'échantillon de milieu poreux est contenu dans un tube de section constante et très étroit dans le sens perpendiculaire au schéma (On ne considérera donc qu'un
écoulement unidimensionnel dans la direction du tube). Dans le réservoir de gauche la charge vaut $h_1$ et dans celui de droite elle vaut $h_2$. Les faces de l'échantillon sont donc aux charge $\Phi(x_1) = h_1$ et $\Phi(x_2) = h_2$.

Notre domaine d'étude est donc un tube de courant unidimensionnel limité par deux faces de charge constante. Le problème a résoudre s'écrit

Problème aux limites de l'expérience de Darcy (figure 3.4a)

$$ \frac{d^2\phi}{dx^2} = 0 \quad (EDP) \qquad (3.21) $$

$$ \phi(x_1) = h_1 \quad (CL1) \qquad (3.22) $$

$$ \phi(x_2) = h_2 \quad (CL2) \qquad (3.23) $$

Nous sommes dans un cas très simple où une solution par intégration directe de (3.21) est possible. On obtient

$$ \phi(x) = ax + b. \qquad (3.24) $$

(3.24) représente la solution générale de l'équation de Laplace pour un écoulement unidimensionnel. On constate que cette équation présente deux constantes d'intégrations $(a, b)$. La donnée des deux conditions aux limites (3.22) et (3.23) est donc nécessaire et suffisante pour définir une valeur de ces deux constantes et obtenir une solution unique à notre problème.

On obtient alors l'expression de la charge dans le tube

$$ \phi = \frac{h_2 - h_1}{x_2 - x_1}(x - x_1) + h1. \qquad (3.25) $$

Notons qu'on retrouve bien

1. que la charge varie de façon linéaire entre les deux réservoirs,

2. que la vitesse est constante et horizontale.

## 3.5 Application à une nappe captive

Cette solution a une application importante. En effet, la démonstration en sera faite plus tard, cet écoulement correspond à l'écoulement dans une nappe captive telle que celle décrite dans la figure 10.1 quand les charges aux extrêmités de la nappe sont constantes. Le tube de l'expérience de Darcy peut être remplacé par un écoulement entre deux plaques espacés d'une distance $w$. Dès lors que l'écoulement est confiné la solution que nous venons d'obtenir s'applique.

## 3.6 Exercices

### 3.6.1 Écoulements simple (1D)

1. De l'eau coule dans une colonne de sable (longueur $L = 120cm$, surface $A = 200cm^2$, $\Delta h = 120cm$). Le sable a une porosité $n = 0.36$ et une perméabilité $K = 20m/j$ quelle sont

(a) le débit au travers de la surface,

(b) la vitesse de Darcy,

(c) la vitesse moyenne de l'eau dans la colonne ?

2. Un traceur radioactif indique une vitesse moyenne de $V = 0.75m/j$ dans un aquifère. la pente de la surface piézométrique est $dh/dx = 0.002$. Déterminer la perméabilité $K$ de l'aquifère si sa porosité est de $n = 0.2$.

### 3.6.2 Vitesse de filtration et vitesse réelle dans un milieu constitué de grains carrés !

on considère une section de roche de côté $A$. Cette section est remplie partiellement de grains (carrés !) de côté $a$. la vitesse moyenne dans les pores est $u$. Elle est orientée perpendiculairement à la section. Dans le cas représenté sur la figure 3.4 calculez la vitesse de filtration $U$ c'est à dire la vitesse moyennée sur la surface de la section. Cette vitesse est aussi connue sous le nom de vitesse de Darcy.
<figure><img src="image_20.png" /><figcaption>FIGURE 3.3 – Section dans un échantillon de roche</figcaption></figure>

### 3.6.3 Vitesse de filtration et vitesse réelle dans un milieu constitué de grains ronds !

on considère une section de roche de côté 10. Cette section est remplie partiellement de grains diamètre 1. la vitesse moyenne dans les pores est $u$. Elle est orientée perpendiculairement à la section. Dans le cas représenté sur la figure 3.4 calculez la vitesse de filtration $U$ (ou vitesse de Darcy) c'est à dire la vitesse moyennée sur la surface de la section.
<figure><img src="image_21.png" /><figcaption>FIGURE 3.4 – Section dans un échantillon de roche</figcaption></figure>
# 4 Écoulement dans une nappe libre

## 4.1 Hypothèse de Dupuit

Par définition une nappe libre présente une surface matérielle, libre d'évoluer dans le temps et l'espace, qui correspond à l'entrée dans la zone saturée du sol où $S = 1$. Les problèmes à surface libre sont plus complexes à résoudre que les problèmes sans surface libre comme la nappe captive que nous venons d'étudier. La figure 4.1a donne une idée de ce à quoi ressemblent les écoulements dans cette situation. En régime stationnaire, la surface libre est une ligne de courant, c'est à dire qu'elle est tangente au vecteur vitesse, et les flux au travers de celle-ci sont nuls. la base de la nappe est elle aussi une ligne de courant, simple car horizontale comme toujours dans ce cours. Stricto sensu la direction et le sens du vecteur vitesse changent donc dans la nappe. Heureusement pour nous, la pente de la surface libre est souvent très petite. On peut alors faire une simplification importante connues sous le nom d'approximation de Dupuit.

Soit $\theta$ l'angle que fait la surface libre avec l'horizontale en un point. La vitesse en ce point s'écrit

$$ \mathbf{U} = \|\mathbf{U}\|\cos\theta\mathbf{i} + \|\mathbf{U}\|\sin\theta\mathbf{k} \quad (4.1) $$

Si $\theta$ est très petit alors $\sin\theta \sim 0$, $\cos\theta \sim 1$ d'où

$$ \mathbf{U} \approx \|\mathbf{U}\|\mathbf{i} \quad (4.2) $$

C'est l'hypothèse de Dupuit. Dans le cas d'un écoulement dans l'espace, l'approximation de Dupuit conduit à

$$ \mathbf{U} \approx U_x \mathbf{i} + U_y \mathbf{j}. \quad (4.3) $$

Dans tous les cas, l'approximation de Dupuit revient donc à négliger la composante verticale de la vitesse : $U_z \approx 0$.

## 4.2 Etude du problème à une dimension

Considérons un écoulement entre deux piézomètres de charge $h_1$ et $h_2$ identiques aux cas de la nappe captive étudié précédemment

### 4.2.1 Problème à résoudre

La vitesse étant horizontale, la charge est verticale et correspond à la hauteur de la surface libre $\phi = h$. On obtient alors

$$ \mathbf{U} = -K \frac{d\phi}{dx} \mathbf{i} = -K \frac{dh}{dx} \mathbf{i}. \quad (4.4) $$

Le débit par unité de largeur de la nappe $q$ est alors obtenu en multipliant la norme de cette vitesse par la hauteur de la nappe

$$ q = \mathbf{U}\mathbf{S} = \mathbf{U}h\mathbf{i} = -Kh\frac{dh}{dx}, \quad (4.5) $$

soit

$$ q = -\frac{K}{2} \frac{dh^2}{dx} \quad (4.6) $$

Si maintenant nous effectuons un bilan de masse sur un élément de volume comme dans le cas de la nappe captive nous arrivons à
<figure><img src="image_22.png" /><figcaption>(a) Écoulement complet</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(b) Écoulement approximé après application de l'hypothèse de Dupuit</figcaption></figure>

FIGURE 4.1 – Écoulement dans une nappe libre. La forme de la surface libre est volontairement exagérée afin de bien comprendre l'intérêt et les liites de l'approximation de Dupuit.

### Problème au limites de la nappe libre 1D

$$ \frac{d^2 h^2}{dx^2} = 0, \qquad (4.7) $$

$$ h(x_1) = h_1, \qquad (4.8) $$

$$ h(x_2) = h_2. \qquad (4.9) $$

## 4.2.2 Solution

L'équation (4.7), non linéaire, diffère de l'expression de l'équation de Laplace 1D à cause de l'approximation de Dupuit. Elle s'intègre cependant simplement. On obtient ainsi la solution du problème :

### Solution du problème au limites de la nappe libre 1D

$$ \phi = h = \sqrt{\frac{h_1^2 - h_2^2}{x_1 - x_2}(x - x_1) + h_1^2}. \qquad (4.10) $$
<figure><img src="image_24.png" /><figcaption>FIGURE 4.2 – Comparaison des profils piézométriques d'une nappe libre (sous l'approximation de Dupuit) et d'une nappe captive. Les conditions aux limites sont identiques.</figcaption></figure>

### 4.2.3 Commentaires

[1] L'approximation de Dupuit est largement vérifiée dans la nature. Elle signifie que l'écoulement dans une nappe libre a essentiellement lieu dans le plan horizontal. Cette conséquence est très importante car c'est elle qui nous permettra de traçer les lignes d'écoulement et les équipotentialles d'une nappe libre à partir de mesure piézométrique de puits.

[2] Dans les deux situations, confinée ou libre, ce que nous venons de faire revient à considérer que la pression dans la nappe est (quasi) hydrostatique. De ce fait la charge est constante sur l'épaisseur de la nappe et correspond à la hauteur de la surface piézométrique. Dans le cas d'une nappe libre cette surface correspond on outre à la surface libre de la nappe.

## 4.3 Terme source et équation de la diffusivité

### 4.3.1 Bilan

La recharge des nappes libres se fait par les précipitations. Il est donc intéressant d'ajouter un terme source à l'équation de conservation de la masse afin de tenir compte de cette recharge.

Soit $P$ une précipitation moyenne ($([L][T]^{-1})$ la conservation de la masse appliquée à une colonne de nappe (figure) donne

$$U(x + dx)h(x + dx) - U(x)h(x) - P = 0. \qquad (4.11)$$

En développant en série limitée au premier ordre le produit $U(x + dx)h(x + dx)$ puis en simplifiant on arrive à

$$\frac{\partial U h}{\partial x} - P = 0. \qquad (4.12)$$

L'application de la loi de Darcy nous permet alors d'obtenir

$$\frac{\partial^2 h^2}{\partial x^2} + \frac{2P}{K} = 0. \qquad (4.13)$$

Cette équation est connue sous le nom d'équation de la diffusivité. Elle est présentée (comme ailleurs dans l'ensemble du cours) dans sa version en régime permanent.

### 4.3.2 Problème du plateau 1D

On imagine un plateau de largeur $L$ séparé par deux rivières de charge $h_1$ et $h_2$ et recevant une précipitation moyenne $P$. On néglige l'écoulement en $y$ en considérant le plateau de profondeur infinie. Le problème à résoudre est alors
<figure><img src="image_25.png" /><figcaption>FIGURE 4.3 – Écoulement d'une nappe libre sous un plateau entre deux rivières.</figcaption></figure>

### Problème aux limites de la recharge d'un plateau

$$ \frac{\partial^2 h^2}{\partial x^2} + \frac{2P}{K} = 0 \qquad (4.14) $$

$$ h(0) = h_1 \qquad (4.15) $$

$$ h(L) = h_2 \qquad (4.16) $$

La solution de ce problème s'obtient par intégration directe de (4.14) et conduit à (exercice !)

$$ h^2 = -\frac{P}{K}x^2 + \left( \frac{h_2^2 - h_1^2}{L} + \frac{PL}{K} \right) x + h_1^2 \qquad (4.17) $$

La figure 4.3 montre le résultat pour un écoulement sans pluie et avec pluie. On note l'apparition d'une convexité de la surface libre. Le sommet de la surface libre correspond à un changement de sens du gradient de charge et donc à un changement de sens de l'écoulement. On voit apparaître une limite de partage des eaux. Les précipitations alimentent les deux rivières.

## 4.4 Exercices

### 4.4.1 Débit d'une expérience de nappe libre

Votre expérience de nappe libre fonctionne bien à l'exception du débitmètre qui ne marche plus ! On rappelle les paramètres de votre expérience :

- Longueur de la nappe : $L = 1 \text{ m}$;

- Largeur de la nappe : $W = 0.02 \text{ m}$;

- Hauteur piézométrique en entrée $h_1(x = 0) = 0.1 \text{ m}$;

- Hauteur piézométrique en sortie : $h_2(x = 1) = 0.01 \text{ m}$;

- Conductivité hydraulique $K = 10^{-3} \text{ m/s}$

1. Quelle est la valeur entière du débit par unité de largeur en ml/mn/m ?

2. Quelle est la valeur entière du débit que vous devriez mesurer au débitmètre ?
<figure><img src="image_26.png" /><figcaption>FIGURE 4.4 – Plateau séparant deux rivières.</figcaption></figure>

4.4.2 Limite de partage des eaux

On considère un plateau séparant deux rivières (fig. 4.4. La distance entre les deux cours d'eau est L = 1000 m. Les altitudes des surfaces libres des rivières 1 et 2 sont respectivement de h₁ = 10 m, et h₂ = 20 m. La précipitation annuelle qui s'infiltre est P = 0.5 m/a. La conductivité du milieu poreux est K = 10⁻⁵ m/s. On suppose que l'approximation de Dupuit est vérifiée.

1. Comment définit-on la limite de partage des eaux dans une nappe libre ?

2. Quelle est la position de cette limite ?

3. Quelle est l'altitude maximale de la nappe ?

4. Sur quelle partie de la nappe l'approximation de Dupuit est-elle effectivement valide ?
# Dispersion et transfert de polluants

## 5.1 Mélange instantané

La situation la plus simple lors de l'étude d'une pollution consiste à considérer un réservoir comme un volume dans lequel les phénomènes de mélanges sont instantanés. En d'autre terme la concentration au sein du réservoir est constante.

Soit un réservoir de volume $V$. Ce réservoir est alimenté en amont par une rivière et se déverse en aval dans une autre rivière. Soit $Q_i$ flux d'eau constant à l'entrée dur réservoir, $Q_o$ flux d'eau constant à la sortie, $V_0$ volume initial du réservoir. Le volume du réservoir à un instant $t$ est donné par

$$V = V_0 + (Q_i - Q_o)t \qquad (5.1)$$

On injecte un traceur dans la rivière alimentant le réservoir (en amont donc). Le flux de masse de traceur entrant dans le réservoir est $C_i Q_i$. Soit M la masse de traceur dans le réservoir, l'instantanéité du mélange implique que la concentration en sortie est $C = M/V$ et le flux sortant $M_o = (Q_o M)/V$.

La variation de masse du traceur dans le temps<fnref n="1" /> est égale à la différence entre flux de masse entrant et flux sortant soit

$$\dot{M} = C_i Q_i - \frac{Q_o M}{V_o + (Q_i - Q_o)t} \qquad (5.2)$$

Soit

$$\dot{M} + \left( \frac{Q_o}{V_o + (Q_i - Q_o)t} \right) M = C_i Q_i \qquad (5.3)$$

### 5.1.1 Cas d'une masse d'eau de volume constant

Si $Q_i = Q_o = Q$, le volume $V$ de la masse d'eau reste constant et (5.3) devient

$$\dot{M} + \left(\frac{Q}{V}\right) M = C_i Q \qquad (5.4)$$

$C_i = 0$ : on part d'une masse $M_0$ et on étudie sa diminution par l'exutoire (pollution ponctuelle par exemple) la solution de (5.4) est obtenue par séparation de variables et donne

$$M = M_0 \exp \left( - \frac{Qt}{V} \right) \qquad (5.5)$$

Le temps caractéristique de rajeunissement de l'eau est le rapport entre débit sortant et masse du réservoir.

$$\frac{dX}{dt} = \dot{X} = d_t X, \quad \frac{d^2 X}{dt^2} = \ddot{X} = d_{tt}^2 X$$

<footnote n="1">1. On notera</footnote>

<footnote>De façon usuelle on utilise les points pour les dérivées temporelles et les primes pour les dérivées spatiales.</footnote>
$C_i \neq 0$ : afin de résoudre (5.3) On peut effectuer un changement de variable

$$Z = M - C_i V \qquad (5.6)$$

(5.4) devient alors

$$\dot{Z} + \frac{Q}{V} Z = 0 \qquad (5.7)$$

Dont la solution en Z est de la forme de (5.5). En revenant à M on obtient

$$M(t) = C_i V \left[ 1 - \exp \left( -\frac{Qt}{V} \right) \right] + M_0 \exp \left( -\frac{Qt}{V} \right) \qquad (5.8)$$

Cette équation très générale sera utilisée dans de nombreux autres cas comme par exemple pour vous décrire l'équilibre dynamique d'une chaîne de montagnes.

### 5.1.2 Puits et sources

En plus de ce qui précède on peut avoir des apports ou soustractions au sein du réservoir. Si un élément décroît naturellement dans le milieu comme un élément radioactif on aura un ou plusieurs termes supplémentaires qui peuvent ou pas être constants. Par exemple si l'on considère un élément qui décroît de façon proportionnelle à sa concentration (5.3) devient

$$\dot{M} + \left( \frac{Q_o + K}{V_o + (Q_i - Q_o)t} \right) M = C_i Q_i \qquad (5.9)$$

## 5.2 Mélange classique

Dans de nombreuses situations, notamment en ce qui concerne les masses d'eau souterraines, le mélange au sein d'une masse d'eau n'est pas instantané. Il faut donc effectuer une conservation de la masse incluant une dimension spatiale.

Soit un volume élémentaire de fluide $dV$ traversé par un flux $F$ d'un élément soluble au travers de sa surface $dS$. La concentration en élément soluble $C$ est supposée petite devant 1, ce qui veut dire qu'une variation de concentration dans l'élément de volume se fait à volume constant.

La conservation de la masse impose que la variation de quantité de l'élément dans le volume au cours du temps est égale à la somme du flux sur sa surface

$$\frac{\partial}{\partial t} \iiint_V C dV = - \iint_S F \cdot dS. \qquad (5.10)$$

Le signe négatif vient du fait que la normale à la surface est dirigée vers l'extérieur. Prenons l'exemple d'un flux entrant. Comme le produit scalaire d'un flux dirigé vers l'intérieur avec une normale dirigée vers l'extérieur est négatif $-F \cdot dS$ sera donc positif et correspondra bien à un flux entrant donc à une augmentation de la concentration de l'élément. Par application du théorème de Gauss

$$\iiint_V \mathrm{div} F dV = \iint_S F \cdot dS, \qquad (5.11)$$

(5.10) devient

$$\iiint_V \left( \frac{\partial C}{\partial t} + \mathrm{div} F \right) dV = 0, \qquad (5.12)$$

Ou encore

$$\frac{\partial C}{\partial t} + \mathrm{div} F = 0. \qquad (5.13)$$

On se limitera ici à des cas à une dimension. Dans de telles situations (5.13) se simplifie en

$$\frac{\partial C}{\partial t} + \frac{\partial F}{\partial x} = 0. \qquad (5.14)$$
## 5.3 Advection et diffusion moléculaire

Les densités de flux F sont de deux types :

1. Advectif, les molécules sont transportées par le courant :

$$F = CU \qquad (5.15)$$

2. Diffusif, les molécules se déplacent sous l'effet du gradient de concentration ou par un effet de dispersion hydrodynamique dans le milieu poreux (c'est souvent ce dernier effet qui prédomine dans les milieux poreux [13]):

$$F = -D\nabla C \qquad (5.16)$$

Dans le cas le plus simple la diffusivité D est constante. En sommant ces deux types de flux et en les intégrant à l'équation de conservation de la masse, on obtient alors

$$\frac{\partial C}{\partial t} + U \frac{\partial C}{\partial x} - D \frac{\partial^2 C}{\partial x^2} = 0. \qquad (5.17)$$

## 5.4 Équation de dispersion et nombre de Péclet

Si on considère un tube de courant de longueur L on peut calculer deux temps caractéristiques

1. un temps d'advection : $\tau_a = L/U$

2. un temps de diffusion : $\tau_d = L^2/D$

le rapport des deux est sans dimension et s'appelle le nombre de Péclet

$$Pe = \frac{\tau_d}{\tau_a} = \frac{UL}{D} \qquad (5.18)$$

On distingue trois cas

1. Pe >> 1 la diffusion est très lente devant l'advection et est donc négligeable. (5.17) se simplifie en

$$\frac{\partial C}{\partial t} + U \frac{\partial C}{\partial x} = 0 \qquad (5.19)$$

2. Pe << 1 l'advection est très lente devant la diffusion. (5.17) se simplifie en

$$\frac{\partial C}{\partial t} - D \frac{\partial^2 C}{\partial x^2} = 0 \qquad (5.20)$$

3. Enfin Pe ~ 1 on ne peut négliger aucun des flux. Il faut résoudre (5.17).

## 5.5 Solutions en régime permanent

### 5.5.1 Pe >> 1

on se place dans le cas le plus simple d'une pollution constante de concentration $C(x = 0) = C_0$ l'équation d'advection devient

$$U \frac{\partial C}{\partial x} = 0, \qquad (5.21)$$

soit

$$C = C_0. \qquad (5.22)$$

La concentration en polluant est constante dans toute la nappe
<figure><img src="image_27.png" /><figcaption>FIGURE 5.1 – Solutions en régime permanent de l'équation de dispersion en fonction du nombre de Péclet.</figcaption></figure>

### 5.5.2 $Pe \ll 1$

On est dans le cadre d'une diffusion et

$$D \frac{\partial^2 C}{\partial x^2} = 0 \qquad (5.23)$$

L'équation différentielle est du second degré donc ont doit disposer de deux conditions aux limites pour la résoudre. Ce qu'on voit bien en écrivant la solution générale

$$C = ax + b \qquad (5.24)$$

On se place dans le cas le plus simple d'une pollution constante de concentration $C(x = 0) = C_0$ et d'une sortie dans une rivière de concentration nulle en polluant $C(x = L) = 0$ et située à une distance $L$ de la source. La solution est alors de la forme

$$C(x) = C_0 \left(1 - \frac{x}{L}\right). \qquad (5.25)$$

### 5.5.3 $Pe \sim 1$

Enfin dans le dernier cas on a

$$U \frac{\partial C}{\partial x} - D \frac{\partial^2 C}{\partial x^2} = 0 \qquad (5.26)$$

Comme l'équation est de degré deux on a deux conditions au limite par exemple $C(0) = C_0$ et $C(L) = C_L$.

Dans le cadre de ce cours $U$ et $D$ sont constantes soit

$$\frac{\partial}{\partial x} \left( C - \frac{D}{U} \frac{\partial C}{\partial x} \right) = 0. \qquad (5.27)$$

ou encore

$$C - \frac{D}{U} \frac{\partial C}{\partial x} = B. \qquad (5.28)$$

La solution du terme est de la forme

$$C = A \exp \left( \frac{Ux}{D} \right) + B \qquad (5.29)$$

$$A = \frac{C_0 - C_L}{1 - \exp(UL/D)} \qquad (5.30)$$

$$B = \frac{C_L - C_0 \exp(UL/D)}{1 - \exp(UL/D)} \qquad (5.31)$$

Les solutions des trois cas sont représentées sur la figure 5.1
<figure><img src="image_28.png" /><figcaption>FIGURE 5.2 – bilans des flux hydriques en $10^3$ km$^3$/an (d'après [15]).</figcaption></figure>

<table><caption>TABLE 5.1 – Masse des différents réservoirs d'eau sur la terre</caption><thead><tr><td>Réservoir</td><td>Masse (10<sup>15</sup>Kg)</td><td>% de la masse</td></tr></thead><tbody><tr><td>Océans</td><td>1400000</td><td>95.9</td></tr><tr><td>Glaces et Neiges</td><td>43400</td><td>2.97</td></tr><tr><td>Subsurface (nappes)</td><td>15300</td><td>1</td></tr><tr><td>Surface (rivières et lacs)</td><td>360</td><td>0.02</td></tr><tr><td>Atmosphère</td><td>15.5</td><td>0.001</td></tr><tr><td>Biosphère</td><td>2</td><td>0.0001</td></tr><tr><td>Total</td><td>1459077</td><td>100</td></tr></tbody></table>
## 5.6 Exercices

### Bilan global et temps de résidence

1. À partir de la figure 5.2 calculer la hauteur d'eau précipitée et évaporée à la surface du globe. Le bilan est il bouclé ?

2. Calculer ces mêmes grandeurs sur les continents. Qu'observez-vous ?

3. Quelle(s) hypothèse(s) pouvez-vous alors formuler ?

4. À partir du volume des réservoirs (table 5.1) et des valeurs de flux (figure 5.2) estimer le temps de résidence de l'eau dans ces différents réservoirs (les regrouper si nécessaire).

5. À quoi correspond ce temps ?

6. Que deviendrait le temps de résidence des continents si le réservoir des glaces disparaissait ?

7. Quelle serait l'évolution de la masse d'eau des continents en fonction du temps si, dû à la fonte des glaciers, le ruissellement était multiplié par deux.

### Pollution ponctuelle d'un réservoir

Un lac artificiel circulaire présente un rayon $R = 100\text{m}$ et une profondeur moyenne $H = 2\text{m}$. A la suite d'une pollution accidentelle, un volume $V = 10000\text{l}$ de saumure contenant une masse $M = 2000\text{kg}$ d'un sel est déversé dans le lac. Le lac est traversé par une petite rivière dont le débit est $Q = 0.1 \text{ m}^3/\text{s}$ et dont la concentration en ce sel est nulle. Trouver la quantité de sel présente dans le lac à tout instant $t$.

### Pollution d'un lac : régimes stationnaire et transitoire

Un lac limité en phosphore d'une surface de $A = 80 \cdot 10^6 \text{ m}^2$ est alimenté par une rivière de débit $Q = 15 \text{ m}^3/\text{s}$ dont la concentration en phosphore est $C_i = 0.01 \text{ mg/l}$. Un effluent ponctuel, issu d'une station de traitement, ajoute un débit massique $Q_M = 1 \text{ g/s}$ de phosphore dans la rivière, en amont du lac. La vitesse de sédimentation du phosphore est $v_s = 10 \text{ m/an}$. La rivière ressort du lac avec le même débit $Q$ et une concentration $C$ inconnue.

1. Estimer la concentration totale en phosphore dans le lac et la rivière en sortie quand le régime stationnaire est atteint.

2. Quelle devrait le flux de masse de l'effluent pour que la concentration du lac reste égale à celle de la rivière.

3. La profondeur moyenne du lac est $H = 10 \text{ m}$ et la concentration en phosphore initiale du lac est nulle. Trouver la quantité de sel présente dans le réservoir à tout instant $t$.

(a) si l'on néglige la sédimentation,

(b) si l'on tient compte de la sédimentation.

4. Quatre ans après la mise en place de l'effluent, un glissement de terrain réduit le débit à l'exutoire du lac à $7 \text{ m}^3/\text{s}$. Quelle est l'évolution de la concentration dans le lac et de la profondeur d'eau du lac dans chacun des deux cas précédents<fnref n="2" />.

<footnote n="2">2. Cette question pour être résolue demande que vous connaissiez l'usage des facteurs intégrants.</footnote>
Deuxième partie

Écoulements dans le plan
# 6 Conservation de la masse

## 6.1 Conservation de la masse dans un élément de volume

Dans le cas d'un écoulement tridimensionnel on a

$$ \mathbf{U} = U_x(x, y, z)\mathbf{i} + U_y(x, y, z)\mathbf{j} + U_z(x, y, z)\mathbf{k} \qquad (6.1) $$

le flux de masse au travers de la surface $d\mathbf{S}_x$ sera donc

$$ F_x(x, y, z) = \rho U_x(x, y, z) dS_x \qquad (6.2) $$

De la même façon on aura

$$ F_x(x + dx, y, z) = \rho U_x(x + dx, y, z) dS_x \qquad (6.3) $$

En effectuant la somme de ces deux flux on aboutit à

$$ F_x(x) + F_x(x + dx) = \rho \frac{\partial U_x}{\partial x} dS_x dx = \rho \frac{\partial U_x}{\partial x} dy dz dx \qquad (6.4) $$

Si on fait de même sur les quatre autres faces on obtient

$$ F_y(y) + F_y(y + dy) = \rho \frac{\partial U_y}{\partial y} dxdzdy \qquad (6.5) $$

et

$$ F_z(z) + F_z(z + dz) = \rho \frac{\partial U_z}{\partial z} dxdydz \qquad (6.6) $$

La somme des trois sommes de flux pendant un temps dt est égale à $-dM$ soit

$$ \frac{dM}{dt} = \frac{d\rho}{dt} dxdydz = -\rho \left( \frac{\partial U_x}{\partial x} + \frac{\partial U_y}{\partial y} + \frac{\partial U_z}{\partial z} \right) dxdydz \qquad (6.7) $$

ou encore

$$ \frac{d\rho}{dt} = -\rho \operatorname{div} \mathbf{U}, \qquad (6.8) $$

où l'opérateur divergence div est définit par

### Divergence

$$ \operatorname{div} = \mathbf{i} \frac{\partial}{\partial x} + \mathbf{j} \frac{\partial}{\partial y} + \mathbf{k} \frac{\partial}{\partial z}. \qquad (6.9) $$

Dans le cadre de ce cours cette équation se réduit à

### Conservation de la masse

$$ \operatorname{div} \mathbf{U} = 0. \qquad (6.10) $$
<figure><img src="image_29.png" /><figcaption>FIGURE 6.1 -</figcaption></figure>

## 6.2 Conservation de la masse en coordonnées cylindriques

Il est important de connaître cette formulation car, en hydrogéologie, un des écoulements les plus importants est celui qui se produit autour d'un puits. Or écrire la conservation de la masse autour d'un puits revient à chercher l'écriture de

$$ \mathrm{div} \mathbf{U} = 0 \qquad (6.11) $$

en coordonnées cylindriques. En effet on peut facilement imaginer qu'autour d'un puits l'écoulement sera orienté de façon radiale et que dans ces conditions il est plus intéressant de travailler en coordonnées cylindriques ou polaires (en particulier dans le cas d'un écoulement plan)

### 6.2.1 Gradient d'une fonction

Les coordonnées cartésiennes et cylindriques sont liées par les relations suivantes

$$ x = r \cos \theta \qquad (6.12) $$

$$ y = r \sin \theta \qquad (6.13) $$

$$ z = z \qquad (6.14) $$

de la même façon les vecteurs unitaires **i**, **j**, **k** définissant un repère cartésien sont reliés aux vecteurs unitaires **e**<sub>r</sub>, **e**<sub>θ</sub>, **e**<sub>z</sub> définissant un repère cylindrique par

$$ \mathbf{i} = \mathbf{e}_r \cos \theta - \mathbf{e}_\theta \sin \theta \qquad (6.15) $$

$$ \mathbf{j} = \mathbf{e}_r \sin \theta + \mathbf{e}_\theta \cos \theta \qquad (6.16) $$

$$ \mathbf{k} = \mathbf{e}_z \qquad (6.17) $$

Le gradient d'une fonction f(x,y,z) s'écrit

$$ \mathrm{grad} f = \frac{\partial f}{\partial x} \mathbf{i} + \frac{\partial f}{\partial y} \mathbf{j} + \frac{\partial f}{\partial z} \mathbf{k} \qquad (6.18) $$
Afin d'obtenir le gradient d'une fonction en coordonnées cylindrique on va chercher le gradient de $g(r, \theta, z) = f(r \cos \theta, r \sin \theta, z)$ on écrira

$$ \frac{\partial f}{\partial r} = \frac{\partial f}{\partial x} \frac{\partial x}{\partial r} + \frac{\partial f}{\partial y} \frac{\partial y}{\partial r}, \qquad (6.19) $$

soit

$$ \frac{\partial f}{\partial r} = \frac{\partial f}{\partial x} \cos \theta + \frac{\partial f}{\partial y} \sin \theta. \qquad (6.20) $$

En effectuant le même calcul pour la dérivée partiel de $f$ par rapport à $\theta$ on obtient

$$ \frac{\partial f}{\partial \theta} = -r \frac{\partial f}{\partial x} \sin \theta + r \frac{\partial f}{\partial y} \cos \theta. \qquad (6.21) $$

Enfin pour $z$ rien ne change.

A partir de (6.20) et (6.21) on peut alors inverser les relations et écrire

$$ \frac{\partial f}{\partial x} = \frac{\partial f}{\partial r} \cos \theta - \frac{1}{r} \frac{\partial f}{\partial \theta} \sin \theta, \qquad (6.22) $$

et

$$ \frac{\partial f}{\partial y} = \frac{\partial f}{\partial r} \sin \theta + \frac{1}{r} \frac{\partial f}{\partial \theta} \cos \theta. \qquad (6.23) $$

En reprenant chaque terme de (6.18) et en le remplaçant par son expression en coordonnées cylindriques on obtient alors

$$ \begin{aligned} \mathrm{grad} f = & \left( \frac{\partial f}{\partial r} \cos \theta - \frac{1}{r} \frac{\partial f}{\partial \theta} \sin \theta \right) (\mathbf{e}_r \cos \theta - \mathbf{e}_\theta \sin \theta) \\ & + \left( \frac{\partial f}{\partial r} \sin \theta + \frac{1}{r} \frac{\partial f}{\partial \theta} \cos \theta \right) (\mathbf{e}_r \sin \theta + \mathbf{e}_\theta \cos \theta) \\ & + \frac{\partial f}{\partial z} \mathbf{e}_z, \end{aligned} \qquad (6.24) $$

qui se simplifie après calcul et donne

$$ \mathrm{grad} f = \nabla f = \frac{\partial f}{\partial r} \mathbf{e}_r + \frac{1}{r} \frac{\partial f}{\partial \theta} \mathbf{e}_\theta + \frac{\partial f}{\partial z} \mathbf{e}_z, \qquad (6.25) $$

où l'opérateur gradient ($\nabla$) en coordonnées cylindrique est défini par

$$ \nabla f = \mathbf{e}_r \frac{\partial}{\partial r} + \frac{\mathbf{e}_\theta}{r} \frac{\partial}{\partial \theta} + \mathbf{e}_z \frac{\partial}{\partial z}. \qquad (6.26) $$

### 6.2.2 Divergence d'un vecteur

Si maintenant on considère une vitesse définie par

$$ \mathbf{U} = U_r \mathbf{e}_r + U_\theta \mathbf{e}_\theta + U_z \mathbf{e}_z \qquad (6.27) $$

la divergence de cette vitesse s'obtiendra par le produit scalaire

$$ \mathrm{div} \mathbf{U} = \nabla \cdot \mathbf{U}. \qquad (6.28) $$

On développe ce produit qui donne

$$ \begin{aligned} \mathrm{div} \mathbf{U} = & \mathbf{e}_r \frac{\partial}{\partial r} (U_r \mathbf{e}_r + U_\theta \mathbf{e}_\theta + U_z \mathbf{e}_z) \\ & + \frac{\mathbf{e}_\theta}{r} \frac{\partial}{\partial \theta} (U_r \mathbf{e}_r + U_\theta \mathbf{e}_\theta + U_z \mathbf{e}_z) \\ & + \mathbf{e}_z \frac{\partial}{\partial z} (U_r \mathbf{e}_r + U_\theta \mathbf{e}_\theta + U_z \mathbf{e}_z) \end{aligned} \qquad (6.29) $$

le premier et le troisième terme se simplifient aisément. La petite difficulté vient de la dérivation angulaire. En effet les dérivées des vecteurs unitaires $\mathbf{e}_r$ et $\mathbf{e}_\theta$ par rapport à $\theta$ ne sont pas nulles. Elles correspondent à une rotation de $\pi/2$ des vecteurs unitaires et valent donc respectivement

$$ \frac{\partial \mathbf{e}_r}{\partial \theta} = \mathbf{e}_\theta, \quad \frac{\partial \mathbf{e}_\theta}{\partial \theta} = -\mathbf{e}_r, \qquad (6.30) $$
Sachant cela on peut alors simplifier (6.30) et il vient

$$ \mathrm{div} \mathbf{U} = \frac{\partial U_r}{\partial r} + \frac{U_r}{r} + \frac{1}{r} \frac{\partial U_{\theta}}{\partial \theta} + \frac{\partial U_z}{\partial z} = 0. \qquad (6.31) $$

## 6.3 Exercices

### 6.3.1 Variation de vitesse dans une nappe

la figure 6.2 représente une section dans une nappe alluviale. Les mesures d'altitude de la surface libre de la nappe, mesurées dans deux puits séparés de 2km sont $z_1 = 100m$ et $z_2 = 90m$. On fait l'hypothèse que l'altitude de la surface libre varie de façon linéaire entre ces deux points. La hauteur de la nappe au droit du premier puits est de $h_1 = 35m$ et divers sondages indiquent que la base de la nappe est horizontale. On a de plus estimé par traçage la vitesse de l'eau dans la roche au point 1 ($u_1 = 10^{-4} m/s$). Enfin diverses mesures sur échantillons révèlent que la porosité efficace de la roche est constante dans l'aquifère et vaut $\omega_{eff} \sim 15\%$

<figure><img src="image_30.png" /><figcaption>FIGURE 6.2 – Coupe dans une nappe libre</figcaption></figure>

1. Quelles est la vitesse de filtration dans la nappe au point 1 ?

2. Quel est le débit par unité de largeur dans la nappe ?

3. Quelle sont la vitesse de filtration et la vitesse réelle au point 2 ?

### 6.3.2 Champs de vecteurs et champs de vitesse

On considère les champs de vecteurs suivants

$$ x\mathbf{i} - y\mathbf{j} \qquad (6.32) $$

$$ x^2\mathbf{i} - y^2\mathbf{j} \qquad (6.33) $$

$$ y\mathbf{i} + x\mathbf{j} \qquad (6.34) $$

$$ y\mathbf{i} - x\mathbf{j} \qquad (6.35) $$

$$ y^2\mathbf{i} + x^2\mathbf{j} \qquad (6.36) $$

$$ y^2\mathbf{i} + x^2\mathbf{j} \qquad (6.37) $$

$$ (x + y)\mathbf{i} + (x + y)\mathbf{j} \qquad (6.38) $$
$$2xyi - y^2j \qquad (6.39)$$

Ces champs peuvent-ils correspondre au champ de vitesses d'un écoulement incompressible ?
# 7 Equation de Laplace

## 7.1 Expression en coordonnées cartésiennes

On a vu que la conservation de la masse d'un fluide incompressible sur une élément de volume aboutissait à

$$ \nabla \cdot \mathbf{u} = 0 \qquad (7.1) $$

En couplant cette équation à l'équation de Darcy (3.20) il vient

### Équation de Laplace

$$ \nabla^2 \phi = \Delta \phi = 0 \qquad (7.2) $$

C'est l'équation de Laplace décrivant l'écoulement d'un fluide incompressible dans un milieu isotrope, homogène et lui même incompressible.

En coordonnées cartésiennes l'équation de Laplace s'écrit

$$ \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} + \frac{\partial^2 \phi}{\partial z^2} = 0 \qquad (7.3) $$

## 7.2 Equation de Laplace en coordonnées Polaires et cylindriques

En coordonnées cylindriques on aura (cette équation se déduit de (6.26))

$$ \frac{\partial^2 \phi}{\partial r^2} + \frac{1}{r} \frac{\partial \phi}{\partial r} + \frac{1}{r^2} \frac{\partial^2 \phi}{\partial \theta^2} + \frac{\partial^2 \phi}{\partial z^2} = 0 \qquad (7.4) $$

On peut bien sûr écrire l'équation de Laplace en coordonnées sphériques mais, si cette forme est très importante en géodésie elle n'a aucun intérêt en hydrogéologie.

## 7.3 Signification

Prenons une courbe définie par l'équation $y = f(x)$. La dérivée première de $f$ donne la tangente au point $x$, la dérivée seconde correspond à la courbure de $f$. Trois cas se représentent

1. $d^2f/dx^2 > 0$ : la courbure est dirigée vers le haut. Si on trace un segment reliant deux points $(x_1, x_2)$ entourant $x$ et équidistants de $x$, ce segment passe au dessus de la courbe et $(y_1 + y_2)/2 > y$. La valeur de $y$ au point $x$ est plus petite que la moyenne des valeurs des points entourant $x$.

2. $d^2f/dx^2 < 0$ : la courbure est dirigée vers le bas. Si on trace le même segment il passe cette fois au dessous de la courbe et $(y_1 + y_2)/2 < y$. La valeur de $y$ au point $x$ est plus grande que la moyenne des valeurs des points entourant $x$.

3. $d^2f/dx^2 = 0$ : la courbure est nulle. Le segment passe par $x$ et $(y_1 + y_2)/2 = y$. La valeur de $y$ au point $x$ est la moyenne des valeurs des points entourant $x$.
<figure><img src="image_31.png" /><figcaption>FIGURE 7.1 – Signification de l'équation de Laplace</figcaption></figure>

L'équation de Laplace signifie donc que la valeur de la charge en un point est égale à la moyenne des charges avoisinantes. Nous nous servirons de cette propriété quand nous aborderons les solutions numériques de l'équation de Laplace.

## 7.4 Conditions aux limites et recherche de solutions

### D'un ensemble de solutions à une solution unique

L'équation de Laplace admet une infinité de solutions comme par exemple les équations suivantes

$$ \phi = (ax + b)(cy + d) \qquad (7.5) $$

$$ \phi = a \ln(r) + b \qquad (7.6) $$

$$ \phi = \sin(x) \cosh(z) \qquad (7.7) $$

$$ \phi = r^n \sin(n\theta) \qquad (7.8) $$

Par contre quand le problème est suffisamment bien contraint il n'existe plus qu'une solution. Il est donc très important lors d'une étude de cas de bien définir les conditions dans lesquelles a lieu l'écoulement. Pour ce faire il faut définir le domaine dans lequel l'écoulement se déroule ainsi que les conditions aux limites de ce domaine.

### 7.4.1 Les différentes types de conditions au limites

Le nombre de conditions aux limites nécessaires pour obtenir une solution d'une équation aux dérivées partielles dépend à la fois du nombre de variables et de l'ordre de l'équation (le degré de la dérivation). Pour chaque dimension on aura besoin d'autant d'un nombre de conditions au limites égal au degré de l'équation différentielle pour cette dimension. dans notre exemple : une dimension $x$ et une équation du second ordre donne $1 \times 2 = 2$ conditions au limites pour résoudre le problème. La résolution de l'équation de Laplace dans l'espace implique donc la connaissance de 6 conditions aux limites deux pour chaque dimension de l'espace. Si on a moins de CL le problème est sous-déterminé et on aura une infinité de solutions, si on en a plus le problème est sur-déterminé et soit une condition est inutile soit le problème n'admet pas de solution.

Il existe trois types classiques de conditions aux limites d'un domaine d'étude.

#### Condition de potentiel imposé (Dirichlet)

$$ \Phi = f(x, y, z, t) \qquad (7.9) $$

Un cas courant en hydrogéologie et qui correspond à notre exemple précédent est $\Phi = Cte$.

#### Condition de flux imposé (Neumann)

$$ u_n = \nabla \Phi \times 1n = \partial_n \Phi = f(x, y, z, t) \qquad (7.10) $$

Ici $n$ indique la direction normale à la surface limite et selon laquelle se fait l'écoulement. Dans le cas d'une limite imperméable $\partial_n \Phi = 0$ (en 2D la limite imperméable est une ligne de courant cf infra).
## Limite mixte ou de Cauchy

$$ \partial_n \Phi + \lambda(x, y, z) \Phi = f(x, y, z) \qquad (7.11) $$

Cette condition s'applique à des écoulements au travers d'une limite semi-perméable par exemple.

**Cas particulier d'une surface libre** Dans le cas d'une surface libre (sans apports) $\Psi = Cte$. la surface fonctionne comme une limite imperméable et correspond à une ligne de courant.

### 7.4.2 Solution analytique versus numérique

Quand elles existent les solutions analytiques d'un problème bien posé présentent l'avantage de donner à voir la façon dont les différentes variables, paramètres et conditions aux limites d'une équation différentielle sont liés entre eux. Elle permet donc de comprendre le rôle de chacun dans la forme d'un écoulement. En outre elles donnent accès à la connaissance exacte de l'écoulement en tout point du plan ou de l'espace. Dans la mesure des moyens disponibles elles sont donc à privilégier. Cependant la recherche de solutions analytiques devient rapidement ardue dans deux situations :

1. quand les conditions aux limites définies précédemment se complexifient, quand la charge ou les flux ne sont plus des constantes par exemple;

2. quand le domaine de l'écoulement devient complexe et ne se réduit pas à un demi plan, un carré ou un cercle.

Les problèmes aux limites n'ont alors rapidement plus de solution analytiques ou, s'ils en ont, leur résolution demande une éducation en mathématiques (notamment l'usage des transformations conforme et des solutions dans le plan complexe) hors de portée de la plupart des étudiants de géosciences. Heureusement pour nous l'équation de Laplace se prête à des solutions numériques que l'on peut mettre en œuvre de façon relativement simple. Posséder une solution analytique permet en outre de tester la méthode de résolution, la vitesse de convergence de l'algorithme etc...Nous utiliserons les différents études d'écoulements présentées dans les chapitres suivants pour présenter des techniques très simples en python.

## 7.5 Additivité des solutions

L'équation de Laplace présente une propriété propre aux applications linéaires : les solutions sont additives, c'est à dire que la somme de deux solutions de l'équation de Laplace est une solution de l'équation de Laplace. Cette dernière est définie sur le domaine commun aux deux premières solutions et caractérisé par la somme des conditions aux limites. Cette propriété très importante est utilisée pour décomposer un problème en plusieurs sous problèmes plus simples comme le montre la figure 7.2.

## 7.6 Exercices

### 7.6.1 Additivité

Montrez que la somme de deux fonctions solutions de l'équation de Laplace est solution de l'équation de Laplace.

### 7.6.2 Orthogonalité

Calculer les intégrales suivantes où $(m, n) \in \mathbb{N}^2$ :

$$ I_1 = \int_0^1 \sin(m\pi x) \sin(n\pi x) \qquad (7.12) $$

$$ I_2 = \int_0^1 \cos(m\pi x) \cos(n\pi x) \qquad (7.13) $$

$$ I_3 = \int_0^1 \sin(m\pi x) \cos(n\pi x) \qquad (7.14) $$

### 7.6.3 Champs de vecteurs et champs de vitesses (suite)

Reprenons l'exercice 6.3.2. Parmi les champs de vecteurs qui correspondent à un écoulement incompressible quels sont ceux qui correspondent à un écoulement potentiel ? Quelle est l'expression du potentiel correspondant ?
<figure><img src="image_32.png" /><figcaption>FIGURE 7.2 – Additivité des solutions de l'équation de Laplace</figcaption></figure>

## 7.6.4 Réfraction des vitesses

On considère un écoulement dans le plan (Oxy). Dans ce plan deux roches aquifères, de conductivités hydrauliques respectives $K_1$ et $K_2$, présentent une interface le long de la droite (Oy) (figure 7.3). L'angle d'incidence de la vitesse de Darcy U dans le milieu 1 par rapport à l'horizontale (Ox) est $\alpha$. Les charges dans chacun des deux milieux sont respectivement $H_1(x, y)$ et $H_2(x, y)$.

1. Que peut on dire de la charge $H_1$ et $H_2$ sur l'interface, dans les milieux 1 et 2 respectivement ?

2. En déduire une relation entre les gradients de $H_1$ et $H_2$ selon la direction tangentielle y.

3. Relier le gradient tangentiel de $H_1$ au gradient de $H_2$ selon la direction de l'écoulement.

4. Ecrire d'autre part la « loi » de Darcy dans les milieux poreux 1 et 2.

5. Le débit qui traverse l'interface est identique de part et d'autre. Écrire la relation correspondante (entre les débits $Q_{1,x}$ et $Q_{2,x}$ perpendiculaires à l'interface). Que peut-on déduire de l'angle $\beta$ que fait l'écoulement dans le milieu 2 avec l'horizontale ?

<figure><img src="image_33.png" /><figcaption>FIGURE 7.3 – Écoulement dans des milieux de perméabilité différente</figcaption></figure>
# Équipotentielles et lignes de courant

## 8.1 Définitions

Une **équipotentielle** est une courbe de charge constante. Étant donné que les écoulements sont potentiels, la vitesse est en tout point perpendiculaire au courbes équipotentielles. Cette propriété permet, à partir de la connaissance des charges, de connaître le sens de l'écoulement et sa norme.

Une **ligne de courant** est une courbe tangente en tout point au vecteur vitesse à un instant $t$ donné. On verra que tracer ces lignes (souvent à la main) permet en hydrogéologie de se faire une idée rapide (et importante) des sens et de l'ordre de grandeur des écoulements dans bien des situations réelles.

## 8.2 Fonction courant

### 8.2.1 Définition

On ne s'intéresse qu'aux écoulements plans (ici dans le plan horizontal $(xy)$ mais les résultats sont valable dans un autre plan $(xz)$ par exemple). Nous allons définir une fonction appelée fonction courant $\Psi(x, y)$ qui possède un certain nombre de propriétés dont la plus importante est que toute ligne de courant peut s'écrire comme

$$ \Psi(x, y) = C \qquad (8.1) $$

où $C$ est une constante.

Sachant que

$$ \operatorname{div} \mathbf{U} = 0, \qquad (8.2) $$

Il existe une fonction vectorielle $\mathbf{A} = A_x \mathbf{i} + A_y \mathbf{j} + A_z \mathbf{k}$ telle que

$$ \mathbf{U} = \operatorname{rot} \mathbf{A} \qquad (8.3) $$

Or

$$ \operatorname{rot} \mathbf{A} = \left( \frac{\partial A_z}{\partial y} - \frac{\partial A_y}{\partial z} \right) \mathbf{i} + \left( \frac{\partial A_x}{\partial z} - \frac{\partial A_z}{\partial x} \right) \mathbf{j} + \left( \frac{\partial A_y}{\partial x} - \frac{\partial A_x}{\partial y} \right) \mathbf{k} \qquad (8.4) $$

En égalant les composantes de $\mathbf{U}$ avec les termes du rotationnel on obtient

$$ U_x = \frac{\partial A_z}{\partial y} - \frac{\partial A_y}{\partial z} \qquad (8.5) $$

$$ U_y = \frac{\partial A_x}{\partial z} - \frac{\partial A_z}{\partial x} \qquad (8.6) $$

et

$$ 0 = \frac{\partial A_y}{\partial x} - \frac{\partial A_x}{\partial y} \qquad (8.7) $$

Une fonction vectorielle $\mathbf{A}$ telle que

$$ A = \Psi(x, y)\mathbf{k} \qquad (8.8) $$

et

$$ U_x = \frac{\partial \Psi}{\partial y}, \quad U_y = -\frac{\partial \Psi}{\partial x} \qquad (8.9) $$
satisfait à la fois à la définition de la vitesse et aux contraintes imposées par la conservation de la masse. La fonction $\Psi(x, y)$ est appelée **fonction courant**

En comparant avec la définition de la vitesse de Darcy on obtient

$$ -K \frac{\partial \phi}{\partial x} = \frac{\partial \Psi}{\partial y}, \quad K \frac{\partial \phi}{\partial y} = \frac{\partial \Psi}{\partial x} \qquad (8.10) $$

### 8.2.2 Propriétés

1. Une conséquence de ce qui précède est que la fonction courant $\Psi$, comme la charge $\phi$ satisfait à l'équation de Laplace

$$ \Delta \Psi = 0 \qquad (8.11) $$

$\Psi$ et $\phi$ sont des fonctions dites harmoniques. Elles ont comme propriétés que les courbes du plan $\Psi(x, y) = C_1$ et $\phi(x, y) = C_2$ sont orthogonales (ici $C_1$ et $C_2$ sont deux constantes quelconques).

2. Une ligne de courant est une courbe sur laquelle la valeur de $\psi$ reste constante. Une ligne de courant est tangente au vecteur vitesse en chacun de ses points ce qui implique que

$$ \frac{dx}{U_x} = \frac{dy}{U_y}. \qquad (8.12) $$

3. Enfin, le débit s'écoulant entre entre deux lignes de courant est donné par (cf exercice 12)

$$ Q = \int_{\Psi_1}^{\Psi_2} dQ = \Psi_2 - \Psi_1 \qquad (8.13) $$

## 8.3 Exemple

<figure><img src="image_34.png" /><figcaption>FIGURE 8.1 – Courbes équipotentielles $K\phi(x) = 10^{-4}(x + 2y) = C_e \text{ m}^2/\text{s}$ (bleu), lignes de courant $\psi(x, y) = 10^{-4}(2x - y) = C_c \text{ m}^2/\text{s}$ (rouge) et vitesse de l'écoulement. la vitesse est définie par $\mathbf{U} = -10^{-4}\mathbf{i} - 2 \cdot 10^{-4}\mathbf{j}$. Elle a pour norme $\|\mathbf{U}\| = \sqrt{5} \cdot 10^{-4} \text{ m/s}$.</figcaption></figure>

Soit un écoulement plan (xy) dont la charge est donnée par

$$ K\phi = ax + by, \qquad (8.14) $$

ou encore

$$ \phi = (ax + by)/K, \qquad (8.15) $$
$$+ h1(x1, y1) = 100m$$

$$+ h2(x2, y2) = 90m$$

<figure><img src="image_35.png" /><figcaption>FIGURE 8.2 -</figcaption></figure>

où $(a, b)$ sont deux constantes ayant la dimension d'une vitesse. On vérifie bien que cette charge correspond à un écoulement potentiel c'est à dire qu'elle vérifie l'équation de Laplace. Une équipotentielle étant une ligne de charge constante on aura

$$\phi(x) = C_e \qquad (8.16)$$

où $C_e$ est une constante quelconque. Soit ici

$$y = \frac{C_e - ax}{b} \qquad (8.17)$$

Comme la vitesse est perpendiculaire à ces équipotentielles, on peut simplement calculer la vitesse en tout point grâce à

$$U = -K \operatorname{grad} \phi = -ai - bj, \qquad (8.18)$$

dont on déduit sa norme

$$\|U\| = \sqrt{a^2 + b^2}. \qquad (8.19)$$

On peut en outre tracer les lignes de courant qui, parce qu'elles sont tangentes en tout point au vecteur vitesses, correspondent aux courbes perpendiculaires en tout point aux équipotentielles. À l'instar des équipotentielles, les lignes de courant sont des courbes le long desquelles la fonction courant prend une valeur constante $C_c$. Pour trouver leur expression nous allons chercher les fonctions $\psi(x, y)$ respectant (8.10). Elles auront pour équation générale

$$\psi(x, y) = K(bx - ay) = C_c. \qquad (8.20)$$

La figure 8.1 montre différentes équipotentielles, des lignes de courant et une représentation du vecteur vitesse pour une charge définie par $\phi(x, y) = x + 2y$ et $K = 10^{-4} \text{ m/s}$.

# 8.4 Exercices

## 8.4.1 Calcul d'un vecteur vitesse

On dispose de trois valeurs de la charge $h$ mesurées dans trois piézomètres et représentées sur la figure 12.16. On considérera que l'écoulement est horizontal (approximation de Dupuit). Évaluez et représentez la vitesse de l'écoulement dans la nappe pour une perméabilité $K = 10^{-5} \text{ m/s}$?

## 8.4.2 Calcul de vitesse dans une nappe captive

Cinq piézomètres ont été installés dans une nappe captive horizontale (figure 8.3) afin d'en étudier l'écoulement. Sur la figure les points rouges représentent la position des piézomètres, le texte donne la mesure de la cote NGF du niveau d'eau dans les piézomètres (on parle de cote piezométrique).

1. Définissez (judicieusement) trois triangles dont chacun des sommets correspond à un piézomètre. Les trois triangles devront comporter un sommet commun et ne se superposeront pas.

2. Pour chaque triangle estimez la vitesse de l'écoulement au sein du triangle et positionnez là sur la carte. Vous justifierez vos approximations éventuelles et votre démarche.
3. Que constatez-vous ?

4. Tracez les équipotentielles correspondant au cotes piézométriques 21,22 et 23 m en conséquence.

<figure><img src="image_36.png" /><figcaption>FIGURE 8.3 – carte des piézomètres</figcaption></figure>

### 8.4.3 Propriétés de la fonction courant

1. Montrer que la fonction courant d'un écoulement potentiel suit l'équation de Laplace.

$$ \Delta\psi = 0 \qquad (8.21) $$

2. Montrer que le débit passant entre de lignes de courant de valeur $Psi_1$ et $Psi_2$ est donné par

$$ Q = \int_{\Psi_1}^{\Psi_2} dQ = \Psi_2 - \Psi_1 \qquad (8.22) $$

3. Montrer que pour qu'une fonction courant existe il faut que $\nabla \cdot \mathbf{u} = 0$
8.4.4 Étude d’écoulements

1. Reprenons l'exercice 6.3.2. Calculez la fonction courant des champs de vecteurs correspondant à un écoulement potentiel.

2. Soit la fonction $\Psi(x, y) = x^3 - y^3$ montrer qu'elle ne représente pas la fonction courant d'un écoulement potentiel.

3. La fonction $\Phi(x, y) = x^2 + y^2$ est elle solution de l'équation de Laplace ? Si oui quelle est la fonction courant correspondante ?

4. Soit la fonction $\Phi(x, y) = ax + by$.

(a) Montrer que $\Phi$ suit l'équation de Laplace.

(b) Calculer la fonction courant $\Psi(x, y)$ correspondante (on supposera que le milieu a une conductivité hydraulique $K$).

(c) Tracer le réseau des équipotentielles et des lignes de courants dans les cas suivants
— $a = 0, b = 1, K = 1\,\text{m/s},$
— $a = 1, b = 2, K = 0.5\,\text{m/s}$
On prendra un espacement $\Delta\Phi = 10\,\text{m}$ et $\Delta\Psi = 10\,\text{m}^2/\text{s}$

(d) Quel est le débit passant entre deux lignes de courant ?

5. La fonction courant d'un écoulement plan est $\Psi(x, y) = b(x^2 - y^2)$.

(a) Déterminer u en un point.

(b) Quelle est l'expression de la charge piézométrique $\Phi$ (on appellera $K$ la conductivité hydraulique du milieu).

(c) Tracer le réseau des lignes de courant et des équipotentielles en prenant un espacement $\Delta\Phi = 10\,\text{m}$ et $\Delta\Psi = 10\,\text{m}^2/\text{s}$ pour $K = 1\,\text{m/s}$ et $b = 1$.

6. Soit la fonction $\Phi(x, y) = x(x^2 - 3y^2)$

(a) Est elle solution de l'équation de Laplace ?

(b) Si oui quelle est la fonction courant correspondante ?

(c) Tracer le réseau des lignes de courant et des équipotentielles en posant $b = 1$ et en prenant un espacement $\Delta\Phi = \Delta\Psi = 1m$ ($K = 1m/s$)

7. Soit la fonction $\Phi(r) = \ln(r)$ en coordonnées polaires $(r, \theta)$

(a) Est elle solution de l'équation de Laplace ?

(b) Si oui quelle est la fonction courant correspondante ?

(c) Tracer le réseau des lignes de courant et des équipotentielles en posant $b = 1$ et en prenant un espacement $\Delta\Phi = \Delta\Psi = 1m$ ($K = 1m/s$)
# Écoulements autour d'un puits

## 9.1 Introduction

Les puits sont essentiels en hydrogéologie pour plusieurs raisons. On s'en sert pour

1. suivre le niveau piézométrique dans une nappe,

2. pomper l'eau de la nappe,

3. contenir une pollution ou dépolluer une nappe.

Avant de nous concentrer sur les puits à proprement parler nous allons nous intéresser aux solutions de l'équation de Laplace en coordonnées polaires.

## 9.2 Écoulements plan en coordonnées polaires

L'équation de Laplace en coordonnées polaires s'écrit

$$ \frac{\partial^2 \phi}{\partial r^2} + \frac{1}{r} \frac{\partial \phi}{\partial r} + \frac{1}{r^2} \frac{\partial^2 \phi}{\partial \theta^2} = 0 \qquad (9.1) $$

Nous allons chercher des solutions de cette équation en séparant les variables c'est à dire de la forme $\phi(r, \theta) = R(r)\Theta(\theta)$. La séparation (exercice !) donne deux systèmes d'équations différentielles ordinaires

$$ r^2 R'' + rR' - kR = 0, \qquad (9.2) $$

$$ \Theta'' + k\Theta = 0 \qquad (9.3) $$

où $k$ est une constante.

Si $k=0$ alors une solution du système formé par (9.2) et (9.3) sera de la forme

$$ \phi(r, \theta) = A \ln(r) + B \qquad (9.4) $$

D'autre part la symétrie radiale impose que $\Theta(\theta) = \Theta(\theta + 2\pi)$ donc $k = \alpha^2 > 0$ avec $\alpha > 0$ et

$$ \Theta = A_{\alpha} \cos(\alpha\theta) + B_{\alpha} \sin(\alpha\theta) \qquad (9.5) $$

Comme on a forcément $\Theta(\theta) = \Theta(\theta+2\pi)$ cela limite les valeurs possibles de $\alpha$ qui doivent être entières $\alpha = n \in \mathbb{N}$ d'où

$$ \Theta = A_n \cos(n\theta) + B_n \sin(n\theta) \qquad (9.6) $$

Pour R on a

$$ r^2 R'' + rR' - n^2 R = 0 \qquad (9.7) $$

dont une solution générale est de la forme

$$ R(r) = C_n r^n + D_n r^{-n}. \qquad (9.8) $$

Dans ces équations $A, B, A_n, B_n, C_n, D_n$ sont des constantes qui dépendent des conditions au limites et de la valeur de $\alpha$. La solution générale aura alors une forme ressemblant à
$$ \phi(r, \theta) = A \ln(r) + B + \sum_{n=0}^{\infty} (C_n r^n + D_n r^{-n}) (A_n \cos(n\theta) + B_n \sin(n\theta)) \quad (9.9) $$

Intéressons-nous à deux situations simples pour comprendre l'importance des conditions aux limites.

## 9.3 Écoulement radial et uniforme dans un anneau

<figure><img src="image_37.png" /><figcaption>FIGURE 9.1 – Écoulement autour d'un puits (a) dans une nappe captive, (b) dans une nappe libre.</figcaption></figure>

Penchons nous sur deux cas qui ont des applications pratiques : celui d'une nappe captive d'épaisseur *D* (figure 9.1a) d'une part, celui d'une nappe libre (figure 9.1b) d'autre part. Nous pouvons à l'aide de l'équation de Laplace en coordonnées polaire étudier l'écoulement sur le domaine annulaire qui s'étend du puits de rayon *R*<sub>p</sub>, se trouvant au centre jusqu'à la limite matérialisée par le cercle de rayon *R*<sub>a</sub> (figure 9.2 pour la clareté de la figure, le rayon du puits est considérablement exagéré.).

### 9.3.1 Cas d'une nappe captive

Dans cas d'un écoulement confiné le problème aux limites est le suivant.
<figure><img src="image_38.png" /><figcaption>FIGURE 9.2 – Domaine de définition et conditions aux limites de charge d'un écoulement dans un anneau.</figcaption></figure>

### Problème aux limites d'un écoulement annulaire confiné

$$ \Delta\phi = 0 \qquad (9.10) $$

$$ \phi(R_a) = \Phi_a \qquad (9.11) $$

$$ \phi(R_p) = \Phi_p \qquad (9.12) $$

Les conditions aux limites permettent de simplifier considérablement la solution générale (9.9) de ce problème. En effet la charge étant constante sur les deux anneaux le problème est axisymétrique (symétrie de rotation autour de l'axe vertical) et ne dépend donc pas de $\theta$. (9.9) se réduit alors à

$$ \phi(r, \theta) = A \ln(r) + B \qquad (9.13) $$

les valeurs des deux constantes sont données par les conditions aux limites si bien qu'à la fin la solution devient

$$ \phi(r) = \left( \frac{\Phi_a - \Phi_p}{\ln(R_a/R_p)} \right) \ln r + \left( \frac{\Phi_p \ln R_a - \Phi_a \ln R_p}{\ln(R_a/R_p)} \right) \qquad (9.14) $$

### 9.3.2 Cas d'une nappe libre

Dans le cas d'une nappe libre (figure9.1b) l'épaisseur de la nappe varie. Dans la cadre de l'approximation de Dupuit on considérera que l'écoulement est partout horizontal. Dans ce cas la charge correspond à la hauteur d'eau $\phi = h$. Le problème aux limites devient alors

#### Problème aux limites d'un écoulement annulaire libre

$$ \Delta h^2 = 0 \qquad (9.15) $$

$$ h(R_a) = h_a \qquad (9.16) $$

$$ h(R_p) = h_p \qquad (9.17) $$

Dans ce cas on peut se ramener à l'équation de Laplace par un changement de variable très simples $X = h^2$ avec $X(R_a) = h_a^2$ et $X(R_p) = h_p^2$. La solution du problème est alors donnée par
$$h(r) = \sqrt{\left(\frac{h_a^2 - h_p^2}{\ln(R_a/R_p)}\right) \ln r + \left(\frac{h_p^2 \ln R_a - h_a^2 \ln R_p}{\ln(R_a/R_p)}\right)}. \qquad (9.18)$$

## 9.4 Loi de Darcy et débit d'un puits

**Nappe captive** L'application de la loi de Darcy à la solution d'un écoulement dans un anneau donne

$$\mathbf{U} = -K \left( \frac{\Phi_a - \Phi_p}{\ln(R_a/R_p)} \right) \frac{1}{r} \mathbf{u}_r, \qquad (9.19)$$

où $K$ est la conductivité hydraulique. Le signe négatif indique que le vecteur vitesse est, dans notre cas, dirigé vers l'intérieur c'est à dire vers le puits.

le débit s'obtient alors par intégration du produit $Q_e = \mathbf{U} \cdot \mathbf{ds}$ représenté figure 9.1a. connaissant l'élément de surface $ds = r d\theta dz$ on en déduit

$$Q_e = -K \left( \frac{\Phi_a - \Phi_p}{\ln(R_a/R_p)} \right) \int_0^D \int_0^{2\pi} dz d\theta, \qquad (9.20)$$

soit

$$Q_e = -2\pi KD \left( \frac{\Phi_a - \Phi_p}{\ln(R_a/R_p)} \right) \qquad (9.21)$$

Rappelons que par convention un flux entrant une surface fermée est négatif. Le flux pompé en surface sera $Q_s = -Q_e$ et sera lui positif soit

$$Q_s = 2\pi KD \left( \frac{\Phi_a - \Phi_p}{\ln(R_a/R_p)} \right). \qquad (9.22)$$

Je laisse le calcul en application mais à partir de cette dernière équation il est possible de retrouver l'expression classique —que l'on trouve dans tous les livres d'hydrogéologie— de la charge autour d'un puits en nappe captive :

$$\phi(r) = \frac{Q_s}{2\pi T} \ln r + C. \qquad (9.23)$$

où $T = KD$ et $C$ est une constante correspondant au second membre de l'équation (9.14).

**Nappe libre** Les calculs à partir de (9.18) sont un peu plus long mais ne présentent aucune difficulté. On obtient successivement

$$\mathbf{U} = -\frac{K}{2r} \left( \frac{h_a^2 - h_p^2}{\ln(R_a/R_p)} \right) \frac{1}{h(r)} \mathbf{e}_r, \qquad (9.24)$$

puis

$$Q_s = \pi K \left( \frac{h_a^2 - h_p^2}{\ln(R_a/R_p)} \right), \qquad (9.25)$$

et enfin l'équation «classique» de la charge autour d'un puits en nappe libre.

$$h^2(r) = \frac{Q_s}{\pi K} \ln r + C. \qquad (9.26)$$

**Notes** Contrairement à la solution logarithmique pour une nappe captive, la solution pour la nappe libre n'est pas solution de l'équation de Laplace à proprement parler (c'est le carré de la charge qui est solution de l'équation de Laplace). Ceci est lié à l'approximation de Dupuit en nappe libre qui rend l'équation différentielle (9.18) non linéaire en $h$ (elle est linéaire en $h^2$) alors que l'équation de Laplace est linéaire en $h$. Nous retiendrons donc que cette solution logarithmique pour une nappe libre est approximative même si dans la majeure partie des situation elle est largement suffisante.
<figure><img src="image_39.png" /><figcaption>FIGURE 9.3 – Pompage dans un île.</figcaption></figure>

Les équations (9.26) et (9.23) sont trompeuses car elles laissent croire que la charge croit infiniment (quand $r \to \infty$, $h \to \infty$). De même quand $r \to 0$ alors $h \to -\infty$. Or c'est impossible physiquement. Un puits ayant un certain rayon non nul, le problème posé par la limite à 0 ne se pose pas. Quand $h$ augmente par contre on est confronté au problème de l'augmentation constante de la charge. Comme nous l'avons vu, cette solution logarithmique correspond à l'écoulement dans un anneau donc dans un domaine limité aux bords duquel une charge fixe est imposée (on pourrait tout aussi bien imposer un flux). c'est le cas par exemple d'une île entourée d'eau où le niveau du lac, ou de la mer, fixe la charge sur le pourtour de l'île (figure 9.3).

Dans la pratique vous entendrez parler d'un rayon d'influence ou rayon d'action d'un puits. Pour résoudre ce problème de la croissance du logarithme les hydrogéologues font l'hypothèse que l'influence du puits devient nulle au delà d'un certain rayon $R$ le rayon d'influence ou d'action du puits pour lequel la charge $\phi = \phi_a$ (resp. $h = h_a$ en nappe libre) est fixée. Cela revient à se ramener à un problème aux limites d'écoulement dans un anneau dont le rayon extérieur correspond à ce fameux rayon d'influence.

## 9.5 Rabattement et importance des puits

### 9.5.1 Principe

La figure 9.4 représente une nappe captive à deux instants. Dans un état initial l'eau s'écoule dans la nappe selon la dynamique étudiée précédemment. la charge varie linéairement (figure 9.4a). On installe dans cette nappe un puits qui pompe un débit $Q$ dans cette nappe. La figure 9.4b représente l'état des écoulements dans la nappe quand ceux-ci atteignent de nouveau un régime stationnaire. On constate deux changements majeurs

1. la surface piézométrique n'est plus linéaire. Elle présente ce que l'on nomme un rabattement qui s'accentue à l'approche du puit.

2. le sens des écoulements s'inverse localement et une limite physique apparaît à gauche du puits qui sépare les eaux qui s'écoulent vers la gauche de celles qui s'écoulent en direction du puits.

Cette seconde constatation montre qu'avec un puits on peut isoler une partie d'une nappe et influer sur le sens des écoulements. On peut donc se servir d'un (ou de plusieurs puits) pour protéger une nappe en cas de contamination.

En effet si une pollution advient, dans le cas d'une nappe sans puits, la pollution pourra se propager avec l'écoulement dans la nappe. Par contre si un puits est mis en place de façon à ce que la zone polluée soit située à droite de la ligne de partage des eaux de la figure 9.4b alors le changement de direction de l'écoulement provoqué par la mise en place du forage permettra d'empêcher la propagation de la pollution.

### 9.5.2 Cas de l'île

l'additivité des solutions pour $h^2$ dans le cas d'un écoulement pour lequel l'hypothèse de Dupuit s'applique nous permet d'obtenir une solution pour la hauteur d'eau dans une nappe libre recevant une précipitation $P$ et dans laquelle on pompe un débit $Q$.

$$h^2 = \frac{P}{2K} (R^2 - r^2) + \frac{Q_s}{\pi K} \ln \left( \frac{r}{R} \right) \qquad (9.27)$$

La figure 9.5 montre le résultat (pas à l'échelle). On voit comment l'équilibre entre précipitations et pompages génère un cône ...
<figure><img src="image_40.png" /></figure>

FIGURE 9.4 – (a) Écoulement dans une nappe captive en régime stationnaire. (b) Rabattement sous l’effet d’un pompage. La courbe en pointillé indique le profil de charge avant l’installation du puits. La ligne en trait plein indique le profil de charge après mise en place du pompage.

<figure><img src="image_41.png" /></figure>

FIGURE 9.5 – Pompage de débit Q dans une île de rayon R, de conductivité hydraulique K recevant une précipitation P
## 9.6 Exercices

### 9.6.1 Écoulements autour d'un puits

Dans chacun des cas suivants proposez une solution approximée au problème posée puis comparez cette
solution au résultat obtenu avec un potentiel décrit par une fonction logarithmique.

1. La pente de la surface piézométrique d'un aquifère confiné d'épaisseur constante $D$, mesurée à $r = 5$ m d'un puits pompant un débit $Q = 200$ m³/h, est $dh/dr = 0.2$ ($h$ et $r$ correspondent respectivement à la hauteur piézométrique et à la distance au puits). Calculer la transmissivité $T = KD$ de l'aquifère.

2. Déterminer le débit d'un puits en nappe captive compte tenu des informations suivantes : Différence de hauteurs piézométriques de $\Delta h = 2.5$ m entre deux piézomètres situés respectivement à $r_1 = 10$ m et $r_2 = 30$ m du centre du puits ; épaisseur de la nappe : $D = 30$ m ; conductivité hydraulique : $K = 0.0001$ m/s. On suppose que les mesures sont réalisées en régime permanent.

3. Déterminer le coefficient de perméabilité dans une nappe captive si les lectures à partir du niveau statique dans deux piézomètres situés respectivement à $r_1 = 20$ et $r_2 = 150$ m du puits sont, dans l'ordre 3.3 et 0.3 m. L'épaisseur de la nappe est de $D = 30$ m. et le débit est de $Q = 0.2$ m³/s.

### 9.6.2 Alimentation en eau potable

Une municipalité compte 15 000 habitants, sa population croit et on s'attend à atteindre une population de 22 000 habitants dans les années qui arrivent. Aujourd'hui la consommation globale par habitant est stable à 375 litres par jour. La municipalité souhaite savoir si l'installation d'alimentation en eau potable actuelle sera toujours suffisante quand la population aura atteint la taille estimée.

La ressource en eau est constituée par une nappe libre, captée au moyen d'un puits de 50 cm de rayon. La charge piézométrique disponible est de 30 mètres et le coefficient de perméabilité est de 0.0003 m/s.

1. Calculer la consommation d'eau potable annuelle totale actuelle de la commune.

2. Calculer le débit actuel du puits

3. Sachant que le rayon d'influence<fnref n="1" /> est de 500 m pour le débit actuel quel est le rabattement,

4. Évaluer la vitesse de filtration à proximité du puits,

5. La nappe peut-elle alimenter 22000 personnes en eau ?

<footnote n="1">1. Le rayon d'influence est une notion courante en hydrogéologie. Il représente le rayon au delà duquel un puits n'a plus d'influence sur la nappe c'est à dire le rayon ou la charge de la nappe n'est pas altérée par le pompage. Sauf en de rares cas cette notion n'a pas de sens même si elle est utile en pratique.</footnote>
# 10 Écoulement dans une nappe captive

<figure><img src="image_42.png" /><figcaption>FIGURE 10.1 – Écoulement dans une nappe captive : 1) piézomètres, 2) couches imperméables délimitant l'aquifère captif, 3) zone d'étude de l'écoulement, 4) Niveau piézométrique hypothétique.</figcaption></figure>

## 10.1 Description du problème et conditions aux limites

Considérons une coupe dans une nappe captive (figure 10.1), d'épaisseur $D$, c'est à dire un aquifère saturé, confiné entre deux couches imperméables (aquicludes). Nous allons chercher à décrire la répartition des charges piézométriques dans la nappe ce qui nous permettra d'en déduire les écoulements. Nous faisons l'hypothèse que la nappe est infinie dans la direction perpendiculaire au plan de coupe et que l'écoulement est bidimensionnel dans le plan. Nous allons étudier l'écoulement dans la portion rectangulaire délimitée par deux côtés verticaux correspondant à deux piézomètres et deux côtés horizontaux correspondant aux limites imperméables (figure 10.1). Cette situation est intéressante car elle permet d'aborder la recherche de solutions analytiques de l'équation de Laplace dans un plan. Afin de chercher des solutions la seule connaissance de l'équation de Laplace ne suffit pas il faut connaitre des conditions aux limites du système que l'on étudie.

Nous en avons deux type différents dans notre cas.

1. Sur les limites latérales (côtés verticaux) la charge $\phi$ est une constante (une mesure donnée par le piézomètre)

$$ \phi(x_1, z) = h_1, \phi(x_2, z) = h_2 \qquad (10.1) $$

Ces conditions sont appelées des conditions de Dirichlet, c'est à dire des conditions sur la charge

2. Les limites supérieures (côtés horizontaux) et inférieures de la nappe sont des limites à flux nul (rien ne passe). Elles correspondent à des lignes de courant dans le cas de la figure 10.1. La composante perpendiculaire de la vitesse est donc nulle ce que l'on peut écrire

$$ U_z = \partial_z \phi = 0 \qquad (10.2) $$

Ces conditions qui porte sur des dérivées de la fonction potentielle, c'est à dire sur des flux ou les vitesses, sont appelées des conditions de Neuman
L'équation de Laplace munie de ses conditions aux limites représente ce que l'on appelle un problème aux limites (PL) ou Boundary Value Problem (BVP) en anglais

## Problème aux limites de la nappe captive

$$ \Delta\phi = 0, \qquad (10.3) $$

$$ \phi(x_1, z) = h_1, \qquad (10.4) $$

$$ \phi(x_2, z) = h_2, \qquad (10.5) $$

$$ \partial_z\phi(z_1) = 0, \qquad (10.6) $$

$$ \partial_z\phi(z_2) = 0. \qquad (10.7) $$

## 10.2 Résolution

Il n'y a pas de méthode générale de résolution d'une équation aux dérivées partielles mais des méthodes qui marchent...ou pas. Ici nous allons utiliser la technique la «plus simple» dite de séparation des variables. On écrira donc

$$ \phi(x, z) = X(x)Z(z), \qquad (10.8) $$

la charge $\phi$ peut s'écrire comme le produit de deux fonctions $X(x)$ et $Y(y)$ qui ne dépendent que d'une variable. L'équation de Laplace se réécrit alors

$$ X''Z + Z''X = 0 \qquad (10.9) $$

On cherche des solutions non nulles X et Y donc on peut réarranger et écrire

$$ \frac{X''}{X} = -\frac{Z''}{Z} = k \qquad (10.10) $$

où $k$ est une constante. Trois types de solutions sont possibles en fonction de la valeur de $k$.

### $k < 0$ on peut écrire $k = -\alpha^2$ où $\alpha > 0$

$$ X = A_x \cos(\alpha x) + B_x \sin(\alpha x), \qquad (10.11) $$

et

$$ Z = A_z \cosh(\alpha z) + B_z \sinh(\alpha z). \qquad (10.12) $$

La vitesse de Darcy étant nulle aux limites supérieures et inférieures de la nappe on a

$$ \partial_z\phi(0) = Z'(0) = \partial_z\phi(D) = Z'(D) = 0 \qquad (10.13) $$

soit

$$ \alpha (A_z \sinh(\alpha z) + B_z \cosh(\alpha z)) = 0, \quad z \in \{0, D\} \qquad (10.14) $$

On peut vérifier que la double condition implique $A_z = B_z = 0$ soit $Z = 0$ ce qui conduit à un champs nulle en tout point donc exit $k < 0$

### $k > 0$ on peut écrire $k = \alpha^2$ où $\alpha > 0$

$$ X = A_x \cosh(\alpha x) + B_x \sinh(\alpha x), \qquad (10.15) $$

et

$$ Z = A_z \cos(\alpha z) + B_z \sin(\alpha z). \qquad (10.16) $$

On a toujours nos conditions de flux nul en haut et en bas de la nappe soit

$$ \alpha (A_z \sin(\alpha z) - B_z \cos(\alpha z)) = 0, \quad z \in \{0, D\} \qquad (10.17) $$

On peut vérifier que la double condition, couplée cette fois à la constance de la charge aux limites en $x_1$ et $x_2$, implique de nouveau $A_z = B_z = 0$ soit $Z = 0$ ce qui conduit à un champs nulle en tout point. Exit donc $k > 0$.
k = 0

$$
X = a_x x + b_x \tag{10.18}
$$

$$
Z = a_z z + b_z \tag{10.19}
$$

L'application de la condition sur les limites supérieures implique $a_z = 0$ soit $Z = b_z$.
En injectant alors dans (10.18) on obtient

$$
\phi(x) = b_z (a_x x + b_x) = ax + b, \tag{10.20}
$$

avec $a = b_z a_x$, $b = b_z b_x$. Les conditions limites permettent alors de retrouver la solution du cas 1D développé au chapitre précédent puis de calculer la fonction courant.

Charge et fonction courant dans une nappe captive

$$
\phi = \left( \frac{h_2 - h_1}{x_2 - x_1} \right) (x - x_1) + h1. \tag{10.21}
$$

$$
\psi = -K \left( \frac{h_2 - h_1}{x_2 - x_1} \right) z + C, \quad (10.22)
$$

où K est la conductivité hydraulique du milieu et C une constante d'intégration arbitraire. La figure 10.2 montre le résultat. Yes!

<figure><img src="image_43.png" /><figcaption>FIGURE 10.2 – Écoulement dans une nappe captive : 1) piézomètres, 2) couches imperméables délimitant l’aquifère captif, 3) carte des équipotentialles et des lignes de courant, 4) Niveau piézométrique.</figcaption></figure>

10.3 Débit de la nappe

Nous avons vu au chapitre précédent la relation existant entre le débit et les lignes de courant (8.13). En effectuant l’intégration sur l’épaisseur de la nappe on obtient

$$
q = \int_{\psi(z_1)}^{\psi(z_2)} d\psi = -KD \left( \frac{h_2 - h_1}{x_2 - x_1} \right) \quad (10.23)
$$

Rappelons que le débit **par unité de largeur** *q* est en (en m²/s)

On peut bien sûr retrouver ce débit en passant par le calul de la vitesse,

$$
\mathbf{U} = -K \frac{d\phi}{dx} \mathbf{i}. \tag{10.24}
$$

En intégrant sur l'épaisseur de la nappe on retrouve (10.23) que l'on peut réécrire
$$q = -T \left( \frac{h_2 - h_1}{x_2 - x_1} \right). \qquad (10.25)$$

La constante $T = KD$ est appelée la **transmissivité** de la nappe.
# Écoulement sous une masse d'eau

<figure><img src="image_44.png" /><figcaption>FIGURE 11.1 – Écoulement sous une masse d'eau : une mare, un lac, une retenue ou une rivière</figcaption></figure>

Ce chapitre va être l'occasion pour nous de travailler plus avant sur les solutions de l'équation de Laplace et de réfléchir à l'impact des conditions au limites sur la forme d'un écoulement et d'un champs potentiel. Evidemment les cas que nous étudierons sont simplifiés à l'extrême mais, on le verra, viser l'épure peut s'avérer riche en enseinements.

## 11.1 Problème(s) à résoudre

On va considérer un problème bidimensionnel (le lac ou la rivière sont infinis dans le sens perpendiculaire à notre feuille). La figure 11.1 permet de réfléchir au problème. Le lac constitue une limite à charge constante donnée par la surface libre du lac. Cependant comme on le voit nous avons besoin de trois autres conditions pour clore le problème. Plusieurs solutions se présentent.

1. Les conditions qui rende le problème le plus simple mathématiquement correspondent à trois côtés de charge constante et nulle.

2. Une condition intéressante et plus réaliste consiste à changer la condition de fond pour la faire correspondre à un milieu imperméable. Dans ce cas on a une condition sur le flux et non plus la charg.

3. Enfin on peut considérer des charges non nulles sur les côtés qui traduisent le fait que notre lac peut-être en partie alimenté par sa nappe sur une partie et l'alimenter sur une autre

Dans la «vrai vie» ces conditions dépendent pour une part essentielle des mesures dont on dispose. Nous détaillerons le premier cas puis nous irons plus vite sur les autres qui peuvent être taités comme autant d'exercices.
## 11.2 Lac alimentant sa nappe

### 11.2.1 Problème à résoudre

Cette situation correspond au (1) le problème à résoudre est alors

Problème aux limites d'un lac alimentant sa nappe (modèle 1)

$$ \Delta\phi = 0, \qquad (11.1) $$

$$ \phi(x, H) = \phi_1, \qquad (11.2) $$

$$ \phi(x, 0) = 0, \qquad (11.3) $$

$$ \phi(0, z) = 0, \qquad (11.4) $$

$$ \phi(L, z) = 0 \qquad (11.5) $$

Le lac est de longueur *L* et la charge s'annule à une profondeur *H*.

### 11.2.2 Solution analytique

#### Séparation des variables

On cherche une solution $\phi(x, z) = X(x)Z(z)$ comme dans le cas de la nappe captive. Là encore l'équation de Laplace se transforme en deux équations différentielles ordinaire

$$ X''(x) = kX(x) \qquad (11.6) $$

$$ Z''(z) = -kZ(z) \qquad (11.7) $$

dont les solutions dépendent de *k*. On recommence donc par étudier chacun des trois cas possibles.

* $k < 0$ On peut alors écrire $k = -\alpha^2$ avec $\alpha > 0$ et

$$ X = A_x \cosh(\alpha x) + B_x \sinh(\alpha x) \qquad (11.8) $$

les conditions aux limites en 0 et L imposent $A_x = B_x = 0$ donc exit ce cas.

* $k = 0$ Dans ce cas $X = ax + b$ et idem $a = b = 0$ donc exit ce cas.

* $k > 0$ On peut alors écrire $k = \alpha^2$ avec $\alpha > 0$ et

$$ X(x) = A_x \cos(\alpha x) + B_x \sin(\alpha x) \qquad (11.9) $$

la première condition donne $A_x = 0$ par contre la seconde donne un résultat intéressant.

$$ X(L) = B_x \sin(\alpha L) = 0 \qquad (11.10) $$

Il faut donc que la constante soit nulle (aucun intérêt sinon on a plus de solution) ou que le sinus s'annule soit

$$ \alpha = \frac{n\pi}{L}, n \in \mathbb{N}_{*} \qquad (11.11) $$

la solution sera donc une somme de sinus puisque alpha peut prendre toutes les valeurs décrites par 11.11. D'la même façon «la» constante $B_x$ est en fait une somme de constantes que nous nommerons $B_n$ une par valeur de *n*.

Connaissant cette contrainte nous pouvons nous attaquer à $Z(z)$. On a

$$ Z_n(z) = C_n \cosh \left( \frac{n\pi}{L} z \right) + D_n \sinh \left( \frac{n\pi}{L} z \right) \qquad (11.12) $$

Comme $Z(x, 0) = 0$ on a $C_n = 0$. Reste donc le sinus hyperbolique. La solution du problème devient donc enfin

$$ \phi(x, z) = \sum_{n} B_{n} D_{n} \sinh\left(\frac{n\pi}{L} z\right) \sin\left(\frac{n\pi}{L} x\right) \qquad (11.13) $$
$B_n$ et $D_n$ étant deux constantes dépendant de $n$ on pose $E_n = B_n D_n$. Pour finir de résoudre le problème nous allons utiliser la dernière condition pour trouver l'expression des constantes $E_n$. Sur le fond du lac $\phi(0, H) = \phi_1$ soit

$$ \sum_n E_n \sinh \left( \frac{n\pi H}{L} \right) \sin \left( \frac{n\pi}{L} x \right) = \phi_1 \quad (11.14) $$

Pour résoudre ce problème on utilise la propritété d'orthogonalité des sinus (cf exercice) c'est à dire que

$$ \int_{0}^{\pi} \sin(nx) \sin(mx) dx = \begin{cases} 0 & si \ n \neq m \\ \pi/2 & si \ n = m \end{cases} \quad (11.15) $$

On trouve alors que

$$ E_n \sinh \left( \frac{n\pi H}{L} \right) \int_{0}^{L} \sin^2 \left( \frac{n\pi}{L} x \right) dx = \int_{0}^{L} \phi_1 \sin \left( \frac{n\pi}{L} x \right) dx = \frac{-2L\phi_1}{n\pi} \quad (11.16) $$

or l'intégrale de gauche est nulle pour les valeurs paires de $n$ et vaut $-L/2$ pour les valeurs impairs. la solution du problème s'écrira donc

### Solution du modèle 1

$$ \phi(x, z) = \sum_{n=1}^{\infty} E_n \sinh \left( \frac{(2n-1)\pi}{L} z \right) \sin \left( \frac{(2n-1)\pi}{L} x \right), \quad (11.17) $$

$$ E_n = \frac{-4\phi_1}{(2n-1)\pi \sinh\left(\frac{(2n-1)\pi H}{L}\right)}. \quad (11.18) $$

<figure><img src="image_45.png" /><figcaption>FIGURE 11.2 – Écoulement sous un lac (11.17). Équipotentielles en bleu, lignes de courant en rouge.</figcaption></figure>

### Interprétation et commentaires

**Écoulement.** Evidemment la première question qui vient est : à quoi ressemble le résultat ! La figure 11.2 montre le résultat. Les équipotentielles (lignes de charge constante) sont représentées en bleu, les lignes de courant en rouge. On voit bien que l'eau s'infiltre depuis le lac vers la nappe avec une séparation centrale entre une partie du lac qui est draînée vers la gauche et une partie drainée vers la droite.

**Nombre de termes dans le calcul.** Pour représenter les écoulements nous n'avons pas réalisé une somme infinie comme l'indique la solution (11.17). Nous n'avons sommé qu'un certain nombre de termes. Pourquoi ? En premier lieu parce que nous ne pouvons pas ! Mais si on pousse un peu plus loin le raisonnement et qu'on observe l'expression de la constante (11.18) on voit qu'elle diminue avec $n$ pour tendre vers 0. A posteriori cela semble évident mais c'est une condition nécessaire pour que la somme ne diverge pas. La question est de combien de termes avons nous besoin.

À cette fin nous calculons le rapport entre la chage du lac $\phi_1$ et sa valeur calculée en fonction du nombre de terms inclus dans le calcul. La figure 11.4 montre le résultat. Même si la valeur des coefficients diminue rapidement il faut tout de même en inclure une centaine afin de voir la limite ressembler à ce qu'elle devrait être : un plateau.
### 11.2.3 Solution Numérique

Repenons notre fonction potentielle $\phi$ et effectuons des développements limités à l'ordre 2 (on peut en guise d'exercice effectuer un développement limité à l'ordre 3). On a

$$ \phi(x+dx, z) = \phi(x, z) + \frac{\partial\phi(x, z)}{\partial x}dx + \frac{1}{2}\frac{\partial^2\phi(x, z)}{\partial x^2}(dx)^2, \qquad (11.19) $$

et

$$ \phi(x - dx, z) = \phi(x, z) - \frac{\partial\phi(x, z)}{\partial x}dx + \frac{1}{2}\frac{\partial^2\phi(x, z)}{\partial x^2}(dx)^2. \qquad (11.20) $$

En sommant (11.19) et (11.20) et en réarrangeant on obtient

$$ \frac{\partial^2\phi(x, z)}{\partial x^2} \sim \frac{\phi(x + dx, z) + \phi(x - dx, z) - 2\phi(x, z)}{(dx)^2}. \qquad (11.21) $$

Le même calcul selon $z$ donne

$$ \frac{\partial^2\phi(x, z)}{\partial z^2} \sim \frac{\phi(x, z + dz) + \phi(x, z - dz) - 2\phi(x, z)}{(dz)^2}, \qquad (11.22) $$

soit enfin en sommant les deux dérivées et en les égalant à 0

$$ \phi(x, y) = \frac{(\phi(x + dx, z) + \phi(x - dx, z))dz^2 + (\phi(x, z + dz) + \phi(x, z - dz))dx^2}{(2dx^2 + 2dz^2)}. \qquad (11.23) $$

<figure><img src="image_46.png" /><figcaption>FIGURE 11.3 – Grille numérique (python) pour le calcul de la solution numérique de l'équation de Laplace. Cercles noirs : points intérieurs, carrés bleus : limites du domaine</figcaption></figure>

L'approximation discrète de (11.23) sur une grille $(i, j)$ de maille carrée et de distance entre les noeuds $dx = dy = \delta$ (figure 11.3) aura pour expression

$$ \phi(i, j) = \frac{\phi(i+1, j) + \phi(i-1, j) + \phi(i, j+1) + \phi(i, j-1)}{4} \qquad (11.24) $$

On retrouve bien ce que nous avions dit plus haut à savoir que la valeur du potentiel en un point correspond à la moyenne des potentiels l'entourant.

Pour résoudre numériquement notre problème nous devrons alors
1. définir une grille de noeuds (i, j), de taille ($N_x, N_y$), de distance internoeuds $\delta$ (figure 11.3),

2. définir les bords de la grille et les conditions au limite sur ces bords,

3. remplir la grille d'un ensemble de valeurs initiales quelconques (ici 0),

4. itérer jusqu'à ce que les calculs convergent vers une solution stable.

le script suivant décrit chacune de ces étapes.

```python
def solnum_PL1():
   
   
   
    # Jacobi iterative method to solve laplace equation with
    # Dirichlet type BCs on a rectangular grid
    #
   
    # Grid size
    Nx=50
    Ny=25
    #
    # Coordinates
    xi = np.linspace(0, Nx, Nx)
    yi = np.linspace(0, Ny, Ny)
    x, y = np.meshgrid(xi, yi)
    delta = 1  # dx=dy=1
    #
    # Initialize potential function
    phi = 0 * x
    #
    # Boundaries
    # Initialize Boundary matrix B
    Boundary_points = np.zeros((Ny, Nx))
    #
    # Declare boundary points
    # Boundary_points =1 => Boundary
    # Boundary_points =0 => inner point where the laplacian is to be calculated
    Boundary_points[: ,0] = 1
    Boundary_points[:,Nx-1] = 1
    Boundary_points[0: ,] = 1
    Boundary_points[Ny-1: ,] = 1
    #
    # Dirichlet condition on the limits (change accordingly)
    phi[Ny-1: ,] = 20
    phi[0: ,] = 0
    phi[: ,0] = 0
    phi[: ,Nx-1] = 0
    #
    #convergence criteria
    epsilon=10
    conv Crit = 0.0005
    #
   
    # Loop to solve Delta Phi = 0
    #
    while epsilon > conv Crit:
        phi_copy = phi.copy()
        for i in range(Ny):
            for j in range(Nx):
                if Boundary_points[i,j] == 0:
                    phi[i,j] = (phi_copy[i-1,j]+phi_copy[i+1,j]+phi_copy[i,j-1]+
                                phi_copy[i,j+1])/4
                    epsilon = np .sqrt(np .matrix((phi-phi_copy)**2) .sum() /(Nx*Ny))
                    print(epsilon)
    #
```
```python
# Calculate velocity vectors
```

```python
U=-np.diff(phi, axis=1)/np.diff(x, axis=1)
```

```python
V=-np.diff(phi, axis=0)/np.diff(y, axis=0)
```

```python
# Plot equipotentials and streamlines
```

```python
fig = plt.figure(figsize=(14/2.54, 7/2.54))
```

```python
ax = fig.add_subplot(111)
```

```python
CS = ax.contour(x, y, phi, levels=range(0, 20, 2), colors='blue')
```

```python
ax.clabel(CS, inline=1, fontsize=10)
```

```python
ax.streamplot(x[1:,1:], y[1:,1:], U[1:,:], V[:1:], color='red', density=0.5)
```

```python
ax.axis('square')
```

```python
plt.show()
```

Vous pouvez jouer avec les valeurs des conditions aux limites pour voir à quoi ressemblerait la solution si d'autres côtés ont un potentiel non nul ce que nous verrons pas la suite analytiquement (mais c'est nettement plus long et exigeant !). Les exercices du chapitre ?? nous permettront d'aborder les problèmes de vitesse et de critère convergence des solutions numériques.

### 11.2.4 Limites

Ce que nous venons de faire ne correspond pas à un «vrai» lac... vous en doutez et pour plusieurs raisons :

1. la profondeur d'un lac n'est pas forcément constante,

2. l'écoulement ne se réduit pas à un réangle sous le lac,

3. les conditions aux limites de charge nulle sont peu vraissemblables.

4. enfin il existe une discontinuité dans les conditions aux limites.

Le premier point peut être résolu en prenant une distribution de hauteur, le second aussi en étendant la zone d'étude au delà des limites du lac. Les conditions aux limites de charge nulle on l'a vu sont des conditions dont l'objectif est double : (1) simplifier le problème au maximum, (2) décomposer un problème plus complexe en sous-problèmes plus simples et utiliser la propriété d'additivité des solutions. Enfin nous ne nous attarderons pas trop le dernier point mais nous avons une discontinuité dans les conditions aux limites. Si le raccord entre les trois bords inférieurs, quand $\phi = 0$, est simple, au niveau du lac par contre on passe de 0 à $H$ de façon brutale. Les points $(0, H)$ et $(L, H)$ sont de facto exclus du domaine de solution comme le montre la figure 11.3 présentant la solution numérique.

<figure><img src="image_47.png" /><figcaption>FIGURE 11.4 – Écoulement sous un lac : rapport entre la charge $\phi_1$ et la charge calculée (11.18) en fonction du nombre de termes inclus dans le calcul.</figcaption></figure>
## 11.3 Ajout d'un niveau imperméable

### 11.3.1 Problème à résoudre

Les nappes ne sont pas infinies. les aquifères sont limités par des niveaux moins perméables voire imperméables. Dans le cas de la nappe captive étudiée précédemment l'aquifère était limité par deux niveaux imperméables au dessus (toit) et en dessous (mur) de la nappe. Nous avions vu que les conditions sur ces limites étaient des conditions de flux nul. Nous allons donc remplacer la condition de charge nulle en $z = 0$ par une condition de flux nul. Le problème devient alors

Problème au limites d'un lac alimentant une nappe de profondeur $H$ (modèle 2)

$$ \Delta\phi = 0, \qquad (11.25) $$

$$ \phi(x, H) = \phi_1, \qquad (11.26) $$

$$ \frac{\partial \phi}{\partial z}(x, 0) = 0, \qquad (11.27) $$

$$ \phi(0, z) = 0, \qquad (11.28) $$

$$ \phi(L, z) = 0. \qquad (11.29) $$

Le lac est toujours de longueur $L$ mais cette fois la nappe est limitée par un aquiclude en $z = 0$ c'est à dire à une profondeur $H$ sous le lac.

### 11.3.2 Solution analytique

La solution $X(x)$ n'est pas modifiée par notre changement de conditions aux limites. Seule la solution en $Z(z)$ en effet en $y=0$ on a

$$ \frac{\partial Z_n}{\partial z}(x, 0) = B_n = 0. \qquad (11.30) $$

La solution générale devient alors :

$$ \phi(x, z) = \sum_{n} C_{n} \cosh\left(\frac{n\pi}{L}z\right) \sin\left(\frac{n\pi}{L}x\right). \qquad (11.31) $$

Nous avons échangé un sinus hyperbolique pour un cosinus hyperbolique. En effectuant le calcul pour la constante (exercice !) on arrive enfin à

Solution du modèle 2

$$ \phi(x, z) = \sum_{n=1}^{\infty} C_{2n-1} \cosh\left(\frac{(2n-1)\pi}{L}z\right) \sin\left(\frac{(2n-1)\pi}{L}x\right), \qquad (11.32) $$

$$ C_n = \frac{-4\phi_1}{(2n-1)\pi \cosh\left(\frac{(2n-1)\pi H}{L}\right)}. \qquad (11.33) $$

### 11.3.3 Solution Numérique

Le script suivant montre les changements induits par la prise en compte d'une condition de flux en bas de notre rectangle

1 \# Dirichlet condition on the upper limits

2 \phi[Ny-1,:]=20

3 \phi[:,0]=0

4 \phi[:,Nx-1]=0

5 \# Neumann condition flux=0 on the lower limit

6 dphi = npduos(Nx)

7

8

9 epsilon=10
<figure><img src="image_48.png" /><figcaption>FIGURE 11.5 – Écoulement sous un lac dans une nappe limitée par un niveau imperméable (11.32). Équipotentielles en bleu, lignes de courant en rouge.</figcaption></figure>

```R
conv Crit = 0.0005

while epsilon > conv Crit:
    phi_copy = phi.copy()
    for i in range(Ny):
        for j in range(Nx):
            if Boundary_points[i,j] == 0:
                phi[i,j] = (phi_copy[i-1,j]+phi_copy[i+1,j]+phi_copy[i,j-1]+phi_copy[i,j+1])/4
            # treatment of the Neumann condition on the lower boundary.
            # Beware of the (i,j) inversion with regard to (x,y)
            # i -> y
            # j -> x
            if i==0 and 0 < j < Nx-1:
                phi[i,j] = phi[i+1,j] + delta * dphi[j]
            epsilon = np.sqrt(np.matrix((phi-phi_copy)**2).sum()/(Nx*Ny))
    print(epsilon)
```

### 11.3.4 Analyse

L'ajout d'une niveau imperméable change significativement la forme de l'écoulement. Alors que l'eau s'infiltre dans toutes les directions dans le cas précédent, ici le mur de la nappe induit une séparation plus marquée entre la partie gauche et droite du lac. les lignes de courants sont nettement plus incurvées avec la restriction du domaine d'écoulement.

## 11.4 Ajout de charges non nulles sur les bords du domaine

### 11.4.1 Problème à résoudre

Le lac est dans un champs avec un écoulement préexistant matérialisé par des charges non nulles sur les cotés. Le problème général est le suivant:

Problème au limites d'un lac alimentant une nappe de profondeur H (modèle 3)

$$ \Delta\phi = 0, \qquad (11.34) $$

$$ \phi(x, H) = \Phi_1, \qquad (11.35) $$

$$ \frac{\partial \phi}{\partial z}(x, 0) = 0, \qquad (11.36) $$

$$ \phi(0, z) = \Phi_0, \qquad (11.37) $$

$$ \phi(L, z) = \Phi_2. \qquad (11.38) $$
<figure><img src="image_49.png" /><figcaption>FIGURE 11.6 – Décomposition du problème aux limites des échanges entre un lac et une nappe dans laquelle se produit un écoulement.</figcaption></figure>

Nous allons faire usage de l'additivité pour résoudre le problème en le décomposant en trois sous problèmes (décrits dans la figure 11.6) que nous sommerons ensuite.

### 11.4.2 Solution analytique

Le problème (1) a déjà été résolu (modèle 2), les problèmes (2a) et (2b) font l'objet d'exercices. Une fois ces trois solutions partielles résolues, leur somme donne la solution complète que l'on peut écrire

Problème au limites d'un lac alimentant une nappe avec écoulement (modèle 3)

$$ \phi(x, z) = \sum_{n=1}^{\infty} \phi_{n,0}(x, z) + \phi_{n,1}(x, z) + \phi_{n,2}(x, z), \qquad (11.39) $$

$$ \phi_{n,0}(x, z) = A_n \sinh(\alpha_n(x - L)) \cos(\alpha_n z), \qquad (11.40) $$

$$ \phi_{n,1}(x, z) = B_n \cosh(\beta_n z) \sin(\beta_n x), \qquad (11.41) $$

$$ \phi_{n,2}(x, z) = C_n \sinh(\alpha_n x) \cos(\alpha_n z), \qquad (11.42) $$

$$ \alpha_n = \frac{(2n - 1)\pi}{2H}, \qquad (11.43) $$

$$ \beta_n = \frac{(2n - 1)\pi}{L}, \qquad (11.44) $$

$$ A_n = \frac{(-1)^{n-1} 4\Phi_0}{(2n-1)\pi \sinh(-\alpha_n L)}, \qquad (11.45) $$

$$ B_n = \frac{4\Phi_1}{(2n-1)\pi \cosh\left(\frac{(2n-1)\pi H}{L}\right)}, \qquad (11.46) $$

$$ C_n = \frac{(-1)^{n-1} 4\Phi_2}{(2n-1)\pi \sinh(\alpha_n L)}. \qquad (11.47) $$

### 11.4.3 Solution numérique

Il suffit de reprendre le script utilisé pour le Problème aux limites précédent et de changer les conditions aux limites des bords gauches et droits. Evidemment une fois qu'on a compris le principe de résolution numérique il n'est pas difficile de voir pourquoi, avec l'arrivée d'ordinateur toujours plus puissant, les solutions analytiques des problèmes aux limites sont devenues rares même quand elles existent et c'est dommage ! Réfléchissez pourquoi.

1. \# Dirichlet condition on the upper limits

2. $\phi_i[\text{Ny}-1,:] = 15$

3. $\phi_i[:0] = 20$

4. $\phi_i[:N_x-1] = 10$

5. \# Neumann condition flux=0 on the lower limit

6. $d\phi_i = \text{np(zeros}(N_x)$
<figure><img src="image_50.png" /><figcaption>FIGURE 11.7 – Solution du problème au limite des interactions entre un lac et une nappe dans laquelle a lieu un écoulement.</figcaption></figure>

### 11.4.4 Analyse

En fonction des valeurs que l'on impose aux charges latérales on voit immédiatement que le lac peut soit être alimenté des deux cotés soit être alimenté d'un côté et alimenter l'autre. Toutes ces situations se rencontrent évidemment dans la nature.

Une question que l'on se pose quand on regarde les figures 11.2, 11.5 ou 11.7; quid de l'environnement ? Comment ces lignes de courant et ses équipotentielles se raccordent-elles avec la nappe autour du lac ? Répondre à cette question dépend de la nature confinée ou non des écoulements autour du lac. En effet comme nous allons le voir dans le chapitre suivant les écoulements non confinés présentant une surface libre leur résolution est plus complexe à mettre en oeuvre.

## 11.5 Le problème du Barrage

Ce problème fait l'objet d'un exercice. On peut chercher une solution analytique mais, en raison notamment du changement de conditions sur les limites supérieures, cette recherche dépasse largement le cadre de ce cours. Par contre une solution numérique est relativement simple à obtenir.

<figure><img src="image_51.png" /><figcaption>FIGURE 11.8 – Problème de l'écoulement sous un barrage séparant deux masses d'eau (lacs)</figcaption></figure>

Le problème à résoudre est le suivant (figure 11.8) : un barrage sépare deux masses d'eau dont les surfaces libres sont situées à deux altitudes différentes correspondant à deux charges dans la nappe $Φ_1 > Φ_2$. Pour simplifier le problème et ne pas se préoccuper d'un éventuel écoulement régional on pose des flux nuls sur les
cotés et le fond de la nappe. On a donc une nappe qui connecte deux lacs séparés par un barrage. La condition à la limite sous le barrage est elle aussi une condition de flux nul : le fond du barrage constitue lui une ligne de courant.

La seule difficulté de la résolution consiste à ajouter une par une les conditions de flux (Neumann) sur les différents bords. Le résultat est obtenu figure 11.9

<figure><img src="image_52.png" /><figcaption>FIGURE 11.9 – Écoulement sous un barrage séparant deux masses d'eau (lacs)</figcaption></figure>

# 11.6 Exercices

## 11.6.1 Problème 2a

Résoudre et discuter le problème au limites suivant. Représenter graphiquement le résultat au moyen d'un script python.

$$ \Delta\phi = 0 \qquad (11.48) $$

$$ \phi(x, H) = 0 \quad , \quad \frac{\partial\phi}{\partial y}(x, 0) = 0 \qquad (11.49) $$

$$ \phi(0, y) = h_0 \quad , \quad \phi(L, y) = 0 \qquad (11.50) $$

## 11.6.2 Problème 2b

Résoudre et discuter le problème au limites suivant. Représenter graphiquement le résultat au moyen d'un script python.

$$ \Delta\phi = 0 \qquad (11.51) $$

$$ \phi(x, H) = 0 \quad , \quad \frac{\partial\phi}{\partial y}(x, 0) = 0 \qquad (11.52) $$

$$ \phi(0, y) = 0 \quad , \quad \phi(L, y) = h_2 \qquad (11.53) $$

## 11.6.3 Solution numérique

Reprendre le problème aux limites 1 des échanges entre un lac et sa nappe.

### 1. Convergence

(a) Établir un critère permettant de mesurer l'écart entre solution numérique et analytique

(b) Faire varier le critère de convergence et représenter l'évolution de l'écart défini précédemment avec le critère de convergence. Que constatez-vous ?

### 2. Vitesse de convergence

(a) faire varier la taille du domaine d'intégration et représenter le nombre d'itérations nécessaires pour atteindre le critère de convergence en fonctin du nombre de noeuds de la grille. Que constatez-vous ?
11.6.4 Le Barrage

Proposer une solution numérique du problème aux Limites défini par la figure 11.8.
## 12.1 Stockage de Fontsante

Il est envisagé de construire un centre de stockage de déchets d'incinération d'ordures ménagères dans le massif de l'Esterel (Var) à proximité du site de l'ancienne mine de fluorine de Fontsante 12.1 et 12.2. Le vallon qu'on envisage d'utiliser se situe le long de la ligne de partage des eaux du Massif : d'un côté le versant ouest qui est drainé vers la retenue de St-Cassien (laquelle constitue une ressource pour la fourniture en eau potable de l'Est Varois), et de l'autre le versant est qui s'écoule vers la Sioule et la mer, sans prélèvements d'eau potable.

Le site choisi dit du "Filon blanc" est drainé en surface vers l'Est et non vers la retenue de St-Cassien qui doit être protégée de tout danger de pollution. L'étude se propose de savoir ce qu'il en est des écoulements souterrains.

### 12.1.1 Écoulements

Le sol est constitué de schistes très peu perméables en profondeur mais diaclasés et altérés en surface. La perméabilité est plus importante dans l'axe des vallons, et plus faible dans les versants. Elle est estimée en moyenne à $K = 10^{-7}m/s$ sur les 20 premier mètres, et $K = 10^{-10}m/s$ en dessous. On a foré 8 piézomètres dont les positions et les niveaux moyens mesurés sont donnés dans la carte jointe 12.2. Deux mesures ont été réalisées donnant à peu de chose près les mêmes valeurs.

1. Quelle est l'échelle de la carte ?

2. Tracez approximativement la carte piézométrique du vallon du Filon Blanc. On prendra par exemple une isopieze tous les 5m. Essayez de distinguer les isopièzes dont vous pensez qu'elles sont probables et celles qui sont extrapolées.

3. En quoi le drainage de surface (limites de versants en blanc et thalwegs en bleu) constitue-t-il une aide ?

4. En déduire les lignes de courant et la direction de l'écoulement.

5. Quelle approximation faîtes vous pour répondre à la question précédente ?

6. Le site présente-t-il un danger potentiel pour la retenue de St-Cassien ?

7. Où préconiseriez vous le forage de deux nouveaux piézomètres ?

8. Estimez, en justifiant votre choix, le flux d'eau souterraine susceptible de s'écouler sous le site de Stockage

### 12.1.2 Pollution par solubilisation

On estime qu'une lame d'eau annuelle de 10mm est susceptible de s'infiltrer à travers le site de stockage (recouvert de terre compactée). En s'infiltrant on estime que l'eau se charge d'éléments en solution à hauteur de 1g/l. La porosité du sous-sol à l'aplomb du site est de l'ordre de 10%

1. Les éléments solubilisés ont une diffusivité $D = 1500 \mu m^2/s$, quel est le temps caractéristique de propagation du polluant sur toute l'épaisseur de la nappe sous la zone de stockage ?

2. Quel est le nombre de Péclet de l'écoulement de polluant ?

3. Y a-t-il un mode de transport dominant ?

4. Quel est les temps caractéristique de propagation du polluant à l'horizontal sur une distance équivalente à l'épaisseur de stockage ?

5. Le gradient de pente nous donne une indication sur le rapport entre vitesse verticale (que nous avons négligée) et la vitesse horizontale des écoulements. Si on considère que le polluant est advecté verticalement par l'eau qui s'écoule sous la zone de stockage quelle est l'échelle de temps correspondante ?
<figure><img src="image_53.png" /></figure>

© IGN 2012 – www.geoportail.gouv.fr/mentions-legales

Longitude : 6° 49' 19.9" E
Latitude : 43° 32' 55.1" N

<caption>FIGURE 12.1 – Carte régionale</caption>

6. Comment ces échelles de temps se comparent-elles ?

7. Quelles conclusions pouvez-vous tirer de ces calculs ?

8. Si vous deviez résoudre ce problème de pollution quel type d'équation serait le plus adapté ?

## 12.2 Etude de l'alimentation de la Loire entre Giens et Blois

La nappe de Beauce est une grande nappe constituée par des formations calcaires (Calcaires de Pithiviers, d'Etampes, de Brie et de l'Orléanais). Les transmissivités moyennes mesurées dans les calcaires de Pithiviers et d'Etampes, sur lesquels coulent la Loire, sont de l'ordre de (T ~ 5 · 10⁻²m²/s). La carte représente les hauteur piézométriques mesurées en 2002 en période de hautes eaux.

1. Tracer la (ou les) limite(s) de partage des eaux souterraines.

2. Tracer quelques lignes de courants montrant les grandes directions d'écoulement.

3. Estimer et représenter le vecteur gradient de charge tous les 10 km entre Blois et 10km au NW de Sully sur Loire.

4. Estimer le débit moyen par unité de largeur apporté par la nappe à la Loire pour chacune de vos mesures.

5. Intégrer ce débit et calculer l'apport total entre Giens et Blois.

6. le débit moyen mensuel de la Loire, calculé de décembre à février, est de 526 m³/s à Giens et de 536 m³/s à Blois. Quelle est la contribution directe de la nappe de la Beauce à l'augmentation de débit de la Loire ?

7. En fait la contribution de la nappe est estimée à environ 26 m³/s. D'où peut provenir l'écart constaté.

## 12.3 Écoulement dans les sables de Fontainebleau au niveau du plateau de Saclay

La figure 12.4 représente la nappe des sables de Fontainebleau au niveau du plateau de Saclay dans les Yvelines. Les sites 1 et 2 sont deux sites de recherches nucléaires. L'objectif de ce travail est de comprendre la façon donc l'eau s'écoule entre les deux rivières qui entaillent le plateau, la Bièvre au nord et l'Yvette au sud, et d'appréhender le devenir d'une éventuelle pollution de l'un des sites de recherche.

Nous allons tout d'abord compléter la carte. Dans un second temps nous nous intéresserons à l'écoulement en plan le long d'une coupe dans la nappe. Enfin nous essaierons de comprendre ce qui se passerait en cas de pollution accidentelle.
<figure><img src="image_54.png" /><figcaption>FIGURE 12.2 – Site d'étude</figcaption></figure>
<figure><img src="image_55.png" /><figcaption>FIGURE 12.3 – Carte piézométrique de la nappe de Beauce (2002), source des données : [17]</figcaption></figure>
<figure><img src="image_56.png" /><figcaption>FIGURE 12.4 – Carte hydrogéologique du plateau de Saclay. Source : [18]</figcaption></figure>
### 12.3.1 Analyse de la carte

1. Tracer les lignes de courant partant des points $X_1, X_2, X_3, X_4, X_5$ et compléter la légende

2. Tracer la ligne de partage des eaux de la nappe de les sables de Fontainebleau entre la Bièvre et l'Yvette. Idem compléter la légende

### 12.3.2 Etude des écoulements le long d'une coupe

1. Représenter le profil des hauteurs piézométriques de la nappe le long du segment de droite [AB].

2. Représenter la position de la Bièvre et de l'Yvette.

3. Tracer la droite passant par les deux rivières. A partir d'arguments géologiques, montrer que cette droite correspond, dans l'hypothèse la plus simple, au mur de la nappe des sables de Fontainebleau

4. Calculer le gradient hydraulique moyen entre le site 2 et la Bièvre.

5. En supposant la perméabilité de la formation égale à $10^{-5}$ m/s calculer la vitesse de Darcy U moyenne sur le segment

6. Calculer le débit par unité de largeur $q$

7. Pour une porosité $\omega$ de ~ 15%, calculer le volume d'eau $V$ contenu dans la zone délimitée par le rectangle CDEF

8. Expliquer pourquoi le flux entrant ou sortant de ce volume par les côtés DE et CF est approximativement nul.

9. Si la nappe est en régime permanent que vaut le flux $Q$ entrant par CD et sortant par EF ?

### 12.3.3 Pollution par infiltration

Le site 2 pollue la nappe en laissant s'infiltrer un flux massique $Q_p$ par unité de temps de fluide fortement radioactif. On néglige les phénomènes d'advection et de diffusion et on considère le volume de nappe compris dans le rectangle CDEF comme un réservoir dont les eaux se mélangent instantanément. On considère en outre que la concentration en fluide polluant est nulle dans les eaux de la nappe situées à l'extérieur du rectangle CDEF.

1. Par où le polluant est-il évacué ? Quelle rivière est susceptible d'être touchée par la pollution ?

2. Écrire l'équation de conservation de la masse $M$ du polluant et la résoudre. On considérera qu'il n'y a ni adsorption ni sédimentation.

3. Quelle est l'échelle de temps $\tau$ problème ? De quels paramètres dépend-elle ?

## 12.4 Alimentation des lacs du Médoc

Entre l'estuaire de la Gironde et l'océan atlantique, deux lacs, les lacs d'Hourtin au nord et de Lacanau au sud, sont alimentés par une nappe libre. Celle-ci s'écoule dans une formation aquifère constituée des sables des Landes, des sables de Castets et d'anciennes dunes quaternaires proches de l'océan. La figure 12.5 représente la carte des isopièzes de la nappe. On se propose ici de s'intéresser à la dynamique d'alimentation de ces deux lacs d'eau douce situés à proximité de l'Océan.

### 12.4.1 Hydrogéologie des formations sableuses

Aux endroits où les données manquent on pourra faire des hypothèses raisonnables afin de remédier au manque de données.

1. Représenter sur la figure 12.5 le bassin d'alimentation des lacs du médoc. On tracera la limite de partage des eaux de la nappe séparant la région qui alimente les lacs des régions qui alimentent la Gironde, l'océan ou la bassin d'Arcachon plus au sud.

2. Représenter quelques lignes de courant qui illustrent les mouvements d'eau dans la nappe

3. Tracer les coupes AB et CD indiquant les hauteurs d'eau dans la nappe sur les figures 12.6 et 12.7. On débutera la coupe à la côte que l'on positionnera au point (0,0)

4. Pour chaque coupe :

(a) Représentez les lignes de partage des eaux par des lignes verticales.

(b) Calculez les gradients hydrauliques moyens dans chaque région ainsi délimitée.

(c) Indiquez par un vecteur les directions d'écoulement.
(d) Dans la suite on fera l'hypothèse que les écoulements sont horizontaux dans le bassin d'alimentation des lacs. Cette hypothèse vous semble-t-elle raisonnable ?

5. plusieurs mesures ont permis de déterminer une transmissivité $T \approx 3.3 \cdot 10^{-3} \text{ m}^2/\text{s}$. Estimez à partir de cette valeur et de la loi de Darcy le débit d'alimentation de chacun des lacs.

### 12.4.2 Biseau salé

Les lacs sont remplis d'eau douce et alimentés par la formation décrite précédemment. Sous cette formation aquifère se trouvent un substratum de conductivité beaucoup plus faible. À la limite continent-océan, le substratum se trouve à une profondeur de $z = 50 \text{ m}$ sous le niveau de la mer. Sa limite remonte ensuite doucement pour atteindre une profondeur de $z = 0 \text{ m}$, 15 km à l'Est de l'océan. La présence de l'océan à proximité induit la formation d'un biseau salé. La densité de l'eau de mer est $\rho_s = 1025 \text{ kg/m}^3$, celle de l'eau douce est $\rho = 1000 \text{ kg/m}^3$.

1. Représentez le profil du biseau sur la Coupe A-B entre l'océan et la rive Est du lac d'Hourtin (on considèrera que le biseau s'arrête quand il touche une formation imperméable)

2. De l'eau salée se trouve-t-elle sous le lac d'Hourtin ?

3. Si oui, à quelle profondeur se situe l'eau de mer sous le lac ?

4. De l'eau de mer est-elle susceptible de polluer le lac ?

## 12.5 Alimentation des lacs de Biscarosse

Dans les Landes au sud du bassin d'Arcachon, les lacs de Cazeaux, Parentis (Biscarosse) et d'Aureillan forment un système connecté alimenté par une nappe libre. Celle-ci s'écoule dans une formation aquifère constituée des sables des Landes, des sables de Castets et d'anciennes dunes quaternaires proches de l'océan. La figure 12.8 représente la carte des isopièzes de la nappe. On se propose ici de s'intéresser à la dynamique d'alimentation de ces lacs d'eau douce situés à proximité de l'Océan.

Aux endroits où les données manquent on pourra faire des hypothèses raisonnables afin de remédier au manque de données.

1. Représenter sur la figure 12.8 le bassin d'alimentation des lacs. On tracera la limite de partage des eaux de la nappe séparant la région qui alimente les lacs des régions qui ne les alimentent pas.

2. Représenter quelques lignes de courant qui illustrent les mouvements d'eau dans la nappe

3. Tracer les coupes AB et CD indiquant les hauteurs d'eau dans la nappe sur les figure 12.9 et 12.10. On débutera la coupe à la côte que l'on positionnera au point (0,0)

4. Pour chaque coupe :

(a) Représentez les lignes de partage des eaux par des lignes verticales.

(b) Calculez les gradients hydrauliques moyens dans chaque région ainsi délimitée.

(c) Indiquez par un vecteur les directions d'écoulement.

(d) Dans la suite on fera l'hypothèse que les écoulements sont horizontaux dans le bassin d'alimentation des lacs. Cette hypothèse vous semble-t-elle raisonnable ?

5. Des mesures ont permis de déterminer une transmissivité $T \approx 6 \cdot 10^{-4} \text{ m}^2/\text{s}$. Estimez à partir de cette valeur et de la loi de Darcy le débit d'alimentation de chacun des lacs.

## 12.6 Nappe de Champigny

Au niveau de la brie, les calcaires dits de Champigny regroupe des formations qui abritent une nappe libre importante pour l'alimentation en eau de l'Île de France.

1. Compléter la légende

2. Entourer les flèches qui vous semblent, de façon évidente, mal dessinées. Expliquer pourquoi.

3. En considérant une transmissivité de $10^{-3} \text{ m}^2/\text{s}$ et en détaillant la procédure de calcul, trouver le débit total de la nappe, dans sa partie alimentant la Seine, au travers de l'isopièze 80m.
<figure><img src="image_57.png" /><figcaption>FIGURE 12.5 – Image Landat 8 de la région des lacs médocains sur laquelle sont superposés les isopièzes de la nappe des sables des landes, de Castets et des sables dunaires. La grille de coordonnées des latitudes et longitudes est donnée en mètres. Les cotes sont données en mètres au dessu du niveau de la mer.sourceforge.net données : [5] et SIGES Aquitaine.</figcaption></figure>
<figure><img src="image_58.png" /><figcaption>FIGURE 12.6 – Coupe A-B. La ligne rouge représente la limite du substratum considéré comme imperméable</figcaption></figure>

<figure><img src="image_59.png" /><figcaption>FIGURE 12.7 – Coupe B-C</figcaption></figure>
<figure><img src="image_60.png" /><figcaption>FIGURE 12.8 – Image Landat 8 de la région de Biscarosse sur laquelle sont superposés les isopièzes de la nappe des sables des landes, de Castets et des sables dunaires. La grille de coordonnées des latitudes et longitudes est donnée en mètres. Les cotes sont données en mètres au dessus du niveau de la mer. [5] et SIGES Aquitaine.</figcaption></figure>
<figure><img src="image_61.png" /><figcaption>FIGURE 12.9 – Coupe A-B</figcaption></figure>

<figure><img src="image_62.png" /><figcaption>FIGURE 12.10 – Coupe C-D</figcaption></figure>
<figure><img src="image_63.png" /><figcaption>FIGURE 12.11 – Nappe des calcaires de Champigny</figcaption></figure>
Troisième partie

Solution partielle des exercices
## Étude d'un sol (§ 1.5.1, pp 11)

1. Marne sablo graveleuse

2. Calculez

la masse volumique du sol

$$ \rho_s = \frac{M_s}{V_s} = \frac{302}{196} = 1.54 \text{ g/cm}^3 = 1540 \text{ kg/m}^3 \qquad (12.1) $$

sa porosité

$$ \phi = 1 - \frac{\rho_s}{\rho_m} = 1 - \frac{1540}{2650} = 0.419 \qquad (12.2) $$

sa teneur en eau

$$ \theta = (0.3318 - 0.3024)/(1000 * 1.96^{-4}) = 0.150 \qquad (12.3) $$

son pourcentage de saturation

$$ S = 0.150/0.419 = 0.358 \qquad (12.4) $$

## Porosité (§ 1.5.2, pp 11)

1.

$$ V_p = V_{tot} - V_g = d^3 \left(1 - \frac{\pi}{6}\right). \qquad (12.5) $$

$$ S_p = 4\pi R^2 = \pi d^2. \qquad (12.6) $$

2.

$$ V_c \sim \pi d^2 \qquad (12.7) $$

c'est une valeur approximative et maximum car, comme les grains sont jointifs, la couche ne peut être uniforme sur toute la surface des grains. On peut retrouver le résultat en approximant le volume compris entre deux coques de diamètre $d$ et $d + \epsilon$ avec ici $\epsilon = 2$ (attention pas 1 car on parle en diamètre !)

$$ V_c = \Delta V = \frac{\pi}{6} [(d + \epsilon)^3 - d^3] \qquad (12.8) $$

3.

$$ V_f = V_p - V_c \sim d^3 \left(1 - \frac{\pi}{6}\right) - \pi d^2 \qquad (12.9) $$

4.

$$ \omega_t = \frac{V_p}{V_t} = 1 - \frac{\pi}{6} \qquad (12.10) $$

$$ \omega_e = \frac{V_f}{V_t} = \left(1 - \frac{\pi}{6}\right) - \frac{\pi}{d} \qquad (12.11) $$

5. figure 12.12

## Piézométrie (§ 2.6.1, pp 18)

1. (a)

$$ \phi_1 = h_1, \phi_2 = h_2 + \frac{\rho_1}{\rho_2}(h_1 - h_2) $$

(b)

$$ P_{tot} = g(\rho_1 h_1 + \rho_2 (h_2 - h_1)) $$

2.

$$ \phi_1 = h_1, \phi_2 = h_2 + \frac{\rho_1}{\rho_2}(h_1 - h_2), \phi_3 = h_3 + \frac{\rho_2}{\rho_3}(h_2 - h_3) + \frac{\rho_1}{\rho_3}(h_1 - h_2) $$
<figure><img src="image_64.png" /><figcaption>FIGURE 12.12 – Porosités totale et efficace d'un sol uniforme</figcaption></figure>

## Nappe côtière, biseau salé (§ 2.6.3, pp 19)

1. $$p_{ed} = \rho_{ed}g(h + z) \qquad (12.12)$$

$$p_{es} = \rho_{es}gz \qquad (12.13)$$

2. $$z = h \frac{\rho_{ed}}{\rho_{es} - \rho_{ed}} = 40h \qquad (12.14)$$

3. $$z = 80m, \tan \theta = 80/200 \rightarrow \theta \approx 22^{\circ} \qquad (12.15)$$

4. Le rabattement maximum est de $1.5m$

## Écoulements simple (1D) (§ 3.6.1, pp 23)

1. (a) $$Q = UA = 0.4 \, \text{m}^3/\text{j} \qquad (12.16)$$

(b) $$U = K \frac{\Delta h}{L} \sim 20 \, \text{m/j} \qquad (12.17)$$

(c) $$v = \frac{U}{n} = \frac{20}{0.36} = 55.5 \, \text{m/j} \qquad (12.18)$$

2. $$K = \frac{Vn}{J} = 75 \, \text{m/j} \qquad (12.19)$$

## Grains carrés ! (§ 3.6.2, pp 23)

$$U = u \frac{S_p}{S_t} \qquad (12.20)$$

$$U = u \left( 1 - 53 \left( \frac{a}{A} \right)^2 \right) \qquad (12.21)$$
## Bilan global et temps de résidence (§ 5.6, pp 36)

1. La surface des océans est donnée par

$$S_{oc} \sim 0.8(4\pi R^2) \sim 4.05 \cdot 10^8 \text{ km}^2. \qquad (12.22)$$

En effectuant le bilan des flux et en divisant par cette surface on obtient $P_{oc} = E_{oc} \approx 1\text{m/a}$ le bilan est bouclé.

2. $P_c = V_c/A_c = 0.72 \text{ m/a}$, $E_c = 0.48\text{m/a}$. Le bilan n'est pas bouclé.

3. Soit le bilan n'est pas nul et alors il y a variation du stock d'eau dans un compartiment, soit il manque un flux sortant ce sont les écoulements de surface et souterrains. Le bilan semble indiquer que c'est la seconde solution qui est la bonne

4. Le temps caractéristique ici est obligatoirement donné par le rapport entre un flux et le volume d'un réservoir.

$$T = \frac{V}{Q} \qquad (12.23)$$

Comme les réservoirs semblent à l'équilibre, on prend le flux que l'on veut. en négligeant les variations de masse en fonction de la température on à $10^3 \text{ km}^3 \approx 10^{15} \text{ kg}$.

Pour l'atmosphère : les précipitation sont de $P = 505 \cdot 10^3 \text{ km}^3/\text{a}$ et le volume d'eau de $V = 15.5 \cdot 10^{15} \text{ kg}$.
On en déduit donc

$$T_{atm} = \frac{15.5}{505} = 0.03 \text{ a} \sim 11\text{j} \qquad (12.24)$$

$$T_{oc} = \frac{1.4 \cdot 10^6}{434} \approx 3200 \text{ a} \qquad (12.25)$$

$$T_{Cont} \approx \frac{43400 + 15300}{107} \approx 550 \text{ a} \qquad (12.26)$$

5. À quoi correspond ce temps ?

Ce temps correspond à un temps caractéristique de résidence ou de renouvellement du réservoir en imaginant que le volume entrant dans le réservoir à chaque instant remplace un volume équivalent qui quitte le réservoir.

6.

$$T_{Cont} \approx \frac{15300}{107} \approx 140 \text{ a} \qquad (12.27)$$

soit un temps de renouvellement près de 4 fois plus rapide qu'aujourd'hui.

7. Si on multiplie le flux de "runoff" par deux on a $R = 72 \cdot 10^3 \text{ km}^3/\text{a}$. Le bilan des flux entrants et sortants donne

$$P - E - R = -36 \cdot 10^3 \text{ km}^3/\text{a} = \frac{dV_{cont}}{dt} \qquad (12.28)$$

où $V_{cont}$ est le le volume d'eau des continents. En négligeant la biosphère et les rivières et lac on obtient après intégration

$$V_{cont} = -36 \cdot 10^3 t + 58 \cdot 10^6 \text{ km}^3 \qquad (12.29)$$

## Pollution ponctuelle d'un réservoir (§ 5.6, pp 36)

$$M(t) = M_0 \exp\left(-\frac{Q}{V}t\right) = 2000 \exp(-1.6 \cdot 10^{-6}t) \qquad (12.30)$$

On peut définir un temps caractéristique de la décroissance

$$\tau = \frac{V}{Q} \sim 623000 \text{ s} \sim 7 \text{ j} \qquad (12.31)$$
<figure><img src="image_65.png" /><figcaption>FIGURE 12.13 – Evolution de la saumure dans le lac en fonction du temps. les points marquent les temps correspondant à t = kτ</figcaption></figure>

Pollution d'un lac : régimes stationnaire et transitoire (§ 5.6, pp 36)

1. Comme on ne connaît pas le volume du lac on peut toujours écrire une équation différentielle représentant le bilan mais on ne peut la résoudre. On va donc se placer en régime stationnaire (quand le temps n'influe plus) pour résoudre ce problème.

$$C = \frac{QC_i + Q_M}{Q + v_s A} = 0.028 \text{ mg l}^{-1} \qquad (12.32)$$

2. On veut $C = C_i$ il suffit donc de remplacer dans les bilans et on obtient

$$Q_M = C_i v_s A = 0.25 \text{ g/s} \qquad (12.33)$$

3. (a) La grande différence avec ce qui précède est qu'ici on vous donne la profondeur donc vous pouvez calculer le volume $V = AH$.

$$M = AHC_e \left[ 1 - \exp \left( -\frac{Q}{AH} t \right) \right] \qquad (12.34)$$

Et

$$C = C_e \left[ 1 - \exp \left( -\frac{Q}{AH} t \right) \right] \qquad (12.35)$$

où $C_e$ correspond à la concentration en régime stationnaire sans tenir compte du terme de sédimention et en supposant nulle la masse initiale de phosphore. Numériquement on obtient

$$C = 0.08 \left[ 1 - \exp \left( -\frac{t}{620} \right) \right] \qquad (12.36)$$

$C$ est donné en mg/l et le temps est compté en jours.

(b)

$$C = \frac{QC_e}{Q + v_s A} \left[ 1 - \exp \left( -\frac{Q + v_s A}{AH} t \right) \right] \qquad (12.37)$$

Numériquement on trouve que le temps caractéristique est

$$\tau = \frac{AH}{Q + v_s A} \sim 230 \text{ j} \qquad (12.38)$$

L'évolution des deux cas est représentée sur la figure suivante. On atteint le régime stationnaire dans les deux cas au bout de ~ 5τ mais on peut voir l'influence considérable de la sédimentation à la fois sur la concentration d'équilibre (plus de deux fois plus faible) et sur le temps caractéristique d'évolution du système (près de trois fois plus court).
<figure><img src="image_66.png" /><figcaption>FIGURE 12.14 – Evolution de la concentration en phosphore du lac en fonction du temps. les points marquent les temps correspondant à t = kτ.</figcaption></figure>

4. (a) Sans sédimentation

On considère tout d'abord qu'au temps $t = 0$ du glissement de terrain l'équilibre précédent est établi (on a dépassé cinq fois le temps d'évolution du système). Donc à $t = 0$, $C = C_e$, $M_0 = V_0 C_e$. L'évolution du volume du lac, quand il n'y pas de sédimentation, est donnée par

$$V = V_0 + \Delta Qt \qquad (12.39)$$

où $\Delta Q = Q_i - Q_o$ la profondeur moyenne est donc

$$H = 10 + 10^{-7}t \qquad (12.40)$$

Soit en gros une hausse de l'ordre du m par an.

L'équation différentielle d'évolution de la masse est

$$\dot{M} + \frac{Q_o}{V_0 + \Delta Qt} M = Q_i C_e \qquad (12.41)$$

le coefficient devant $M$ dépend de $t$ on ne peut plus séparer les variables comme avant. On doit donc calculer un facteur intégrant soit

$$I = \exp \int \frac{Q_o}{V_0 + \Delta Qt} \qquad (12.42)$$

soit

$$I = \exp \left( \frac{Q_o}{\Delta Q} \ln(V_0 + \Delta Qt) \right) = (V_0 + \Delta Qt)^{\alpha} \qquad (12.43)$$

avec $\alpha = Q_o/\Delta Q$ l'équation différentielle devient alors

$$(I\dot{M}) = Q_i C_e (V_0 + \Delta Qt)^{\alpha} \qquad (12.44)$$

En intégrant des deux côtés on obtient

$$(V_0 + \Delta Qt)^{\alpha} M = \frac{Q_i C_e}{(\alpha + 1)\Delta Q} (V_0 + \Delta Qt)^{\alpha+1} + c \qquad (12.45)$$

où $c$ est une constante d'intégration que l'on trouve grâce à la condition initiale Pour avancer on remarque que $(\alpha + 1)\Delta Q = Q_i$. On a alors

$$M = C_e(V_0 + \Delta Qt) + \frac{c}{(V_0 + \Delta Qt)^{\alpha}} \qquad (12.46)$$

La condition initiale donne $c = 0$ ce permet d'obtenir la solution de notre problème. On voit bien l'importance de la condition initiale (ou des conditions aux limites pour les problèmes dans l'espace) car c'est elle qui conditionne la forme de la solution. Ici

$$M = C_e(V_0 + \Delta Qt). \qquad (12.47)$$
Et

$$C = C_e \qquad (12.48)$$

Tout cela pour ça me direz vous...Oui parce-que les conditions initiales sont très simples et égales au régime stationnaire.

## (b) Avec sédimentation

On considère tout d'abord qu'au temps $t = 0$ du glissement de terrain l'équilibre précédent est établi (on a dépassé cinq fois le temps d'évolution du système). Donc à $t = 0$, $C_0 = QiC_e/(Q_i + v_sA)$, $M_0 = V_0C_e$. On reprend l'équation de conservation en intégrant la sédimentation au numérateur du coefficient de la masse $M$ soit

$$\dot{M} + \left( \frac{Q_o + v_s A}{V_0 + \Delta Qt} \right) M = Q_i C_e \qquad (12.49)$$

Le facteur intégrant devient

$$I_s = \exp \left( \frac{Q_o + v_s A}{\Delta Q} \ln(V_0 + \Delta Qt) \right) = (V_0 + \Delta Qt)^{\alpha_s} \qquad (12.50)$$

où

$$\alpha_s = \frac{Q_o + v_s A}{\Delta Q} \qquad (12.51)$$

L'intégration mène toujours à

$$(V_0 + \Delta Qt)^{\alpha_s} M = \frac{Q_i C_e}{(\alpha_s + 1)\Delta Q} (V_0 + \Delta Qt)^{\alpha_s + 1} + c \qquad (12.52)$$

Soit

$$M = \frac{Q_i C_e}{(\alpha_s + 1)\Delta Q} (V_0 + \Delta Qt) + \frac{c}{(V_0 + \Delta Qt)^{\alpha_s}} \qquad (12.53)$$

Cette fois $(\alpha_s + 1)\Delta Q = v_sA + Q_i$ soit

$$M = \frac{Q_i C_e}{Q_i + v_s A} (V_0 + \Delta Qt) + \frac{c}{(V_0 + \Delta Qt)^{\alpha_s}} \qquad (12.54)$$

Ici la constante $c$ est toujours nulle du fait du changement des conditions initiales soit

$$C = \frac{Q_i C_e}{Q_i + v_s A} \qquad (12.55)$$

et

$$M = \frac{Q_i C_e}{Q_i + v_s A} (V_0 + \Delta Qt) \qquad (12.56)$$

Pour finir si on part d'une condition initiale différente $C_0$ on aura

$$C_0 V_0 = \frac{Q_i C_e}{Q_i + v_s A} V_0 + \frac{c}{V_0^{\alpha_s}} \qquad (12.57)$$

soit

$$c = V_0^{\alpha_s + 1} \left( C_0 - \frac{Q_i C_e}{Q_i + v_s A} \right) \qquad (12.58)$$

L'équation finale ressemble donc à

$$M = \frac{Q_i C_e}{Q_i + v_s A} (V_0 + \Delta Qt) + \left( C_0 - \frac{Q_i C_e}{Q_i + v_s A} \right) \frac{V_0^{\alpha_s + 1}}{(V_0 + \Delta Qt)^{\alpha_s}} \qquad (12.59)$$

Pour la concentration on aura

$$C = \frac{Q_i C_e}{Q_i + v_s A} + \left( C_0 - \frac{Q_i C_e}{Q_i + v_s A} \right) \left( \frac{V_0}{V_0 + \Delta Qt} \right)^{\alpha_s + 1} \qquad (12.60)$$

Même si elle est moins évidente à première vue, cette formule montre que l'on tend vers la concentration asymptotique en décroissant depuis la valeur de concentration initiale par une décroissance un peu plus compliquée qu'une simple décroissance exponentielle.
Variation de vitesse dans une nappe (§ 6.3.1, pp 41)

1.

$$U_1 = u_1 \omega_{eff} = 1.5 \cdot 10^{-5} \text{ m/s} \qquad (12.61)$$

2.

$$q_1 = U_1 h_1 = 5.25 \cdot 10^{-4} \text{ m}^2/\text{s} \qquad (12.62)$$

3.

$$U_2 = \frac{q_1}{h_2} = 2.1 \cdot 10^{-5} \text{ m/s,} \qquad (12.63)$$

et

$$u_2 = \frac{U_2}{\omega_{eff}} = 1.4 \cdot 10^{-4} \text{ m/s} \qquad (12.64)$$

Champs de vecteurs et champs de vitesse (§ 6.3.2, pp 41)

<figure><img src="image_67.png" /><figcaption>FIGURE 12.15 – Écoulements solutions de l'équation de conservation de la masse.</figcaption></figure>

Le principe est le même dans tous les cas. On teste la condition

$$\mathrm{div} \mathbf{U} = \left( \frac{\partial}{\partial x} \mathbf{i} + \frac{\partial}{\partial y} \mathbf{j} \right) \cdot \mathbf{U} = 0 \qquad (12.65)$$

$$x\mathbf{i} - y\mathbf{j} \qquad (12.66)$$

1 - 1 = 0 donc OUI

$$x^2\mathbf{i} - y^2\mathbf{j} \qquad (12.67)$$

2x - 2y = 0 quels que soient x et y donc NON

$$y\mathbf{i} + x\mathbf{j} \qquad (12.68)$$

0 + 0 = 0 donc OUI

$$y\mathbf{i} - x\mathbf{j} \qquad (12.69)$$

idem 0 - 0 = 0 OUI

$$y^2\mathbf{i} + x^2\mathbf{j} \qquad (12.70)$$

idem 0 + 0 = 0 OUI
<figure><img src="image_68.png" /><figcaption>FIGURE 12.16 -</figcaption></figure>

4. $$U_1 = -K_1 \nabla H_1 \quad (12.83)$$

idem pour U2

5. $$Q_{1,x} = -K_1 \frac{dH_1}{dx} = -K_2 \frac{dH_2}{dx} = Q_{2,x} \quad (12.84)$$

6. $$\tan \alpha = \tan \beta \frac{K_1}{K_2} \quad (12.85)$$

## Calcul d'un vecteur vitesse (§ 8.4.1, pp 49)

figure 12.16

## Calcul de vitesse dans une nappe captive (§ 8.4.2, pp 49)

On applique la technique vue précédemment la fig 12.17 donne la solution

## Propriétés de la fonction courant (§ 8.4.3, pp 50)

1. Utiliser l'identité vectorielle suivante : $$\nabla \times \nabla f = 0, \quad (12.86)$$

le rotationnel d'un gradient d'une fonction $f$ est nul lui aussi.

2. Dans un écoulement plan on relie deux points de deux lignes de courant par une courbe quelconque. un élément de longueur de cette courbe s'écrit $d\mathbf{S} = d\mathbf{x}i + d\mathbf{y}j$ et le débit passant au travers de ds sera donné par la norme de $\mathbf{U} \times d\mathbf{S} = \|\mathbf{U}\|\|d\mathbf{S}\| \sin(\alpha)$ où $\alpha$ est l'angle que font les deux vecteurs. Les calculs suivent...

3. Le fait que c'est une condition suffisante a déjà été vu dans le chapitre de cours correspondant. pour montrer que c'est nécessaire on prend un contre exemple.

## Étude d'écoulements (§ 8.4.4, pp 51)

1. Premier cas : $\mathbf{U} = \mathbf{x}\mathbf{i} - \mathbf{y}\mathbf{j}$.

$$\psi(x, y) = xy \quad (12.87)$$

Second cas : $\mathbf{U} = \mathbf{y}\mathbf{i} + \mathbf{x}\mathbf{j}$

$$\psi(x, y) = \frac{1}{2}(y^2 - x^2) \quad (12.88)$$
2. Il suffit de vérifier que le laplacien n'est pas nul.

3. Idem le laplacien n'est pas nul.

4. (a) Les dérivées second s'annulent donc $\Delta\phi = 0$

(b)

$$ \psi = K(bx - ay) \qquad (12.89) $$

<figure><img src="image_69.png" /></figure>

<figure><img src="image_70.png" /></figure>

(d) $Q = 1 \text{ m}^2/\text{s}$

5. (a) On dérive $\psi$ soit

$$ \partial_x\psi = -U_y = 2bx \text{ et } \partial_y\psi = U_x = -2by \text{ d'où} $$

$$ U = -2byi - 2bxj \qquad (12.90) $$

(b) on intègre la vitesse pour obtenir $\phi$ cette fois soit $U_x = -2by = -K\partial_x\phi$ d'où

$$ \phi = \frac{2b}{K}yx \qquad (12.91) $$

On vérifie avec l'autre composante que c'est bon! Comme toujours les résultats sont définis à une constante près.

(c)

<figure><img src="image_71.png" /></figure>

6. (a) Oui.

(b)

$$ \psi(x, y) = -3Kx^2y + Ky^3 \qquad (12.92) $$
<figure><img src="image_72.png" /></figure>

## Écoulements autour d'un puit (§ 9.6.1, pp 58)

1.

$$ T = \frac{Q}{2\pi r J} = 32 \, \text{m}^2/\text{h} \qquad (12.93) $$

La Transmissivité permet de ne pas avoir à connaître la profondeur et la conductivité de la nappe.

2.

$$ U = K \nabla h = 1.25 \cdot 10^{-5} \, \text{m/s} \qquad (12.94) $$

$$ Q = (2\pi r D)U = 0.047 \, \text{m}^3/\text{s} \sim 170 \, \text{m}^3/\text{j} \qquad (12.95) $$

3. La position du cylindre sur lequel on effectue le bilan de masse est $r = 0.5(20 + 150) = 85 \, \text{m}$ et la conductivité $K$ est donnée par

$$ K = \frac{Q}{2\pi r D \nabla h} = 5.4 \cdot 10^{-4} \, \text{m/s} \qquad (12.96) $$

## Alimentation en eau potable (§ 9.6.2, pp 58)

1. La consommation annuelle $C_a$ est obtenue en effectuant le produit de la consommation journalière par le nombre d'habitant et ce sur une année soit

$$ C_a = 15000 \times 0.375 \times 365.25 \approx 2 \cdot 10^6 \, \text{m}^3/\text{a}. \qquad (12.97) $$

2. Le débit $Q$ est donné par la consommation journalière divisée par la durée d'une journée soit

$$ Q = \frac{15000.375}{86400} \approx 6.5 \cdot 10^{-2} \, \text{m}^3/\text{s}. \qquad (12.98) $$

3. Le rabattement est calculé au puit soit à une distance $r = 0.5 \, \text{m}$ de l'origine. Le rayon d'influence $R = 500 \, \text{m}$ sert de condition aux limites. On sait que la charge y est $h_R = 30 \, \text{m}$. La nappe étant libre, on se sert de l'équation donnant la charge d'une nappe libre en fonction du rayon (??). Pour répondre à la question on doit d'abord calculer la charge en $r, h_r$, et en déduire la différence

$$ h_R - h_r \approx 9.4 \, \text{m} \qquad (12.99) $$

4. Afin d'évaluer la vitesse de Darcy il faut calculer le gradient de charge. On le fait grâce à l'équation donnant la charge dans la nappe (??) que l'on dérive.

$$ U = \frac{Q}{\pi} \frac{1}{2r\sqrt{h_R^2 - \frac{Q}{\pi K} \ln \frac{R}{r}}} \approx 0.001 \, \text{m/s} \qquad (12.100) $$

5. Pour le savoir il suffit de calculer le débit et le rabattement correspondant qui doit rester inférieur à 30m.

$$ h_R - h_r \approx 15.6 \, \text{m} \qquad (12.101) $$

Le pompage nécessaire dans la nappe permet donc bien de couvrir la croissance de la municipalité.
## Problème 2a (§ 11.6.1, pp 73)

### PL

$$ \Delta\phi = 0 \qquad (12.102) $$

$$ \phi(x, H) = 0 \quad , \quad \frac{\partial\phi}{\partial y}(x, 0) = 0 \qquad (12.103) $$

$$ \phi(0, y) = \Phi_0 \quad , \quad \phi(L, y) = 0 \qquad (12.104) $$

Pour la représentation graphique on prendra $\Phi_0 = 20m$

### Solution

$$ \phi = \sum_{n=1}^{\infty} A_n \sinh(\alpha_n(x - L)) \cos(\alpha_n y) \qquad (12.105) $$

$$ \alpha_n = \frac{(2n - 1)\pi}{2H} \qquad (12.106) $$

$$ A_n = \frac{(-1)^{n-1} 4\Phi_0}{(2n - 1)\pi \sinh(-\alpha_n L)} \qquad (12.107) $$

### Script

```python
def solution_ecoulement_gauche():
    # Length L, Depth H, Number of terms in the fourier series N
    L=200
    H=100
    N=100

    # Initialize Matrices
    xi= np.linspace(0,L,100)
    yi= np.linspace(0,H,100)

    x,y = np.meshgrid(xi,yi)
    phi=0*x

    # Non zero BC
    h0=20

    # Initialize fourier constants
    c= np(zeros(N))

    # Loop to calculate Phi
    for i in range(N):
        alpha = (2*i+1)*np.pi/(2*H)
        c[i] = (-1)**(i) * 4* h0/((2*i+1) * np.pi * np.sinh(-alpha*L))
        print(2*i+1, c[i])
        phi += c[i]*np.sinh(alpha*(x-L))*np.cos(alpha*y)

    # Calculate velocity (here K=1 m/s)
    U=-np.diff/np diff(x, axis=1)
    V=-np diff(, axis=0)/np diff(y, axis=0)

    # Plot the result
    fig = plt.figure(figsize=(14/2.54,7/2.54))
    ax = fig.add_subplot(111)
    # bins = np. arange(h0)+1
    CS = ax. contour(x,y,phi, levels=[1,2,4,6,8,10,12,14,16,18], colors='blue')
```
39

ax.clabel(CS, inline=1, fontsize=10)

40

ax.streamplot(x[1:,1:],y[1:,1:],U[1:,:],V[:1:],color='red',density=0.5)

41

ax.axis('square')

42

ax.set_xlim(0,L)

43

ax.set_ylim(0,H)

44

45

ax.set_xlabel('Distance')

46

ax.set_ylabel('Depth')

47

48

plt.show()

## Problème 2b (§ 11.6.1, pp 73)

### PL

$$ \Delta\phi = 0 \qquad (12.108) $$

$$ \phi(x, H) = 0 \quad , \quad \frac{\partial\phi}{\partial y}(x, 0) = 0 \qquad (12.109) $$

$$ \phi(0, y) = 0 \quad , \quad \phi(L, y) = \Phi_{i2} \qquad (12.110) $$

Pour la représentation graphique on prendra $\Phi_2 = 10m$

### Solution

$$ \phi = \sum_{n=1}^{\infty} A_n \sinh(\alpha_n x) \cos(\alpha_n y) \qquad (12.111) $$

$$ \alpha_n = \frac{(2n - 1)\pi}{2H} \qquad (12.112) $$

$$ A_n = \frac{(-1)^{n-1} 4\Phi_2}{(2n-1)\pi \sinh(\alpha_n L)} \qquad (12.113) $$

### Script

C'est à très peu de choses près le même que le précédent.
<figure><img src="image_73.png" /><figcaption>FIGURE 12.17 – carte des isopièzes</figcaption></figure>
<figure><img src="image_74.png" /><figcaption>FIGURE 12.18 – Solution du problème 2a</figcaption></figure>

<figure><img src="image_75.png" /><figcaption>FIGURE 12.19 – Solution du problème 2b</figcaption></figure>
Quatrième partie

Annexes
# Précision et incertitudes de mesures

## 13.1 Objectif

Il s'agit ici de différencier la précision des mesures, qui correspond aux limites des instruments, de l'incertitude des mesures liées à des fluctuations naturelles ou induites par des dépendances non explorées comme l'influence de l'humidité et de la température du laboratoire sur les remontées capillaires dans les tubes par exemple. Pour une explication détaillée je vous conseille le remarquable ouvrage de John Taylor [16]

## 13.2 Précision des mesures

### 13.2.1 Introduction

la précision d'une mesure $\delta$ est déterminée par l'instrument. Pour vous, durant ces TP d'hydrogéologie, l'instrument le plus utilisé sera la règle de mesure graduée dont la plus petite graduation est le millimètre. Une estimation optimiste de la précision correspond à 1/2 graduation soit 0.5mm. Une estimation plus simple et naturelle correspond à dire que la précision de mesure correspond à la limite basse de l'instrument soit une graduation ou 1mm.

On notera que cette valeur est toujours positive

$$ \delta h > 0 $$

En l'absence d'autres informations, la précision sert d'incertitude c'est à dire que, quand vous mesurez une hauteur d'eau piézométrique $h$, vous obtenez une valeur $h_0$ connue à la précision $\delta_h$ près. La valeur «vrai» de $h$ est donc

$$ h = h_0 \pm \delta_h \qquad (13.1) $$

### 13.2.2 Sommes de variables

Les imprécisions se propagent c'est à dire que la somme (ou la différence) de deux mesures aura une précision égale à la somme des précisions individuelles

$$ \delta(h_1 + h_2) = \delta_{h_1} + \delta_{h_2} \qquad (13.2) $$

### 13.2.3 Produit par des constantes

le produit d'une valeur par un nombre constant se retranscrit dans la précision de façon simple

$$ h = Ah_0 \pm |A\delta_h| \qquad (13.3) $$

L'apparition de la valeur absolue permet de généraliser le résultat à des valeurs négatives de $A$. On peut alors combiner nos deux premier résultats. Ainsi pour une somme de variables

$$ C = \sum_i \alpha_i X_i \qquad (13.4) $$

la précision $\delta C$ sera

$$ \delta_C = \sum_i |\alpha_i| \delta X_i = \sum_i |\alpha_i \delta X_i|, \qquad (13.5) $$
avec, rappelons le, $\delta_{X_i} \rangle 0, \forall i$.

Avec ces calculs on voit que La précision d'une moyenne de plusieurs valeurs de $h_0$ ne change pas, elle est toujours de $\delta_h$.

$$ \langle h \rangle = \frac{1}{N} \sum_{i} (h_i \pm \delta_h), \qquad (13.6) $$

soit après simplifications

$$ \langle h \rangle = \frac{1}{N} \sum_{i} h_i \pm \delta_h. \qquad (13.7) $$

Attention on verra plus loin que cette précision diffère de l'incertitude associée au calcul de $\langle h \rangle$. Rappelons que c'est une précision.

### 13.2.4 Produits de variables

Le cas de produits de variables est un peu plus délicat à traiter car les précisions ne se somment plus simplement. considérons la loi de Darcy sous sa forme

$$ Q = A \frac{H_1 - H_2}{L}. \qquad (13.8) $$

Dans cette équation le débit est le produit de trois variables : A la section de l'écoulement, $H_1 - H_2$ la différence de hauteur d'eau entre deux piézomètre, enfin L la distance qui les séparent.

Chaque grandeur est connue avec une certaine précision $\delta A$, $\delta H$, $\delta L$. Pour calculer la précision de la mesure de Q on transforme ce produit en somme en passant au logarithme

$$ \log Q = \log A + \log(H_1 - H_2) - \log L, \qquad (13.9) $$

et on dérive ce qui donne

$$ \frac{dQ}{Q} = \frac{dA}{A} + \frac{d(H_1 - H_2)}{H_1 - H_2} - \frac{dL}{L}. \qquad (13.10) $$

En remplaçant les différentielles par des précisions et en prenant la somme des valeurs absolues on obtient

$$ \delta Q = Q \left( \left| \frac{\delta A}{A} \right| + 2 \left| \frac{\delta H}{H_1 - H_2} \right| + \left| \frac{\delta L}{L} \right| \right). \qquad (13.11) $$

Question : d'où vient le 2 devant $\delta H$ ?

### 13.2.5 Fonctions de variables

Dernier cas, la précision d'un résultat obtenu par une fonction d'une mesure. Considérons la variable Y dépendant de h par l'intermédiaire d'une fonction F continue et dérivable.

$$ Y = F(h), \qquad (13.12) $$

où h est la grandeur mesurée

$$ h = h_0 \pm \delta h. \qquad (13.13) $$

Le calcul de la précision s'obtient par le développement limité de F au voisinage de $h_0$

$$ F(h_0 + \delta h) = F(h_0) + \frac{dF}{dh}(h_0)\delta h. \qquad (13.14) $$

La précision de Y est alors donnée par le second membre, toujours pris positivement,

$$ \delta Y = \left| \frac{dF}{dh}(h_0) \delta h \right| \qquad (13.15) $$

L'extension à des fonctions de plusieurs variables est laissée à votre sagacité.
## 13.3 Variables aléatoires et Incertitude

### 13.3.1 Écart type et incertitude

Si vous avez la possibilité de reproduire les mesures plusieurs (de nombreuses) fois alors la grandeur que vous mesurez peut être considérée comme une variable aléatoire. Ainsi si on répète les mesures de hauteur N fois $h_1, h_2, ..., h_i, ..., h_N$, on pourra calculer une moyenne

$$ \langle h \rangle = \frac{1}{N} \sum_{N} h_i \qquad (13.16) $$

et un écart-type, ou pluôt voire deux,
Le classique :

$$ \sigma_h = \sqrt{\frac{1}{N} \sum_{N} (h - \langle h \rangle)^2}, \qquad (13.17) $$

et le corrigé :

$$ \sigma_h = \sqrt{\frac{1}{N-1} \sum_{N} (h - \langle h \rangle)^2}. \qquad (13.18) $$

Si N est grand les deux sont très proches. Ils ne divergent que quand N est petit auquel cas il faut préférer la version corrigée.

Si la distribution est normale (ce qu'elle est le plus souvent au delà d'une trentaine de mesures) on peut écrire

$$ h_{\text{vraie}} = \langle h \rangle \pm \sigma_h \qquad (13.19) $$

qui signifie qu'avec 68% de certitude la valeur vraie de $h$, celle qu'on cherche, se trouve dans un intervalle de plus ou moins $\sigma_h$ de la moyenne calculée.

si on veut plus de certitude on prendra $2\sigma$ et on aura une certitude à 95% près. Ici donc, plutôt que de prendre la précision de la mesure on préférera prendre l'écart type de la distribution soit

$$ \delta h = \sigma_h \qquad (13.20) $$

### 13.3.2 Somme de VA et produit par une constante

La fonction de densité de probabilité de la loi normale à la forme suivante

$$ G_{X,\sigma}(x) = \frac{1}{\sigma\sqrt{2\pi}} e^{\frac{(x-X)^2}{2\sigma^2}} \qquad (13.21) $$

Il est alors intéressant de retrouver ce qu'on a vu précédemment pour la précision. Ainsi la variable $Ax$ où $A$ est une constante a pour moyenne et écart type

$$ \mu(Ax) = AX, \quad \sigma(Ax) = A\sigma_x. \qquad (13.22) $$

La somme de deux variables $x$ et $y$ a pour moyenne et écart type

$$ \mu(x + y) = X + Y, \quad \sigma(x + y) = \sqrt{\sigma_x^2 + \sigma_y^2} \qquad (13.23) $$

### 13.3.3 Fonction de VA

De façon générale une fonction $q(x, y)$ de deux variables indépendantes et normales $x$ et $y$ représente une nouvelle variable aléatoire normale dont on trouvera par développement limité autour du point $(X, Y)$ la moyenne et l'écart type. Le développement limité permet en effet d'approximer la fonction $q(x, y)$ par une somme de trois fonctions pour lesquelles il est, relativement, simple de trouver moyenne et écart type.

$$ q(x, y) = q(X, Y) + \frac{\partial q}{\partial X}(X, Y)(x - X) + \frac{\partial q}{\partial Y}(X, Y)(y - Y). \qquad (13.24) $$

Si on se concentre sur le second membre de l'égalité, $q(X, Y)$ est une constante et représente la moyenne de la distribution car les deux autres termes sont centrés donc de moyenne nulle. $\partial_x q(X, Y)(x - X)$ est une distribution centrée multiplié par une constante, elle a donc un écart type égal à $\partial_x q(X, Y)\sigma_x$. Idem pour le troisième et dernier terme.
En appliquant les principes de sommation on obtient donc

$$ \langle q(x, y) \rangle = q(X, Y), \quad \sigma_q = \sqrt{\left(\frac{\partial q}{\partial x}(X, Y)\sigma_x\right)^2 + \left(\frac{\partial q}{\partial y}(X, Y)\sigma_y\right)^2} \quad (13.25) $$

Prenons l'exemple du produit $q(x, y) = xy$. Le développement limité donne

$$ xy = XY + Y(x - X) + X(y - Y), \quad (13.26) $$

dont on déduit la moyenne et l'écart type du produit

$$ \langle xy \rangle = XY, \quad \sigma_{xy} = \sqrt{(Y\sigma_x)^2 + (X\sigma_y)^2} \quad (13.27) $$

en divisant les deux termes de l'écart type par XY on obtient

$$ \frac{\sigma_{xy}}{XY} = \sqrt{\left(\frac{\sigma_x}{X}\right)^2 + \left(\frac{\sigma_y}{Y}\right)^2} \quad (13.28) $$

## 13.4 Erreur systématique

C'est la plus complexe car elle ne suit aucune loi et en règle générale est très difficile à identifier quand on ne sait pas la forme exacte de ce qu'on cherche. Elle peut venir d'un défaut des instruments de mesure, d'un biais systématique dans la mesure (erreur de parallaxe dans les mesures optiques) ou d'une influence non aléatoire d'une tierce variable qui modifie ou décale, la valeur d'une constante.
# TP 1 : Porosité, capillarité, piézométrie

**Note**

L'objectif de cette première séance est de nous intéresser à la porosité, la teneur en eau, la saturation ainsi qu'aux effet de la capillarité. Il correspond à la lecture des chapitres 1 et 2 du cours.

## 14.1 Dispositif

On dispose du matériel suivant :

- deux récipients composés d'un bécher et d'un tube connecté à la base du bécher par une raccord coudé (piézomètre),

- d'un statif, d'une pince et d'une noix de serrage,

- d'un cristallisoir, de seringues, d'entonnoirs

- de tubes et pipettes de différents diamètres,

- d'une balance,

- de billes de verre ($\rho = 2500 \text{ kg/m}^3$, $d = 1 \text{ mm}$, 3 mm ou 4 mm),

- d'eau déminéralisée colorée.

## 14.2 Capillarité

- Mesurez le diamètre des réservoirs et des tubes

- Remplissez le petit cristallisoir d'eau colorée. Plonger un des tubes dans le cristallisoir. Mesurez la hauteur d'eau dans le récipient et dans le tube.

**Question 1**

La forme de la surface de l'eau dans le tube

1. présente une concavité tournée vers le bas;

2. présente une concavité tournée vers le haut;

3. ne présente pas de concavité.

- Recommençez l'expérience avec les tubes de diamètre différent

**Question 2**

Représentez graphiquement la hauteur de l'eau dans le tube en fonction du rayon de celui-ci. Représentez sous la forme d'une courbe la hauteur prédite par la loi de Jurin. Vos mesures

1. sont en accord avec la loi de Jurin;

2. sont en désaccord avec la loi de Jurin;

3. Nous ne pouvons pas conclure.
14.3 Porosité

— Pesez les récipients à vide !

— En utilisant l'entonnoir, remplissez un des deux récipients de 1000 ml de sables.

— Pesez la masse de sédiments.

— Tapotez le récipient (doucement) avec une barre.

Question 3

La baisse du niveau dans le récipient est due

1. à une perte de masse;

2. à un réarrangement des grains.

— Remplissez le récipient en répétant la procédure (remplir-tapoter). À la fin le volume doit rester à 1000ml, même quand vous tapotez.

Question 4

1. Quel est le volume occupé par les sédiments ?

2. Quel est le le volume des vides ?

3. Quelle valeur de porosité pouvez-vous déduire de ces mesures ?

— Ajoutez 200 ml d'eau dans le récipient vide et dans le récipient rempli de sables. Pour le récipient rempli de sédiment on injectera l'eau dans le piézomètre.

— Comparez la hauteur dans le piézomètre avec la hauteur dans le récipient sans sédiment une fois que le niveau d'eau dans le piézomètre reste statique.

Question 5

Vous pouvez maintenant effectuer une nouvelle estimation de la porosité du milieu. Quelle est sa valeur ?

— Comparez la hauteur d'eau dans les piézomètres et dans les sédiments (on estimera celle-ci par le front moyen d'humidité).

Question 6

Quelle est la taille des pores que l'on peut déduire des observations ?

— Ajoutez progressivement l'eau, par volume de 50ml, jusqu'à saturation complète du milieu.

— À chaque ajout comparez la hauteur d'eau dans le récipient vide ($h_v$) et dans le récipient rempli de sédiments ($h_s$).

Question 7

Représentez les données sur un graphique en portant les valeurs de $h_v$ en abscisse et de $h_s$ en ordonnée ?

Une relation semble se dégager entre ces deux grandeurs. Celle-ci est

1. lineaire,

2. logarithmique,

3. quadratique.
Question 8

La porosité moyenne de la colonne peut-être obtenue à partir de ce graphique car elle correspond

1. à la pente de la régression linéaire,

2. à l'inverse de la pente de la régression linéaire,

3. à la courbure de la régression polynomiale d'ordre 2
# TP 2 : Expérience de Darcy, Nappe Captive

<figure><img src="image_76.png" /><figcaption>FIGURE 15.1 – Nappe captive expérimentale</figcaption></figure>

**Note**

L'objectif de cette deuxième séance est de nous intéresser à des écoulements dans des nappes captives. Cette séance correspond à la lecture des chapitres 3, 6, 7 et 9. Le chapitre 10 permettra aux plus courageux d'aller plus loin.

## 15.1 Dispositif

le dispositif représente un cylindre de plexiglas dans lequel se trouve un aquifère composé de billes de verre (figure 15). Aux extrémités de l'aquifère, deux bouchons équipés de raccords permettent l'arrivée et la sortie de fluides. une grille empêche le passage des billes dans le circuit d'eau. L'alimentation se fait par l'intermédiaire d'une pompe 12V reliée à une alimentation. L'eau déminéralisée est pompée dans un réservoir, alimente un château d'eau qui lui-même alimente l'aquifère. En sortie l'eau passe dans un second château puis retourne au réservoir de pompage. Des tubes de plexiglas de 4 mm, les piézomètres, plongent dans l'aquifère

- Actionnez la pompe dès le début de la séance.

- Faîtes un schéma du dispositif expérimental.

- Mesurez les dimensions de l'aquifère.

- Mesurez la distance entre les piézomètres
Question 1

Les réservoirs d'alimentation sont construits afin de permettre de réguler

1. le débit des réservoirs,

2. la hauteur d'eau dans les réservoirs,

15.2 Fonctionnement d'une nappe captive

Faîtes bien attention de purger toutes les bulles d'air des piézomètres avant les premières mesures. Faîtes aussi attention que l'eau ne déborde pas dans les tubes situés sur le couvercle de la cuve. Enfin surveillez la vidange du débordement de la cuve d'alimentation.

Question 2

Pour que la nappe soit captive les deux réservoirs doivent

1. être à la même hauteur,

2. avoir une surface libre située au dessus du tube,

3. avoir une surface libre située en dessous du tube,

4. avoir leur base située au dessus du tube,

5. avoir leur base située en dessous du tube.

— Une fois la hauteur $h_2$ du réservoir de sortie fixée on ne la bouge plus.

— Fixez à une hauteur $h_1$ le réservoir d'entrée.

— Observez la montée d'eau dans les piézomètres.

— Suivez le débit.

Question 3

Le débit

1. augmente,

2. diminue,

3. reste constant.

— Quand le niveau est stable, Mesurez la hauteur dans les piézomètres.

— En parallèle mesurez le débit de façon régulière durant la mesure piézométrique. Cette mesure permettra de calculer le débit moyen $Q$

Question 4

La hauteur d'eau dans les piézomètres

1. est partout la même,

2. diminue de l'amont vers l'aval,

3. diminue de l'aval vers l'amont.

— Représentez la courbe $h(x)$. On choisira judicieusement l'origine de l'axe des x.

Question 5

La courbe qui relie les hauteurs mesurées dans chacun des piézomètres à la forme

1. d'une droite,

2. d'une parabole,

3. d'un polynome d'ordre trois.
Question 6

Quelle est la valeur de la pente S de cette courbe ? Cette pente correspond

1. à la vitesse de l'écoulement,

2. à la pression dans le tube,

3. au gradient de pression dans le tube.

— Répétez les mesures pour trois hauteurs $h_1$ différentes.

— À chaque changement de hauteur observez le changement dans les piézomètres.

Question 7

Représentez graphiquement pour chaque hauteur d'entrée le débit mesuré en fonction de la pente S reliant les hauteurs piézométriques. La relation qui apparaît entre Q et S semble

1. linéaire,

2. logarithmique,

3. quadratique.

— Calculez la surface A traversée par l'écoulement.

Question 8

Le rapport Q/A correspond

1. à la vitesse réelle de l'eau,

2. à la vitesse moyenne de l'eau.

3. à la vitesse de Darcy,

Question 9

À partir de vos mesures et de la loi de Darcy calculez la conductivité hydraulique de votre aquifère (en m/s) :

— Fermer le robinet d'alimentation et arrêter l'alimentation électrique et observez l'évolution de la hauteur d'eau dans les piézomètres.

Question 10

La hauteur d'eau dans la nappe tend progressivement vers

1. 0,

2. la hauteur de la surface libre du réservoir d'entrée,

3. la hauteur de la surface libre du réservoir de sortie.

15.3 Pompage dans une nappe

— Créez une nappe statique dont le niveau sera la plus haut possible. Les deux réservoirs devront être alimentés.

— Ouvrez le robinet central et récupérer l'eau dans un bécher

— Mesurez toutes les deux minutes le volume d'eau accumulée pendant un intervalle de temps que vous choisiserez.

— En parallèle mesurez les hauteurs piézométriques et le débit au débitmètre
Question 11

Le débit de pompage correspond à

1. la valeur,

2. la moitié de la valeur,

3. le quart de la valeur

du débit mesuré au débitmètre

Question 12

Représentez les données $h$ de chaque piézomètre en fonction de $r$ la distance d'un piézomètre au puits.

La relation entre $h$ et $r$ semble

1. linéaire,

2. logarithmique,

3. quadratique.

Question 13

Ce résultat est-il en accord avec la forme $h(r)$ proposée au chapitre 9 pour un pompage dans une nappe captive ?

1. Oui car le domaine de définition et les symétries du dispositif sont identiques,

2. Oui car les conditions aux limites sont semblables,

3. Non car le domaine de définition et les symétries du dispositif ne sont pas identiques.

4. Non car les conditions aux limites diffèrent.
# TP 3 : Nappe libre, transfert de polluant

<figure><img src="image_77.png" /><figcaption>FIGURE 16.1 – Nappe libre expérimentale</figcaption></figure>

**Note**

L'objectif de cette dernière séance est de comprendre le fonctionnement d'une nappe libre. Afin de préparer la séance vous devrez lire des chapitres 4 et 5. Les chapitres 8 et 11 permettront aux plus courageux d'aller plus loin

## 16.1 Dispositif

le dispositif représente un aquifère composé de billes de verres limité par des plaques de plexiglas sur les côtés et sur le bas (figure 16). Aux extrémités de l'aquifère deux réservoirs sont limités par des grilles qui permettent le passage de l'eau mais empêchent le mouvement des grains composant l'aquifère. L'alimentation en eau se fait par l'intermédiaire d'une pompe 12V reliée à une alimentation. L'eau déminéralisée est pompée dans un réservoir, alimente une premier château d'eau qui lui-même alimente le réservoir amont de l'aquifère. En aval de l'aquifère l'eau est récupérée dans le second réservoir puis passe dans le second château avant de retourner au réservoir de pompage. des piézomètres sont connectés à la base de l'aquifère.

— Veillez à actionner la pompe dès le début de la séance.

— Faîtes un schéma du dispositif expérimental.

— Mesurez les dimensions de l'aquifère.

— Mesurez la distance entre les piézomètres

## 16.2 Etude d'une nappe libre

On notera *L* la longueur totale de la nappe, *w* la largeur de la nappe, *x* la distance horizontale entre un piézomètre et l'entrée d'eau.

— Positionnez les réservoirs d'entrée et de sortie afin de vous assurer que la nappe sera libre.

— Observez la montée de la nappe.

— Mesurez la hauteur dans les piézomètres et le débit.
Question 1

La surface A traversée par l'écoulement au niveau de chaque piézomètre est

1. $h * w$,

2. $h * x$,

3. $h * (L - x)$,

4. $h^2$.

— Calculez la vitesse de Darcy $U_D$ (cf TP 2) au niveau de chacun des piézomètres ?

— Représentez graphiquement $U_D$ en fonction de la distance à l'entrée.

Question 2

La vitesse de Darcy

1. augmente avec la distance à l'entrée,

2. diminue avec la distance à l'entrée,

3. reste constante quelque soit le piézomètre.

Question 3

La forme de la surface libre est

1. linéaire,

2. logarithmique,

3. polynomiale d'ordre 2.

— Faîtes trois fois l'opération en modifiant les conditions aux limites

Question 4

En utilisant les calculs de conductivité hydraulique du TP2 et votre cours représentez sur le même graphique les mesures piézométriques et la forme analytique de la surface piézométrique. La surface théorique coïncide-t-elle avec vos mesures ?

1. Oui,

2. Non.

Question 5

Pour arriver à cette forme analytique on doit faire l'hypothèse que

1. l'écoulement est horizontal,

2. l'écoulement est vertical,

3. l'écoulement suit la surface piézométrique.

Cette hypothèse est-elle toujours valable dans votre expérience ? Si non où et pourquoi ne l'est-elle pas ?

16.3 Propagation d'un polluant

— À l'aide d'un seringue injectez du colorant dans la nappe en amont de l'expérience.

— Suivez l'évolution de la tache de colorant toutes les 30s.

— Estimez la taille de la tache et calculer sa vitesse moyenne de déplacement.

— Tracez l'évolution de la vitesse de la tache en fonction du temps.

— Tracez l'évolution de la taille de la tache en fonction du temps.
Question 6

La vitesse moyenne de la tache est

1. inférieure à la vitesse de Darcy,

2. égale à la vitesse de Darcy,

3. supérieure à la vitesse de Darcy.

Question 7

À partir de l'observation de l'évolution de la grosseur de la tâche on peut dire que le mécanisme principal
qui contrôle le déplacement du colorant est

1. l'advection du colorant,

2. la diffusion moléculaire du colorant,

3. la descente de la tache sous l'effet de la gravité.

Question 8

Le mécanisme principal peut se déduire de la valeur du nombre de Peclet (Pe = UL/D ou U est la vitesse réelle de l'écoulement, L la longueur caractéristique de celui-ci et D a diffusivité du colorant. chapitre 5). En considérant une diffusivité moléculaire de l'ordre de D = 10⁻⁹ m²/s on peut dire que

1. Pe << 1

2. Pe ≃ 1

3. Pe >> 1

Question 9

Ce résultat est il cohérent avec votre hypothèse (question 7) ?

1. Oui.

2. Non.
# Changelog des versions sur Hal

**versions 1 - 5** Première forme en deux parties I Notes de cours et II applications, correction d'erreurs et réécriture de paragraphes

**versions 6-7** Changement de plan en 4 parties 2 de cours + 2 d'applications, solutions de l'équation de Laplace et solutions des exercices

**version 8** Post covid : les cours deviennent principalement des TP, changement de plan I hydrogéologie 1D (colle avec les TP), II écoulements dans le plan (suppléments) + applis et solution des exercices

**versions 9** Calcul d'incertitude en annexe.

**version 10** Extension de l'annexe sur les incertitudes, énoncés des TP fournis en annexe, corrections et réécriture de paragraphes.

**version 11** Corrections diverses + ajout d'exercices sur les nappes libres.

**version 12** Corrections diverses + modification des énoncés de TP (QCM); ajout concernant la capacité de rétention;
Bibliographie

[1] J. Bear. *Dynamics of Fluids in Porous Media*. Dover, New York, 1972.

[2] J. Boussinesq. Recherches théoriques sur l'écoulement des nappes d'eau infiltrées dans le sol et sur le débit des sources. *Journal de mathématiques pures et appliquées*, 10 :5-78, 1904.

[3] W. Brutsaert. *Hydrology and introduction*. Cambridge University Press, 2005.

[4] H.S. Carslaw and J.C. Jaeger. *Conduction of heat in solids*. Clarendon press, 1992.

[5] P Corbier, G Karnay, B Bourgine, and M Saltel. Gestion des eaux souterraines en région aquitaine-reconnaissance des potentialités aquifères du mio-plio-quaternaire des landes de gascogne et du médoc en relation avec les sage. Rapport final BRGM/RP-56475-FR, BRGM, 2008.

[6] H. Darcy. *Les fontaines publiques de la ville de Dijon : exposition et application des principes à suivre dans les questions de distribution d'eau*. Victor Dalmont, Paris, 1856.

[7] S. L. Dingman. *Physical hydrology*. Prentice Hall, second edition, 2002.

[8] S. J. Farlow. *Partial differential equations for scientists and engineers*. Dover, 1993.

[9] R.A. Freeze and J.A. Cherry. *Groundwater*. Prentice-Hall, Englewood Cliffs, NJ, 1979.

[10] E. Gilli, C. Mangan, and J. Mudry. *Hydrogéologie-4e éd. : Objets, méthodes, applications*. Dunod, 2016.

[11] A.I. Johnson et al. Specific yield : compilation of specific yields for various materials. Technical Report 1662-D, 1967.

[12] G. Marsily (de). *Hydrogéologie quantitative*. Masson, Paris, 1981.

[13] Owen M Phillips. *Flow and reactions in permeable rocks*. Cambridge University Press, 1991.

[14] P. Polubarinova-Kochina. *Theory of ground water movement*. Princeton University Press, 1962.

[15] F. Press and R. Siever. *Understanding Earth*. W.H. Freeman and Company, New York, 1998. 2nd edition.

[16] John Taylor. *Introduction to error analysis, the study of uncertainties in physical measurements*. University Science Books, 1997.

[17] F Verley, F Brunson, P Verjus, and M Cholez. Nappe de beauce-piézométrie hautes eaux 2002. Technical report, Tech. rep., DIREN Centre et Ile-de-France, Orléans, France, 2003.

[18] JF Vernoux, J Barbier, M Donsimoni, JJ Seguin, and J Vairon. Etude hydrogeologique du plateau de saclay (essone). *Rapport BRGM SGR/IDF R*, 40840, 1999.