# Modélisation

La modélisation de l’écoulement et du transport hydrogéologique qui sous-tend le code THYRSIS, est le couplage d’un modèle d’écoulement et transport verticaux en zone non saturée et d’un modèle d’écoulement et transport horizontaux en zone saturée, eux-mêmes fondés sur un champ d’écoulement en régime permanent, construit indépendamment, pour chaque site étudié.

Ce champ d’écoulement peut être importé dans THYRSIS à partir de calculs réalisés indépendamment, ou bien élaboré à l’aide du plugin, en créant un maillage 2D avec le mailleur GMSH [GR09] et en déterminant les différentes grandeurs caractéristiques de l’écoulement avec la procédure d’inversion HYNVERSE [RJ08], [RT16].

## Écoulement de la nappe en régime permanent

L’écoulement de la nappe en régime permanent doit être préalablement défini, car il est utilisé pour les calculs d’écoulement transitoire et de transport en zone saturée, mais aussi pour construire les simulations d’écoulement et transport verticaux en zone non saturée, notamment pour déterminer, en chaque zone d’injection :

- les épaisseurs de colonne non saturée,

- la perméabilité à saturation de chaque colonne.

L’écoulemen de la nappe en régime permanent est un calcul 2D dont les paramètres sont :

- l’infiltration, constante en régime permanent, mais pas nécessairement homogène,

- l’épaisseur de nappe, variable sur tout le domaine modélisé,

- le champ de perméabilité.

Les sorties du calcul sont :

- le champ de potentiel, ou champ piézométrique, i.e. l’altitude du toit de la nappe dans le cas d’une nappe libre.

- le champ des vitesses de Darcy.

En outre, c’est le même maillage qui est utilisé pour l’écoulement et le transport en zone saturée.

Toutes ces informations sont stockées dans des bases de données Spatialite.

L’élaboration du modèle d’écoulement de la nappe en régime permanent est une phase longue et complexe de la modélisation hydrogéologique. Le modèle est construit en exploitant un grand nombre de données concernant la géologie (extension verticale de l’aquifère, par détermination du son niveau imperméable), la topographie, la météorologie (détermination de l’infiltration) et l’hydrogéologie. Ces données sont obtenues par étude bibliographique et analyse des forages réalisés sur le site et aux environs, le domaine à modéliser s’étendant souvent bien au delà des limites du site étudié. La détermination du champ de perméabilité est particulièrement difficile, et repose le plus souvent sur des techniques d’inversion, alliant géostatistique et optimisation (code HYNVERSE [RJ08], [RT16]).

## Écoulement de la nappe en régime transitoire

Le modèle d’écoulement de la nappe en régime transitoire est fondé sur le modèle élaboré en régime permanent, en y ajoutant une valeur d’emmagasinement, égale à la porosité efficace dans le cas d’une nappe libre.

Une simulation d’écoulement en régime transitoire requiert la connaissance de la variation de l’infiltration dans le temps, qui est fournie sous forme d’un fichier externe, selon une syntaxe adaptée (cf. [Fichiers utilisateur]).

## Écoulement et transport en zone non saturée

La modélisation de la zone non saturée est fondée sur la construction de colonnes verticales 1D. Le nombre de colonnes est égal au nombre de zones d’injection, et leur épaisseur s’obtient par différence entre la topographie et la piézométrie en régime permanent.

L’écoulement est régi par l’équation de Richards, associée à la formulation de Van Genuchten [VG80], dans sa version modifiée par [IVB06]:

\[\begin{split} S_{e}=\begin{cases} \frac{1}{S_{c}}\left[1+\left(\alpha h\right)^{n}\right]^{-m} & h>h_{e}\\ 1 & h\leq h_{e} \end{cases} \end{split}\]
\[\begin{split} K=\begin{cases} K_{s}\sqrt{S_{e}}\left[\frac{1-\left(1-\left(S_{e}S_{c}\right)^{1/m}\right)^{m}}{1-\left(1-S_{c}^{1/m}\right)^{m}}\right]^{2} & S_{e}<1\\ K_{s} & S_{e}\geq1 \end{cases} \end{split}\]

où :

- \(h\) : potentiel de pression (L)

- \(h_{e}\) : pression d’entrée (L)

- \(K\) : perméabilité (\(L.T^{-1}\))

- \(K_{s}\) : perméabilité à saturation (\(L.T^{-1}\))

- \(S_{e}=\frac{\theta-\theta_{r}}{\theta_{s}-\theta_{r}}\) : teneur en eau réduite, ou saturation effective

- \(S_{c}=\left[1+\left(\alpha h_{e}\right)^{n}\right]^{-m}\)

- \(\theta\) : teneur en eau

- \(\theta_{r}\) : teneur en eau résiduelle

- \(\theta_{s}\) : teneur en eau à saturation

- \(m=1-\frac{1}{n}\)

- \(n\) and \(\alpha\) (\(L^{-1}\)) sont des paramètres d’ajustement du modèle.

Notons qu’il y a correspondance biunivoque entre la teneur en eau \(\theta\) et la saturation (ou degré de saturation) \(s\) via la porosité \(\omega\) : \(\theta=s\omega\) .

Lorsque \(h_{e}=0\) , on a \(S_{c}=1\) et on retrouve la formulation de Van Genuchten originale.

### Maillage des colonnes ZNS

Les colonnes sont décomposées en mailles d’épaisseur centimétrique, afin d’une part de garantir la stabilité numérique au niveau de l’injection, même avec de faibles valeurs de dispersivité, et d’autre part de garantir des bilans de masse exacts.

Plus précisément, l’épaisseur de la première maille, en haut de la colonne, est fixée à :

\[ min\left(1cm,\frac{e_{ZNS}}{100}\right) \]

où \(e_{ZNS}\) désigne l’épaisseur de la colonne insaturée, exprimée en cm.

Cette épaisseur de la première maille est donc le plus souvent égale à 1 cm, sauf pour des épaisseurs de colonne inférieures à 1 m. Ensuite, afin d’accélérer les calculs, les épaisseurs de maille augmentent en progression géométrique jusqu’au milieu de la colonne, puis diminuent symétriquement pour atteindre en bas de colonne une valeur égale à celle du haut de colonne. La progression géométrique est déduite du rapport entre l’épaisseur de la maille au milieu de la colonne et celle de la première maille, rapport qui est paramétrable par l’utilisateur, et égal par défaut à 5 (coef_maillage dans la table [simulations]).

De plus, l’épaisseur d’une maille est contrainte de rester inférieure à \(0.1/\alpha\) (coefficient de van Genuchten) de manière à garantir la convergence dans le cas de valeurs élevées de \(\alpha\) value.

### Paramètres ZNS

A la différence de la zone saturée, l’écoulement et le transport en zone non saturée sont traités en une même simulation, qui débute par un calcul d’écoulement en régime permanent, suivi d’un calcul de transport (et d’écoulement le cas échéant) en régime transitoire

Les paramètres du modèle de transport en zone saturée sont :

- la porosité cinématique : \(\omega_{c}\) (de fracture et/ou de matrice, selon le site étudié),

- la porosité totale \(\omega_{t}\),

- la masse volumique \(\rho_{s}\) du composé,

- les coefficients \(\alpha\) et n de van Genuchten,

- les saturations résiduelle et maximale,

- la dispersivité longitudinale (la dispersivité transverse n’existe pas en 1D),

- DK, Coefficient de partage (m³/kg),

- la limite de solubilité du composé \(L_{s}\) (cf. [Injection en masse]).

Parmi ces paramètres, la porosité totale et la masse volumique sont utilisées pour calculer le coefficient de retard R intervenant dans l’équation de transport, selon la formule :

\[ R=1+\frac{\left(1-\omega_{t}\right)}{\omega_{c}}\rho_{s}K_{d} \]

De même qu’en zone saturée, la simulation en zone non saturée d’une infiltration transitoire requiert la connaissance de la variation de l’infiltration dans le temps, qui est fournie sous forme d’un fichier externe, selon une syntaxe adaptée [Fichier utilisateur]).

### Dépôt réparti sur tout le domaine

Il est possible de prendre en compte un dépôt réparti sur tout le domaine modélisé, soit en simulant autant de colonnes verticales insaturées qu’il y a de noeuds du maillage, ce que permet la parallélisation du plugin, soit en simulant des groupes de colonnes insaturées d’épaisseur et de perméabilité voisines, de manière à réduire le nombre de simulations. Dans ce cas, deux paramètres doivent être renseignés (table [simulations]) :

- pas_echant_zns désigne le pas d’échantillonnage des colonnes. Il définit la taille des colonnes d’un même groupe. Ainsi, un pas de 5 m définit des groupes de colonnes de taille [0,5m], [5m,10m], [10m,15m], etc. Le nombre de groupes est d’autant plus réduit que le pas est élevé.

- rapport_max_permeabilite désigne le rapport maximal de perméabilité qui peut exister entre les perméabilités (à saturation) de deux colonnes d’un même groupe.

## Transport en zone saturée

Le transport en zone saturée est fondé sur le modèle d’écoulement en régime [permanent] ou [transitoire]

La condition d’injection dans la nappe est :

- soit un flux, une concentration ou une masse, définis explicitement, via l’interface graphique, si la migration en zone non saturée ne doit pas être modélisée (cf. [injection]),

- soit un flux, défini implicitement à partir des flux de composé calculés au pied des colonnes non saturées. Dans ce cas, le flux est injecté aux cellules (OpenFOAM) ou aux noeuds (METIS) du maillage correspondant aux différentes zones d’injection, de manière à couvrir la surface d’injection.

Les paramètres du modèle de transport en zone saturée sont :

- la porosité (de fracture et/ou de matrice, selon le site étudié),

- les dispersivités longitudinale et transverse,

- le coefficient de rétention du composé.

Dans le cas d’un régime transitoire d’écoulement, le coefficient d’emmagasinement est égal à la porosité (efficace).

## Caractéristiques de l’injection

L’injection d’un composé est définie par les caractéristiques suivantes :

- coordonnées (x, y) du centre de la zone d’injection,

- profondeur,

- surface S concernée par l’injection,

- flux F, exprimé en kg/s ou Bq/s,

- durée de l’injection \(\Delta T\).

Par commodité, l’injection peut également être définie en terme de concentration ou de masse, comme indiqué ci-après. Une fuite, caractérisée par un volume d’eau, peut aussi être modélisée.

### Injection en concentration

Dans le cas d’une injection définie par une concentration `C`, le flux `F` s’obtient par :

\[ F=ISC \]

où `I` désigne le flux d’infiltration d’eau de pluie, défini avec le modèle d’écoulement de la nappe en régime [permanent].

### Injection en masse

“Dans le cas d’une injection définie par une masse `M`, le flux injecté `F` et la durée d’injection \(\Delta T\) sont calculés à partir de la connaissance de la limite de solubilité \(L_{s}\) du composé injecté :

\[ F=ISL_{s} \]
\[ \Delta T=\frac{M}{F} \]

### Fuite

Il est possible d’accompagner l’injection de matière par une injection d’eau en précisant le volume d’eau \(V\) injecté, le flux \(F\) ou la concentration \(C\). Dans ce cas, la vitesse d’injection \(v\) est calculée par :

\[ v=\frac{V}{S \Delta T} \]

Cette vitesse ne peut être supérieure à la perméabilité à saturation \(K_{s}\) Si tel est le cas, ou plus exactement si :

\[ v > 0.9 K_{s} \]

alors la durée d’injection est modifiée par :

\[ \Delta T=\frac{V}{0.9 S K_{s}} \]

et un message est émis signalant qu’il est conseillé d’augmenter la surface d’injection en prenant :

\[ S=\frac{V}{0.9 K_{s} \Delta T} \]

Si la concentration \(C\) est fournie, le flux s’en déduit par :

\[ F=vSC \]

Le tableau ci-dessous synthétise les différents calculs de flux, selon les données disponibles :

Données | \(V=0\) | \(V>0\)
\(F, \Delta T\) | \(F\) | \(v=min(\frac{V}{S \Delta T}, 0.9 K_{s})\), \(\Delta T = max(\Delta T, \frac{V}{0.9 S K_{s}})\)
\(C, \Delta T\) | \(F=ISC\) | \(v=min(\frac{V}{S \Delta T}, 0.9 K_{s})\), \(\Delta T = max(\Delta T, \frac{V}{0.9 S K_{s}})\), \(F=vSC\)
\(M, L_{s}\) | \(F=ISL_{s}\) | Pas de calcul

## Analyses d’incertitudes

THYRSIS offre la possibilité de faire varier l’ensemble des paramètres de l’écoulement et du transport en zones non saturée et saturée, selon des lois de probabilités standard (uniforme, log-uniforme, normale, log-normale), de manière à estimer des intervalles de confiance sur les concentrations aux points de calcul, ou des cartes de probabilité de dépassement de seuil. Il suffit pour cela de préciser un nombre de simulations à réaliser (supérieur à 1) et d’indiquer les lois de probabilité à utiliser pour chaque paramètre, ainsi que leurs valeurs minimale et maximale. Un jeu de paramètres est alors généré par [Hypercube Latin], et THYRSIS réalise les simulations demandées (réparties sur le nombre de processeurs spécifié dans les [Préférences]), et agrège les résultats (moyenne, écart-type).

Les représentations graphiques tiennent compte de ces simulations en traçant les intervalles de confiance ou les cartes de probabilité.

## Double milieu

Le plugin THYRSIS implémente, via le code METIS, le modèle d’écoulement double milieu de Gerke & Van Genuchten [GVG93], Celui-ci est fondé sur l’hypothèse que l’équation de Richards s’applique à la fois au milieu matriciel et au milieu fracturé. Pour un système vertical mono-dimensionnel, les équations d’écoulement dans un volume élémentaire représentatif s’écrivent, pour chacun des milieux, fracture (f) et matrice (m) respectivement :

Dans la porosité de fracture :

\[ C_{f}\frac{\partial h_{f}}{\partial t}=\frac{\partial}{\partial z}\left(K_{f}\frac{\partial h_{f}}{\partial z}-K_{f}\right)-\frac{\Gamma_{w}}{w_{f}}-S_{f} \]

Dans la porosité de matrice :

\[ C_{m}\frac{\partial h_{m}}{\partial t}=\frac{\partial}{\partial z}\left(K_{m}\frac{\partial h_{m}}{\partial z}-K_{m}\right)-\frac{\Gamma_{w}}{1-w_{f}}-S_{m} \]

avec :

\[ C_{f,m}=S_{w}S_{s}^{f,m}+\omega_{f,m}\frac{\partial S_{w}}{\partial h} \]
\[ \Gamma_{w}=\alpha_{w}\left(h_{f}-h_{m}\right) \]
\[ \alpha_{w}=\alpha_{w}^{*}K_{a} \]
\[ \alpha_{w}^{*}=\frac{\beta}{a^{2}}\gamma_{w} \]

et :

\(K_{f}\) and \(K_{m}\) (\(L.T^{-1}\)) : perméabilités de fracture et de matrice,

\(h_{f}\) and \(h_{m}\) (L) : potentiel hydraulique dans la fracture et la matrice,

\(\Gamma_{w}\) (\(T^{-1}\)) : terme d’échange décrivant le transfert d’eau entre la fracture et la matrice,

\(S_{f}\) and \(S_{m}\) (\(T^{-1}\)) : termes source dans la fracture et la matrice,

\(w_{f}\) : fraction volumétrique de fracture dans le système total,

\(C_{f}\) and \(C_{m}\) (\(L^{-1}\)) : capacités spécifiques de l’eau dans la fracture et la matrice,

\(S_{w}\) : degré de saturation en eau,

\(S_{s}\) (\(L^{-1}\)) : coefficient d’emmagasinement spécifique dans la fracture et la matrice,

\(\omega_{f}\) and \(\omega_{m}\) : porosité de la fracture et de la matrice,

\(K_{a}\) (\(L.T^{-1}\)) : perméabilité de transition entre la fracture et la matrice,

\(a\) (L) : demi-distance moyenne entre les fracture,

\(\beta\) : facteur de structure représentant la géométrie du milieu fracturé (\beta=3 dans le cas d’un ensemble de fractures orthogonales régulièrement espacées),

\(\gamma_{w}\) : facteur empirique égal à 0,4.