Caractérisation stochastique des hétérogénéités hydrostratigraphiques appliquée à la modélisation hydrogéologique régionale

L'une des principales sources d'incertitude des modèles hydrogéologiques porte sur les hétérogénéités hydrostratigraphiques. Ces dernières jouent un rôle majeur sur l'écoulement de l'eau souterraine et le transport des contaminants. Pourtant la caractérisation hydrogéologique inclut très rarement l'évaluation de l'effet de l'incertitude des hétérogénéités. Cette situation s'explique en partie par la difficulté à reproduire, de manière satisfaisante, ces caractéristiques physiques en modèles équivalents. Peu ou pas de méthodologie complète et efficace permet d'effectuer la modélisation hydrogéologique régionale dans un cadre stochastique tout en quantifiant l'incertitude des hétérogénéités inter et intra unités. La difficulté de simuler efficacement les hétérogénéités inter unités de ce type de contexte est un frein à la modélisation hydrogéologique stochastique régionale. Également, la difficulté associée à la paramétrisation stochastique régionale intra unités (conductivité hydraulique ou K ) de ces modèles représente aussi un obstacle. Ces difficultés constituent une problématique majeure quand vient le besoin de quantifier adéquatement l'incertitude des modèles. Les méthodes de simulation géostatistique moderne, telles que la méthode plurigaussienne et multipoint, permettent de produire des modèles équivalents des hétérogénéités inter unités, appelés réalisations. Malgré leurs récents développements, ces méthodes conservent différentes lacunes telles que le manque de réalisme géologique dans certains contextes. Par exemple, la méthode plurigaussienne ne peut facilement incorporer des transitions asymétriques entre unités, de même que la méthode multipoint qui éprouve aussi des difficultés avec la non- stationnarité ainsi que les forts contrastes d'épaisseurs. Par contre, la méthode MCP peut s'affranchir de ces difficultés notamment grâce au transiogramme et sa propriété de forçage 0/1 des probabilités. Les hétérogénéités intra unités des modèles hydrogéologiques sont définies à partir des don- nées de K scalaire à une échelle locale. Cette échelle n'est généralement pas représentative de la taille des éléments d'un modèle hydrogéologique régional. La caractérisation stochastique de la K équivalente (tenseur-K ) à l'échelle de l'élément est un défi en soi. Elle nécessite plu- sieurs étapes incluant des contraintes à respecter dont la préservation de la structure spatiale et la contrainte d'inégalité Kverticale ≤ Khorizontale des tenseurs ainsi que les corrélations inter composantes. À cela s'ajoute, la mise à l'échelle de la K qui tient compte de l'effet des hétérogénéités locale régionalisée, la régionalisation des composantes corrélées des tenseurs-K des différentes unités ainsi que leur calage avec les variables d'état (p. ex. charge hydraulique) dans un contexte multivariable contraint. Une méthodologie complète devrait également fournir des réalisations significativement différentes afin de quantifier adéquatement l'incertitude associée aux tenseurs-K. L'objectif général de la thèse consiste à développer une démarche complète de caractérisation stochastique de l'incertitude des hétérogénéités inter et intra unités hydrostratigraphiques appliquées à la modélisation hydrogéologique régionale. Plus spécifiquement, la thèse vise à développer une méthodologie de simulation géostatistique des unités hydrostratigraphiques (UHS) afin de fournir des modèles équivalents d'aspect réaliste pour un contexte de sédimentation directionnelle (asymétrie) et à développer une méthodologie de caractérisation stochastique, de mise à l'échelle et de calage des tenseurs-K équivalents à partir de la K scalaire tout en préservant la structure spatiale des tenseurs, la contrainte d'inégalité ainsi que les corrélations inter composantes. La méthodologie globale s'oriente autour de trois étapes principales : i) la compilation et la préparation des données, ii) la simulation de modèles hydrostratigraphiques directionnels et iii) la caractérisation stochastique des tenseurs-K. La première étape a porté sur la compilation et le traitement des données (hydrostratigraphie et K) pour les développements méthodologiques et les différents tests méthodologiques appliqués sur le sous-bassin Innisfil Creek (comté de Simcoe, Ontario). La deuxième étape est basée sur le développement d'une démarche de simulation par la méthode MCP afin de générer des modèles équivalents dans un contexte hydrostratigraphique 3D complexe dont la séquence des unités est ordonnée verticalement. La méthodologie proposée inclut l'estimation des probabilités bivariables par transformée de Fourier à partir du modèle déterministe et des critères de qualité de simulation. La méthode est testée pour montrer sa capacité à quantifier la variabilité locale de la solution et son incertitude. La troisième étape fournit une méthodologie complète de caractérisation stochastique des hétérogénéités intra UHS qui inclut : la mise à l'échelle d'hétérogénéités locales (K) en tenseur-K équivalent, l'utilisation du modèle linéaire de corégionalisation pour la simulation des composantes des tenseurs-K équiprobables, la préservation des corrélations entre les composantes des tenseurs-K par nouvelle transformation gaussienne bivariée et l'utilisation de la méthode de déformation graduelle (GDM) pour caler les tenseurs-K en fonction des observations de charge hydraulique. Les résultats montrent la capacité de la méthode MCP à simuler efficacement un contexte hydrostratigraphique régional 3D complexe avec de multiples unités. De plus, les résultats montrent la capacité de la méthode à reproduire l'ordination verticale des unités et tout ça malgré la forte hypothèse d'indépendance conditionnelle. Les résultats des simulations conditionnelles et non conditionnelles pour le cas synthétique montrent clairement l'absence de tout biais substantiel. Dans l'exemple plus complexe du contexte réel, le biais est moins important pour le cas non conditionnel et pour l'échantillonnage aléatoire des données que pour le biais observé dans les données dû à l'échantillonnage préférentiel. La surreprésentation de données conditionnantes n'a pas d'impact significatif sur les proportions moyennes simulées mais l'ajout de données supplémentaires réduit la variabilité. Soulignons que la méthode MCP fournit de bons résultats malgré une image d'entraînement fortement non stationnaire. Avec une autre méthode (PGS ou MPS), rien n'aurait garanti un tel succès dans le contexte étudié. Les résultats de la méthodologie de caractérisation de la K montrent que la nouvelle méthode proposée fournit des modèles équivalents calés des tenseurs-K. La méthode de mise à l'échelle de la K scalaire en tenseur-K équivalent tient compte de la structure des hétérogénéités locales et des propriétés statistiques de chacune des UHS. La régionalisation des composantes des tenseurs-K effectuée par simulation géostatistique à partir de modèles linéaires de corégionalisation (LCM) rend le temps de calcul raisonnable. La simulation reproduit aussi, pour chacune des unités, les caractéristiques statistiques des différentes composantes de K. Pour rétablir les corrélations non linéaires des tenseurs-K non-gaussiens ainsi que les histogrammes des composantes souhaités, une nouvelle transformation bivariée est appliquée sur les réalisations. Les résultats du calage des tenseurs-K montrent que la méthode des déformations graduelles permet le calage multivariable contraint tout en préservant la covariance spatiale et croisée de UHS. Il a été démontré que la méthode conserve les relations non linéaires entre les composantes du tenseur-K grâce à l'intégration de la transformation bivariée. L'approche proposée a ainsi permis de préserver les propriétés spatiales des tenseurs ainsi que les relations entre les composantes. Les résultats démontrent que la GDM s'applique à un contexte multivariable contraint et pas seulement au contexte univariable. Soulignons que les inversions convergent rapidement en termes du nombre d'itérations grâce à l'optimisation du paramètre de déformation de la GDM. Il a aussi été montré que la méthode permet de caler la recharge séparément ou simultanément avec les tenseurs-K. Les exemples de calage montrent qu'il est avantageux d'inclure la recharge stochastique dans l'inversion. Par ailleurs, il a été observé qu'un nombre élevé d'observations de charge n'améliore pas significativement les résultats. Une comparaison, portant sur la définition de zones de captage autour de deux puits de pompage, entre la méthode GDM proposée et la méthode d'inversion PEST montrent des tenseurs-K géologiquement plus réalistes avec GDM que ceux obtenus avec PEST. L'ensemble des résultats permet de conclure que la méthode MCP simule efficacement un contexte régional complexe d'une succession stratigraphique directionnelle d'un bassin sédimentaire glaciaire. La méthode prend efficacement en charge la forte non-stationnarité spatiale des unités ainsi que d'imposantes contraintes comme la déposition verticale ordonnée en tenant compte de l'asymétrie stratigraphique à partir des probabilités de transition entre les unités. La méthode MCP s'est révélée sans biais dans le cas non conditionnel. La méthodologie proposée de caractérisation des tenseurs-K permet le calage individuel ou simultané des tenseurs-K de différentes UHS et/ou de recharge. La méthodologie inclut une mise à l'échelle des conductivités hydrauliques quasi ponctuelles vers un tenseur équivalent en tenant compte de la structure spatiale à l'échelle locale. Bien que la mise à l'échelle soit intensive en temps de calcul, il a été démontré qu'il est possible de limiter au minimum le nombre d'éléments à mettre à l'échelle grâce à la régionalisation des tenseurs-K par simulation géostatistique. La thèse a apporté quelques contributions originales. D'abord, il s'agit d'une première appli- cation de la méthode MCP pour la simulation 3D. À ma connaissance, aucune autre méthode de simulation stochastique ne semble être en mesure de solutionner avec autant de satisfaction un contexte hydrostratigraphique complexe 3D comme celui du comté de Simcoe. La méthode des déformations graduelles permet le calage efficace des composantes des tenseurs- K avec préservation des corrélations non linéaires des composantes et de leurs fonctions de covariance. À ma connaissance, il s'agit de la première approche utilisant la GDM pour un contexte multivariable avec contraintes. Enfin, cette thèse contribue à l'avancement des connaissances et des techniques de la caractérisation stochastique des hétérogénéités inter et intra hydrostratigraphiques. Dans une plus grande perspective, la thèse contribue à fournir des outils adaptés pour mieux quantifier l'incertitude associée aux hétérogénéités présente dans les modèles hydrogéologiques régionaux.

Abstract

One of the main sources of uncertainty in hydrogeological models relates to hydrostrati- graphic heterogeneities. They play a major role in quantification of groundwater dynamics and mass transport. Yet hydrogeological characterization rarely includes the assessment of the effect of uncertainty and spatial variability of hydrostratigraphic units and hydrogeologic parameters. This is partly explained by the difficulty of reproducing these physical charac- teristics satisfactorily in equivalent models. There is at present little or no comprehensive and efficient methodology that allows regional hydrogeological modelling to be carried out in a stochastic framework while quantifying the uncertainty of inter- and intra-unit hetero- geneities. The difficulty to effectively simulate inter-unit heterogeneities under such condi- tions is a hindrance to stochastic regional hydrogeological modelling. As well, the difficulty associated with the regional stochastic intra-unit parameterization (hydraulic conductivity or K ) of these models further aggravates the problem. These difficulties constitute a major issue when it comes to the need to adequately quantify model uncertainty. Modern geostatistical simulation methods, such as plurigaussian and multipoint simulations, make it possible to produce equivalent models of inter-unit heterogeneities, referred to as realizations. Despite their recent developments, these methods retain various shortcomings such as the lack of geological rigour in certain settings. For example, the plurigaussian method cannot easily incorporate asymmetrical transitions between units, whereas the multipoint method has dif- ficulties with non-stationarity and strong thickness contrasts. On the other hand, the MCP (Markov-type categorical prediction) method has the capacity to overcome these difficulties thanks to the transiogram and its 0/1 probability forcing property. The intra-unit heterogeneities of the hydrogeological models are defined from scalar K data at a local scale. This scale is generally not representative of the size of the elements used in regional hydrogeological models. Stochastic characterization of the equivalent K (K -tensor) at the element scale is a challenge in itself. It requires several steps including constraints to be respected such as the preservation of the spatial structure, the natural constraint Kvertical ≤ Khorizontal, and the inter-component correlations. In addition, the upscaling of K takes into account the effect of local heterogeneities, the regionalization of the correlated components of the K -tensors of the different units and their calibration with state variables (e.g., hydraulic head) in a constrained multivariate context. A comprehensive methodology should also provide significantly different realizations in order to adequately quantify the uncertainty associated with the K -tensors. The general objective of this thesis is to develop a comprehensive approach to stochastic characterization of the uncertainty of inter- and intra-hydrostratigraphic unit heterogeneities applied to regional hydrogeological modelling. More specifically, the thesis aims to develop a methodology for geostatistical simulation of hydrostratigraphic units (HSU) in order to provide equivalent models of realistic appearance for a directional sedimentation context (asymmetry) and to develop a methodology for stochastic characterization, upscaling and calibration of equivalent K -tensors from the scalar K while preserving the spatial structure of the tensors, the inequality constraint as well as the inter-component correlations. The overall approach revolves around three main steps: (i) data compilation and preparation, (ii) simulation of directional hydrostratigraphic models, and (iii) stochastic characterization of K -tensors. The first step focuses on data compilation and processing (hydrostratigraphy and K ) for methodological developments and various methodological tests applied to the Innisfil Creek sub-watershed (Simcoe County, Ontario). The second stage is based on the development of a simulation approach using the MCP method to generate equivalent models in a complex 3D hydrostratigraphic context with a vertically ordered sequence of units. The proposed methodology includes the estimation of bivariate probabilities by Fourier transform from the deterministic model and the simulation quality criteria. The method is tested to show its ability to quantify the local variability of the solution and its uncertainty. The third step provides a complete methodology for stochastic characterization of intra-HSU heterogeneities which includes: upscaling of local heterogeneities (K ) into equivalent K - tensors, use of the linear co-regionalization model for simulation of equiprobable K -tensor components, preservation of correlations between K -tensor components by new bivariate Gaussian transformation and use of the Gradual Deformation Method (GDM) to calibrate the K -tensors according to hydraulic head observations. The results show the ability of the MCP method to effectively simulate a complex 3D regional hydrostratigraphic with multiple units. Furthermore, the results show the ability of the method to reproduce the vertical ordering of the units, despite the strong assumption of conditional independence. The results of the conditional and unconditional simulations for the synthetic case clearly show the absence of any substantial bias. In the more complex example of field hydrogeological settings, the bias is less important in the unconditional case and in the random sampling case than the bias observed in the data due to preferential sampling. Over-representation of conditioned data does not have a significant impact on the mean simulated proportions, but the addition of supplementary data reduces variability. The MCP method provides good results despite a highly non-stationary training image. With another method (PGS or MPS), nothing would have guaranteed such success in the studied context. The results of the K characterization methodology show that the presented original method provides equivalent calibrated models of the K -tensors. The scalar K upscaling method to equivalent K -tensor takes into account the structure of local heterogeneities and the sta- tistical properties of each of the HSU. The regionalization of the K -tensor components by geostatistical simulation using Linear Co-Regionalization Models (LCMs) keeps the compu- tation time reasonable. The simulation also reproduces, for each of the units, the statistical characteristics of the different K components. To restore the non-linear correlations of the non-Gaussian K -tensors as well as the histograms of the desired components, a new bivariate transformation is applied to the realizations. The results of the calibration of the K -tensors show that the method of gradual deformations allows the multivariate constrained calibration while preserving the spatial and cross-covariance of HSU. The method has been shown to preserve the non-linear relationships between the K -tensor components through the integra- tion of the bivariate transformation. The proposed approach has thus preserved the spatial properties of the tensors and the relationships between the components. The results show that the GDM applies to a constrained multivariate context and not only to the univariate context. It should be noted that the inversions converge rapidly in terms of the number of iterations thanks to the optimization of the deformation parameter of the GDM. It has also been shown that the method allows to set the recharge separately or simultaneously with the K -tensors. The calibration examples show that it is advantageous to include stochastic recharge in the inversion. It has also been observed that a high number of head observations does not significantly improve the results. A comparison of the definition of capture zones around two pumping wells between the proposed GDM method and the PEST inversion method shows geologically more realistic K -tensors with GDM than those obtained with PEST. The overall results suggest that the LCM method effectively simulates a complex regional context of a directional stratigraphic succession in a glacial sedimentary basin. The method effectively supports the strong spatial non-stationarity of the units as well as imposing con- straints such as ordered vertical deposition by taking into account stratigraphic asymmetry from the transition probabilities between the units. The MCP method proved to be unbiased in the unconditional case. The proposed K -tensor characterization methodology allows in- dividual or simultaneous calibration of K -tensors from different HSU and/or recharge units. The methodology includes an upscaling of quasi-point hydraulic conductivities to an equiva- lent tensor taking into account the spatial structure at the local scale. Although the upscaling is computationally intensive, it has been shown that it is possible to minimize the number of elements to be scaled by regionalizing the K -tensors by geostatistical simulation. The thesis provided substantial original contributions. This is a first application of the MCP method for 3D simulations. To the author's knowledge, no other stochastic simulation method seems to be able to solve as satisfactorily a 3D hydrostratigraphic context as complex as that of Simcoe County. The gradual deformation method allows for efficient calibration of the K - tensor components while preserving the non-linear correlations of the components and their covariance functions. To the author's knowledge, this is also the first time where GDM has been used in a multivariate context with constraints. Finally, this thesis contributes to the advancement of knowledge and techniques for stochastic characterization of inter- and intra- hydrostratigraphic heterogeneities. In a broader perspective, the thesis contributes original methods and tools adapted to better quantify the uncertainty associated with heterogeneities present in regional hydrogeological models.