Developpement d’une sonde aquatique autonome pour la cartographie des drains karstiques noyes

Approche distribuée 

L’approche distribuée permet de tenir compte des hétérogénéités spatiales des aquifères. Elle nécessite la discrétisation spatiale des propriétés hydrodynamiques du système modélisé. Ceci revient à résoudre des modèles d’écoulements basés sur des équations aux dérivées partielles (E.D.P.) à l’aide de méthodes numériques (succinctement décrites ci-après). Il existe trois échelles d’observation des écoulements (Sahimi, 2012) : (i) à l’échelle microscopique le pore est la longueur caractéristique, les écoulements sont décrits par les équations de Stokes, (ii) à l’échelle locale la taille d’échantillon caractéristique varie depuis la taille d’une carotte jusqu’à celle du puits, les écoulements sont décrits par des équations homogénéisées, comme par exemple la loi de Darcy (1856), à l’échelle macroscopique la taille de l’échantillon est celle du réservoir et les équations qui régissent les écoulements sont pour la plupart des équations de bilan. Pour la modélisation des écoulements dans le karst, Teutsch et Sauter (1991) proposent une classification des méthodes de modélisation distribuée selon (i) leur description des hétérogénéités spatiales, (ii) leur applicabilité et (iii) l’effort d’investigation à mener. Ces méthodes de modélisation évoluent depuis des modèles simples et homogènes (milieux poreux équivalents) vers des modèles déterministes complexes (réseaux de conduits ou de fractures) en passant par des modèles de complexités intermédiaires (double continuum ou milieu hybride). Ces modèles sont décrits cidessous. 

 Milieu poreux equivalent 

Si le milieu souterrain est composé d’un réseau de conduits ou de fractures assez dense et qu’il admet une connectivité hydraulique importante, ou à l’inverse que le réseau est faiblement karstifié, un VER peut être défini. Alors, le comportement hydrodynamique du milieu peut être modélisé à l’aide d’un modèle continu, de type milieu poreux équivalent. Le flux traversant l’aquifère suit ainsi la loi de Darcy (1856). Cette loi est valable pour un fluide monophasique incompressible et s’exprime par : V = − k µ ∇ (P + ρgz), (1.3) avec V(x, t) [LT−1 ] est la vitesse de filtration (vit. moyenne qui vérifie l’équation de continuité) du fluide à travers le milieu, k [L2 ] le tenseur de perméabilité intrinsèque du milieu, ∇ est l’opérateur divergence (multiplié à un tenseur, le résultat est un vecteur), µ [ML−1T −1 ] la viscosité dynamique du fluide, P la pression [ML−1T −2 ], ρ [ML−3 ] la masse volumique du fluide, g [LT−2 ] l’accélération de la pesanteur et z [L] la cote altimétrique. En hydrogéologie, deux hypothèses sont souvent faites. On suppose que (i) le fluide (eau) est incompressible (c.-à-d. que le terme ρg/µ est invariant) et dans ce cas on définit la conductivité hydraulique K¯ (x) = ρgk/µ), (ii) que les vitesses réelles sont lentes et que l’on peut négliger le terme de charge dynamique défini dans la charge hydraulique, pour exprimer la charge piézométrique h(x, t) [L]. La charge peut alors s’exprimer comme le rapport de la pression hydrostatique par le poids par unité de volume d’eau auquel s’ajoute la côte z (c.-à-d. h = P/(ρg) + z). La côte piézométrique est alors confondue avec la charge. Sauf mention, nous adoptons dans la suite du manuscrit ces deux hypothèses. La loi de Darcy peut ainsi s’écrire : V = −K∇h. 16 Chapitre 1. Eléments d’hydrogéologie karstique Il paraît important ici d’approfondir quelques points liés à cette loi, puisque ceux-ci seront repris dans la suite du manuscrit, en particulier dans les Sect. 2 et 3. En combinant la loi de Darcy, avec l’équation de conservation de la masse et l’équation d’état du fluide, on obtient une équation de type diffusion. Cette équation permet de calculer l’évolution du flux et de la pression du fluide dans le milieu, selon une unique variable indépendante : pression ou charge piézométrique. Suivant le type d’aquifère considéré (libre ou confiné), la forme de l’équation sera différente. Dans le cas d’aquifères confinés, la pression du fluide est plus élevée que la pression atmosphérique. Le niveau de l’eau dans une telle situation est supérieur à celui du toit de l’aquifère, éventuellement il peut aussi être supérieur à celui de la surface du sol (artésianisme). L’aquifère est ainsi constamment saturé et l’évolution de la nappe d’eau dans ce type d’aquifère est décrite par l’équation suivante : Ss ∂h ∂t + q = ∇ ·  K∇h  , (1.4) avec K(x) le tenseur de conductivité hydraulique, q(x, t) [L3T −1/L 3 ] est le terme source. Ce débit volumique défini en chaque point est conventionnellement positif en cas de prélèvement dans le milieu et négatif en cas d’apport. Ss(x) [L −1 ] est le coefficient d’emmagasinement spécifique de l’aquifère. Ce coefficient correspond au volume d’eau libéré par unité de volume du milieu poreux sous l’effet d’un abaissement de charge h(x, t) unitaire. Cette quantité est définie en fonction des paramètres du milieu, la porosité totale φt , les coefficients de compressibilité du sol cs, du liquide cl et du milieu cp et en fonction de la masse volumique du fluide (ρg), elle s’écrit : Ss = ρgφt(cl − cs + cp/φt). L’Equation 1.4 exprime que la variation par unité de temps de la charge piézométrique à l’intérieur d’un volume V , multipliée par la capacité que ce volume possède d’emmagasiner du fluide (premier terme du membre de gauche) est égale au flux de la vitesse qui entoure le volume V (membre de droite). Le terme q (second membre de gauche) ne représente qu’un terme d’apport ou de prélèvement. Cette équation de bilan équilibre les flux entrants et sortants, y est ajouté un abaissement de charge par unité de temps. La résolution de cette équation nécessite de connaître les paramètres du milieu et du fluide et les conditions initiales (CI) et aux limites (CL). La résolution peut se faire de façon analytique (on cherche une solution exacte), semi-analytique (on détermine une solution en domaine non temporel, p. ex. domaine Laplace) ou numérique (solution approchée). Une fois l’équation résolue, h est décrit en tout point du domaine (sur x) et pour tout t. De plus, du fait que cette équation est linéaire, le principe de superposition peut être appliqué. Pour rappel, ce principe physique dit que si (h1, q1) et (h2, q2) sont deux solutions particulières de Eq. 1.4 qui vérifient des CL données, alors toutes combinaisons linéaires des deux solutions, sont aussi aussi solution de l’équation. Ce principe, important en physique, permet de considérer une variable qui permet une écriture simplifiée des CI du problème : le rabattement s(x, t). Il est défini comme la différence entre la solution du problème stationnaire h0(x) et la solution du problème transitoire h(x, t) qui, vérifie les même CL que h0(x) et de CI h0(x) [s(x, t) = h0(x) − h(x, t)]. Dans le cas d’aquifères libres, l’aquifère n’est pas saturé dans son intégralité. La surface piézométrique fluctue et représente la limite entre les zones non saturées et saturées du milieu. Dans ce cas, la section d’écoulement n’est pas constante. Sur la longueur les propriétés hydrodynamique du milieu (en particulier sa transmissivité) sont alors modifiées. Toutefois, on peut formuler les deux hypothèses suivantes (reprises par la suite dans la Sect. 3). Premièrement, les vitesses d’écoulements (lentes) selon la composante verticale z peuvent être négligées face aux vitesses d’écoulement horizontales (plus importantes), ces dernières sont toutes parallèles entre elles sur la même verticale (hyp. de Dupuit-Forscheimer). Deuxièmement, le tenseur de conductivité hydraulique admet la verticale comme une de ses directions principales. L’équation qui régit l’écoulement s’écrit : Sy ∂h ∂t + Q = ∇xy ·  K (h − ι) ∇xyh  , (1.5) avec Sy(x, y) (L −1 ) le coefficient d’emmagasinement spécifique (qui rend compte de la porosité efficace φe) du milieu, les indices x et y à l’opérateur ∇ indique les dérivées selon ces axes (on omet z qui est négligé), K le tenseur de conductivité hydraulique est réduit à une matrice de quatre éléments et h(x, y, t) [L] représente la charge hydraulique sur une verticale (et de fait la cote de la surface libre), ι(x, y) représente une fonction qui décrit le mur de l’aquifère, et enfin Q(x, y, t) [L3T −1/L2 ] est le débit prélevé par unité de surface de nappe. Il est supposé que les échanges entre la nappe et l’extérieur s’effectuent sur toute son épaisseur. On note que cette équation est non-linéaire, ce qui rend sa résolution complexe. Par la suite l’utilisation d’une telle formulation est discutée en Sect. 3. Pour les aquifères karstiques, cette approche a été utilisée pour modéliser des flux à l’échelle régionale d’un aquifère (p. ex. dans les travaux de Angelini et Dragoni (1997); Josnin et al. (2000); Scanlon et al. (2003)). Toutefois, (Teutsch, 1990) indique que celle-ci devrait être évitée pour des aquifères très karstifiés, notamment à cause de l’hypothèse faite quant aux régimes d’écoulements (non turbulents), ce qui, dans le cas des aquifères karstiques peut être discuté. 

 Double porosité – double perméabilité 

Les modèles double-porosité ont été initialement développés pour la simulation des transferts de flux et de pression dans des réservoirs de grès. L’approche continue considère deux régions définies en tout point du domaine, deux continua. Dans chaque continuum, l’écoulement se comporte comme dans un milieu poreux homogène qui occupe tout le réservoir. Ainsi dans cette approche les deux continua se recouvrent. En chaque point, deux valeurs moyennes de pression locale sont évaluées pour un volume élémentaire représentatif du continuum considéré. Le modèle à double-porosité et double-perméabilité de Barenblatt et al. (1960) considère deux milieux : le premier continuum est à perméabilité faible et forte porosité, le second continuum est à forte perméabilité et faible porosité. Les flux qui traversent chaque continuum peuvent être décrits par deux équations de diffusivité liées par un terme d’échange. Ce terme est traité dans chaque domaine par un terme puits/source. Les équations qui régissent l’écoulement peuvent s’écrire :    Ss1 ∂h1 ∂t = ∇ ·  K1∇h1  + q1 + Q21, Ss2 ∂h2 ∂t = ∇ ·  K2∇h2  + q2 + Q12, (1.6) où les indices 1 et 2 rapportent respectivement au premier ou second continuum, Ssi(x) [L −1 ] le coefficient d’emmagasinement spécifique du milieu i, Ki [LT −1 ] le tenseur de conductivité hydraulique du milieu i, hi(x, t) la charge [L] locale moyenne du milieu i, qi est le terme source pour le milieu i et Qij = −Qji [L3T −1L 3 ] est le flux échangé entre les deux milieux i et j par unité de volume. Le flux d’échange entre les milieux selon Barenblatt et al. (1960) est une fonction linéaire de la différence de pression moyenne entre les deux milieux, il correspond donc à un état d’échange pseudo-stationnaire, Q12 = λ ·(h1 − h2). Pour un milieu à double-porosité et simple-perméabilité, Warren et Root (1963) proposent de simplifier le modèle de Barenblatt et al. (1960) en considérant qu’un seul domaine participe à l’écoulement. La Fig. 1.5 présente schématiquement le concept de double-porosité. L’aquifère discontinu est composé d’un ensemble de vides, fractures et matrice. Celui-ci est conceptuellement représenté sous forme de fractures qui délimitent une superposition de blocs de matrice. (a) Schéma du réservoir réel (b) Réservoir schématisé Figure 1.5 – Modèles schématiques d’un réservoir réel à porosités multiples et représentation conceptuelle adoptée pour la modélisation. Le réservoir réel contient différentes porosités (vide, fracture, matrice) représentée Fig. 1.5a. Le modèle conceptuel du réservoir représenté Fig. 1.5b, montre la simplification en blocs à faible porosité (matrice) séparés par un ensemble de fractures. Modifié d’après Warren et Root (1963). L’écoulement est alors décrit de la manière suivante. Dans ce modèle, le premier milieu est à forte perméabilité et à faible porosité. Le second milieu est à très faible perméabilité et forte porosité. La faible perméabilité du second milieu la rend négligeable (k2 = 0 et donc K2 = 0 ) et aucun écoulement ne s’y tient (S 2 0 = 0). Ainsi, selon le modèle de Warren et Root (1963), le premier continuum assure une fonction transmissive alors que le second assure une fonction capacitive. Sous ces hypothèses le système Eq. 1.6 est simplifié et s’écrit (Warren et Root, 1963) :    Ss1 ∂h1 ∂t = ∇ ·  K1∇h1  + q1 + Q21, Ss2 ∂h2 ∂t = Q12. (1.7) Ces modèles ont été étudiés analytiquement de façon importantes et la théorie des écoulements associés à ce type de représentation a été établie en particulier pour tenir compte des pertes de charges à l’intermédiaire de chaque continuum (Moench, 1984). Analytiquement, cette représentation des écoulements dans un double continuum a été développée selon des hypothèses fortes, en particulier, de symétrie radiale et de propriété supposées homogènes. Pour la modélisation des écoulements dans le karst, Teutsch (1988, 1990) et Sauter (1992) montrent que cette adaptation de la représentation du milieu en un double milieu peut être appliquée pour modéliser les processus d’écoulement et de transport à l’échelle régionale. En comparaison avec le modèle à simple porosité, cette approche permet de mieux rendre compte du comportement hydrodynamique grâce à la prise en compte des hétérogénéités de conductivité hydraulique et d’emmagasinement qui existent entre le réseau de conduit et la roche environnante. Nous rappelons toutefois que les modèles double-porosité sont basés sur des principes d’homogénéisation des processus d’écoulement ou de transport, depuis une échelle où chaque composante hydraulique est identifiée dans l’espace jusqu’à une échelle (plus grande) où ces composantes ne sont plus identifiables individuellement (V.E.R.). Enfin, ces modèles supposent que les écoulements associés au réseau de conduit s’effectuent selon un régime d’écoulement laminaire (donc non turbulents) dans leur continuum respectif, ces deux hypothèses ne permettent pas de décrire 1.3. Modélisation et hydrodynamique karstique 19 la géométrie du champ des paramètres hydrauliques du milieu.

 Milieu hybride (approche couplée milieu discret et milieu continu)

 Les milieux hybrides permettent de prendre en compte au sein d’un milieu considéré comme poreux et homogène des singularités hydrauliques locales connues (p. ex fracture ou conduit) (Kiraly, 1998; Liedl et al., 2003). Dans un aquifère et pour une dimension donnée, toutes les singularités ne sont pas nécessairement connues à l’échelle régionale, néanmoins les principales peuvent être décrites à l’échelle locale. Ces modèles supposent donc qu’en un ou plusieurs points bien définis, l’écoulement est influencé par un objet hydraulique connu (p. ex. un conduit). Ces modèles résolvent les équations qui décrivent les écoulement dans ces objets hydrauliques, comme par exemple les équations de Saint-Venant. En général, un schéma de calcul en différences, volumes ou éléments finis est utilisé. Conceptualisant les écoulements au sein des conduits karstiques et les interactions avec le milieu fracturé environnant, de Rooij (2008), Reimann (2009) et Reimann et al. (2011) proposent des modèles d’écoulements qui rendent compte des écoulements au sein d’un conduit en charge ou à surface libre et des échanges de flux avec le milieu continu environnant. Le module C.F.P. (pour Conduit Flow Processes) inclut dans le code MODFLOW (McDonald et Harbaugh, 1984; Harbaugh, 2005) initialement développé par Shoemaker et al. (2008) permet de prendre en compte des écoulements dans un conduit en charge. Les deux milieux interagissent en fonction des différences de charge entre chaque milieu. Ce module (C.F.P.) a récemment été modifié (C.F.P.2) dans le but d’interpréter un essai par pompage mené sur la source karstique des Cent-Fonts en France (Reimann et al., 2013). Une modification majeure apportée à la version initiale est l’introduction d’un paramètre d’emmagasinement associé au domaine de conduit. Ce nouveau paramètre permet en particulier d’éviter qu’une condition à la limite appliquée au conduit ne se propage directement dans celui-ci. Ainsi les effets d’amortissement de la variation de pression liés à l’emmagasinement des flux pour des temps de pompages courts sont pris en compte. Najib (2007) propose de modéliser le flux produit à chaque fracture dans des drains artificiels forés sous un bâtiment disposé sur l’épikarst. A l’intersection de chaque fracture avec les drains artificiels est appliqué une solution analytique de puits maintenu à chaque constante. Chaque débit produit, pour chaque fracture, est calculé à l’aide d’une solution analytique de milieux à double porosité. Les charges piézométriques sont calculées en appliquant le principe de superposition. L’association de la résolution analytique du problème et de la prise en compte de conditions aux limites spécifiques induit l’utilisation d’un ensemble de puits images, disposés sur le pourtour du domaine. La solution analytique utilisée pour chaque puits maintenu à charge imposée est basée sur les travaux de Hamm (1994). Bien que des débits ainsi que des hauteurs piézométriques puissent être calculées avec cette approche analytique, les résultats souffrent particulièrement des imprécisions associées à la méthode (puits images et principe de superposition) limitée qui est utilisée pour contraindre les conditions aux limites du modèle. 

 Réseaux de fractures/conduits discrets 

En supposant que les transferts de flux s’opèrent uniquement dans les structures les plus perméables (réseaux de conduits ou de fractures) d’un aquifère, une représentation plus réaliste des écoulements peut-être obtenue (Andersson et Dverstorp, 1987; Long et al., 1982; De Dreuzy et al., 2002). Les principales contraintes liées à ces approches sont relatives à la localisation et au choix du modèle mathématique employé pour décrire l’écoulement dans le milieu. Snow (1969) définit le flux traversant une fracture plane par une loi cubique (la perméabilité dépend du cube de 20 Chapitre 1. Eléments d’hydrogéologie karstique l’ouverture de la fracture). Pour décrire les écoulements dans le karst, les lois privilégiées devraient rendre compte des régimes d’écoulement laminaire ou turbulent et des pertes de charges associées. Les lois d’écoulement communément utilisées dans ces cas sont les lois de Darcy-Weisbach, de Hagen-Poiseuille, de Louis (1968) ou des approximations des équations de Saint-Venant (Diersch, 2002; Reimann et al., 2011). Atkinson (1977) à partir de la formulation du flux proposée par la loi de Darcy-Weisbach propose de calculer les pertes de charges quadratiques dans les conduits karstiques. Jeannin et Maréchal (1995) proposent une revue et une comparaison détaillées des lois et paramètres d’écoulement en conduite. Jeannin (2001) propose de modéliser le comportement hydraulique du système épikarstique de la grotte de Hollöch (Suisse) à l’aide du modèle de Louis (1968), l’auteur montre notamment qu’il est possible de calculer les pertes de charge associées aux conduits. Bien que ces modèles puissent permettre de modéliser le comportement hydraulique de réseaux de fractures ou de conduits, ceux-ci restent néanmoins limités par la nécessité de caractériser précisément (i) les paramètres hydrauliques qui régissent les écoulement au sein du milieu et (ii) la localisation des objets (fissures, conduits) considérés dans la modélisation du système karstique. 

 Méthodes numériques de résolutions des équations d’écoulement

 La modélisation à l’aide des modèles distribués nécessitent d’utiliser des méthodes de résolution numérique, que ce soit pour résoudre les équations de transport ou d’écoulement. L’utilisation de telles méthodes plutôt que des solutions analytiques appropriées peut-être résumée ainsi (De Marsily, 1986) : — lorsque le domaine de résolution est borné et de forme complexe, la méthode des images ne peut-être utilisée avec une solution analytique (p. ex. trop complexe pour représenter la limite), — si le problème est non linéaire (p. ex. transmissivité variable en aquifère libre), — les propriétés hydrauliques varient dans l’espace selon des géométries complexes, — l’expression de la solution analytique est connue mais sa résolution requiert un important effort de calcul. Pour résoudre les E.D.P. qui décrivent les écoulements de fluides en milieu souterrain, différentes méthodes de résolution numérique sont communément utilisées. On peut notamment citer les méthodes suivantes : — éléments finis (M.E.F., FEM en anglais), — différences finies (M.D.F., FDM en anglais), — volumes finis (M.V.F., FVM en anglais), — éléments et intégrales de frontières (M.E.I.F., BEM en anglais), — éléments analytiques (M.E.A., AEM en anglais). Toutes ces méthodes nécessitent une discrétisation d’un domaine de calcul. A l’exception de la M.E.I.F. et de la M.E.A, qui ne nécessitent qu’une discrétisation de la frontière du domaine, toutes les autres méthodes nécessitent de discrétiser une surface ou un volume du domaine. En hydrogéologie, les codes de calculs à mailles couramment utilisés pour les calculs de flux (et de potentiel) sont par exemple : MODFLOW (McDonald et Harbaugh, 1984) qui utilise la méthode des différences finies, BIGFLOW (Ababou et Bagtzoglou, 1993), FEFLOW (Diersch, 2005), HydroGeoSphere (Therrien et al., 2006) qui utilisent la méthode des éléments finis. Certains de ces codes incorporent aussi la résolution des équations de transport de masse. De Marsily (1986) indique qu’il n’existe pas de critère précis pour décider quelle méthode numérique employer par rapport à une autre. Toutefois, certaines considérations (p. ex. les limites, la précision) peuvent orienter le choix vers une des méthodes numériques citées. Les méthodes de résolution qui nécessitent une discrétisation de la surface ou du volume en de multiples sous-éléments (p. ex. M.E.F, M.D.F., M.V.F.) résolvent l’équation de diffusivité en 1-D, 2-D ou 3-D. Les M.D.F. sont généralement limitées par la forme des mailles associées à la discrétisation du domaine (c.-à-d. maillages simples carrés ou rectangulaires). Les M.E.F. peuvent être complexes à mettre en oeuvre (code) et à utiliser (fichiers d’entrée spécifiques et maillage approprié). Bien qu’elles permettent de discrétiser un domaine précisément selon des méthodes de maillages particulières (p. ex. mailles triangulaires), les systèmes d’équations résolues peuvent devenir complexes et des erreurs de convergence peuvent apparaître pour certains types de calcul (p. ex. en régime transitoire) (Brezzi et Fortin, 1991).