Neste estudo, consideramos um problema transiente de fluxo unidimensional vertical de água em meio poroso insaturado, modelado pela equação Richards não-linear. As reações constitutivas de Van Genuchten são empregadas para representar a capacidade hidráulica e a condutividade. A fórmula da solução é otimizada e avaliada usando a equação governante em um critério de autoconsciente. Os resultados são apresentados para alguns tipos de solo e seus parâmetros relacionados, que são mencionados em literatura. / In this study, we consider a transiente vertical one-dimensional flow problem of water in unsaturated porus media, modelled by the non-linear Richards equation. Constitutive relations of Van Genutchten are employed to represent the hydraulic capacity and conductivity. The solution formula is optimized and evaluated using to governing equation for a self-consistency criterion. The results are presented for some oil types and its related soil parameters, that are reported in the literature.

Neste estudo, consideramos um problema transiente de fluxo unidimensional vertical de água em meio poroso insaturado, modelado pela equação Richards não-linear. As reações constitutivas de Van Genuchten são empregadas para representar a capacidade hidráulica e a condutividade. A fórmula da solução é otimizada e avaliada usando a equação governante em um critério de autoconsciente. Os resultados são apresentados para alguns tipos de solo e seus parâmetros relacionados, que são mencionados em literatura. / In this study, we consider a transiente vertical one-dimensional flow problem of water in unsaturated porus media, modelled by the non-linear Richards equation. Constitutive relations of Van Genutchten are employed to represent the hydraulic capacity and conductivity. The solution formula is optimized and evaluated using to governing equation for a self-consistency criterion. The results are presented for some oil types and its related soil parameters, that are reported in the literature.

O bulbo molhado formado na irrigação por gotejamento pode ser modelado a partir informações das propriedades físico-hídricas do solo. Mas apesar do grande progresso na modelagem, a aplicação de modelos na irrigação e drenagem ainda não foi implementada em nenhum nível de maneira substancial. Assim, o objetivo do presente trabalho foi desenvolver um modelo numérico utilizando a técnica dos volumes finitos para estimar a forma e as dimensões do volume de solo molhado sob irrigação por gotejamento, e verificar a validade do modelo através da comparação com dados recolhidos em condições experimentais. O ensaio foi conduzido na área experimental do Departamento de Engenharia de Biossistemas da Escola Superior de Agricultura Luiz de Queiroz-ESALQ/USP. As sondas de TDR foram confeccionadas segundo procedimentos descritos por Souza et al. (2006). O perfil do bulbo molhado formado no solo abaixo de um emissor do tipo gotejador foi determinado para 3 tipos de material: areia, solo arenoso e solo argiloso e para 3 vazões nominais do emissor: 2, 4 e 6 L/h. Cada uma das combinações entre tipo de solo e vazão foi repetida três vezes em caixas (1000 L) diferentes, totalizando 27 medições. Um modelo numérico foi desenvolvido para simular a distribuição da água no perfil do solo abaixo de uma fonte pontual. A solução da equação diferencial para o movimento da água em solo não saturado foi realizada pela discretização do espaço-tempo utilizando a técnica dos volumes finitos. Um algoritmo em linguagem Visual Basic foi escrito para implementar o conjunto de equações e simular a evolução do bulbo molhado no tempo. Foi simulada a formação do bulbo molhado nas mesmas condições do experimento, e gráficos de isolinhas de umidade foram traçados no software Surfer. O modelo numérico proposto foi capaz de simular a formação do bulbo molhado em diferentes condições de tipo de solo e vazão do emissor. / The soil water patterns in drip irrigation can be modeled from physical and hydraulic properties of soil. But despite the great progress in modeling, the models have not been implemented at any level in irrigation and drainage. The objective of this study was to develop a numerical model using finite volume technique to estimate the shape and dimensions of the wetted soil volume under drip irrigation and verify the validity of the model by comparing data collected under experimental conditions. The trial was conducted at the experimental site of the Department of Biosystems Engineering, Escola Superior de Agricultura Luiz de Queiroz- ESALQ/USP. The TDR probes were manufactured according Souza et al. (2006). The soil water patterns under a drip emitter type was determined for three types of material: sand, sandy soil and clay soil and three nominal flows of the emitters: 2, 4 and 6 L/h. Each of the combinations of soil type and flow rate was repeated three times in boxes (1000 L), totaling 27 measurements. A numerical model was developed to simulate the distribution of water in the soil profile below a point source. The solution of the differential equation for the movement of water in unsaturated soil was carried out by the discretization of space-time using the technique of finite volume. An algorithm in Visual Basic language was written to implement the set of equations and simulate the evolution of wetted soil volume in time. Contour plots of soil water content were drawn in Surfer software. The proposed numerical model was able to simulate wetted soil volume under different conditions of soil type and flow of the emitter.

O bulbo molhado formado na irrigação por gotejamento pode ser modelado a partir informações das propriedades físico-hídricas do solo. Mas apesar do grande progresso na modelagem, a aplicação de modelos na irrigação e drenagem ainda não foi implementada em nenhum nível de maneira substancial. Assim, o objetivo do presente trabalho foi desenvolver um modelo numérico utilizando a técnica dos volumes finitos para estimar a forma e as dimensões do volume de solo molhado sob irrigação por gotejamento, e verificar a validade do modelo através da comparação com dados recolhidos em condições experimentais. O ensaio foi conduzido na área experimental do Departamento de Engenharia de Biossistemas da Escola Superior de Agricultura Luiz de Queiroz-ESALQ/USP. As sondas de TDR foram confeccionadas segundo procedimentos descritos por Souza et al. (2006). O perfil do bulbo molhado formado no solo abaixo de um emissor do tipo gotejador foi determinado para 3 tipos de material: areia, solo arenoso e solo argiloso e para 3 vazões nominais do emissor: 2, 4 e 6 L/h. Cada uma das combinações entre tipo de solo e vazão foi repetida três vezes em caixas (1000 L) diferentes, totalizando 27 medições. Um modelo numérico foi desenvolvido para simular a distribuição da água no perfil do solo abaixo de uma fonte pontual. A solução da equação diferencial para o movimento da água em solo não saturado foi realizada pela discretização do espaço-tempo utilizando a técnica dos volumes finitos. Um algoritmo em linguagem Visual Basic foi escrito para implementar o conjunto de equações e simular a evolução do bulbo molhado no tempo. Foi simulada a formação do bulbo molhado nas mesmas condições do experimento, e gráficos de isolinhas de umidade foram traçados no software Surfer. O modelo numérico proposto foi capaz de simular a formação do bulbo molhado em diferentes condições de tipo de solo e vazão do emissor. / The soil water patterns in drip irrigation can be modeled from physical and hydraulic properties of soil. But despite the great progress in modeling, the models have not been implemented at any level in irrigation and drainage. The objective of this study was to develop a numerical model using finite volume technique to estimate the shape and dimensions of the wetted soil volume under drip irrigation and verify the validity of the model by comparing data collected under experimental conditions. The trial was conducted at the experimental site of the Department of Biosystems Engineering, Escola Superior de Agricultura Luiz de Queiroz- ESALQ/USP. The TDR probes were manufactured according Souza et al. (2006). The soil water patterns under a drip emitter type was determined for three types of material: sand, sandy soil and clay soil and three nominal flows of the emitters: 2, 4 and 6 L/h. Each of the combinations of soil type and flow rate was repeated three times in boxes (1000 L), totaling 27 measurements. A numerical model was developed to simulate the distribution of water in the soil profile below a point source. The solution of the differential equation for the movement of water in unsaturated soil was carried out by the discretization of space-time using the technique of finite volume. An algorithm in Visual Basic language was written to implement the set of equations and simulate the evolution of wetted soil volume in time. Contour plots of soil water content were drawn in Surfer software. The proposed numerical model was able to simulate wetted soil volume under different conditions of soil type and flow of the emitter.

Les eaux souterraines constituent une réserve d’eau potable non négligeable, leur alimentation se fait en majeure partie par les précipitations, appelée recharge des nappes. Du fait de leur grande importance, la compréhension du fonctionnement de ces ressources en eau est plus que jamais indispensable. Celle-ci passe par l’élaboration de modèles mathématiques. Ces outils nous offrent une meilleure appréhension et une bonne prévision des phénomènes physiques. Les systèmes hydrogéologiques sont généralement très complexes et caractérisés par des dynamiques hydriques très variables dans le temps et dans l’espace. Cette complexité a attiré l’attention de nombreux hydrogéologues et un grand nombre de modèles très sophistiqués a été développé afin de décrire ces systèmes avec précision. Cependant, la prise en compte de la recharge de ces réservoirs reste toujours un défi dans la modélisation hydrogéologique. En effet, le plus souvent, les modèles hydrogéologiques simulent l’écoulement dans la nappe tout en considérant la recharge comme une constante sur le domaine et indépendante du système. De plus, elle est souvent calculée de façon simplifiée. Or, la recharge traduisant la quantité des précipitations atteignant les nappes est une composante hydrologique complexe et variable car elle interagit avec les nappes et dépend des conditions climatiques, du couvert végétal et du transfert de l’eau dans le sol. Ce présent travail vise à intégrer cette recharge variable et complexe aux modèles hydrogéologiques. À cet effet, un modèle couplé a été développé. Une première partie de ce modèle permet de calculer la recharge des nappes en modélisant les interactions précipitations-sol et l’hydrodynamique dans le sol. Cette modélisation a été effectuée en utilisant des modèles conceptuels simples basés sur des lois empiriques (Gardénia, Nash) et des modèles physiques résolvant l’équation de Richards. La recharge ainsi calculée est intégrée à la deuxième partie du modèle simulant l’hydrodynamique dans les nappes décrite par l’équation de diffusivité. Des méthodes numériques précises et robustes ont été utilisées pour résoudre les équations du modèle mathématique : les éléments finis non conformes ont été utilisés pour résoudre l’équation de diffusivité et l’équation de Richards est résolue sous sa forme mixte par une méthode itérative en temps. En somme, ce modèle couplé permet de décrire les variations de niveaux de nappe à partir des données météorologiques connaissant les paramètres caractéristiques de cet aquifère. [...] / Groundwater is the main available water resource for many countries; they are mainly replenished by water from precipitation, called groundwater recharge. Due to its great importance, management of groundwater resources is more essential than ever, and is achieved through mathematical models which offer us a better understanding of physical phenomena as well as their prediction. Hydrogeological systems are generally complex thus characterized by a highly variable dynamic over time and space. These complexities have attracted the attention of many hydrogeologists and many sophisticated models that can handle these issues and describe these systems accurately were developed. Unfortunately, modeling groundwater recharge is still a challenge in groundwater resource management. Generally, groundwater models are used to simulate aquifers flow without a good estimation of recharge and its spatial-temporal distribution. As groundwater recharge rates show spatial-temporal variability due to climatic conditions, land use, and hydrogeological heterogeneity, these methods have limitations in dealing with these characteristics. To overcome these limitations, a coupled model which simulates flow in the unsaturated zone and recharge as well as groundwater flow was developed. The flow in the unsaturated zone is solved either with resolution of Richards equation or with empirical models while the diffusivity equation governs flow in the saturated zone. Robust numerical methods were used to solve these equations: we apply non-conforming finite element to solve the diffusivity equation and we used an accurate and efficient method for solving the Richards equation. [...]

After the seminal works by Freeze and Harlan (1969), the scientific community realized that groundwater and vadose zone equation were breaking up. Hydrologists split into three communities following the motto “you are my boundary condition”: groundwater people, vadose zone scientists and surface water hydrologists. This compartmentalization of the scientific community fostered a deepening of knowledge in single branches, allowing to break things down into simple parts. However, this division represented an obstacle to the comprehension of the complexity that characterises the interactions between them. Eventually, this separation of the communities continued into software code. As a matter of fact, the boundary conditions were hard-wired, but they offered a poor representation of the physics in the interaction between different domains. Recently, there has been a renewed interest in studying the big picture, the interactions between different domains. This it is evident in the development of a new research field named the Earth’s Critical Zone (CZ). It is defined as the “ heterogeneous, near surface environment in which complex interactions involving rock, soil, water, air, and living organism regulate the natural habitat and determine the availability of life-sustaining resources” (National Research Council, 2001). Further interest in the studying the CZ is given by the ever-increasing pressure due to the growth in human population, wealth, and climatic changes. This thesis focuses on the CZ while recognising the central role of having a solid set of tools for modeling the water movements in all conditions. Recently, Prentice et al. (2015) identified Reliable, Robust, and Realistic, the three R’s, as the three characteristics that numerical models should have. Soil moisture is one of the key components to simulate the processes in the critical zone. The governing equation to describe the water flow in a porous material is know as the Richards equation and it dates back to 1931.The numerical solution of the Richards equation is far from trivial because of its mildly nonlinearity and it is often discarded in favour of more empirical models. After the pioneering work by Celia et al. (1990), a lot of work has been done in this direction and several model, for instance Hydrus, GEOtop, Cathy, Parflow adopted variants of the Newton algorithm to allows global convergence. Since Casulli and Zanolli (2010), anticipated by Brugnano and Casulli (2008), a new method called nested Newton has been found to guarantee convergence in any situation, even under the use of large time steps and grid sizes. The research presented in this thesis used this integration algorithm. Besides the numerical aspect, another issue was the correct definition of the boundary condition at the soil surface. As a matter of fact, the definition of the surface boundary condition is necessary to capture the generation of surface run-off. In the literature several approaches were proposed to couple surface and subsurface flow, and in this work the approach presented by Gugole (2016) has been used. The novelty regarded the discretization of the shallow water equation and the Richards equation in an unique algebraic system that was solved in a conservative manner. Richards equation was criticized from many points of view, but it is difficult to criticize its core mass conservation. The definition of the hydraulic properties of the soil, including both the soil water retention function (SWRC) and the hydraulic conductivity models, often uses simplified representation of the pore system describing it as bundle of cylindrical capillaries where the largest ones drain first and are filled last. As pointed out by Bachmann et al. (2002), “physical effects, like surface water film adsorption, capillary condensation and surface flow in liquid films, as well as volumetric changes of the pore space are often ignored”. Thus, the capillary bundle concept is a rough, even if still useful approximation of soil reality. From these observations, during the research the code has been designed to offer the opportunity to easily implement new soil water hydraulic models that might be proposed in the future. The Richards’ equation alone is not anymore sufficient to model the water flow in soils. In fact, soil temperature affects the water flow in soils. This is evident in cold regions where soil water is subject to freezing and thawing processes, but also in unfrozen soil, where temperature modifies water properties such as viscosity, the surface tension, and the contact angle. These microscopic variations of the water physical properties have significant impacts in the mass and energy budget within the CZ. For instance, it has been observed that the infiltration rates between the stream and the vadose zone show a clear diurnal pattern: infiltration rates are highest in late afternoon, when stream temperature is greatest, and they are lowest in early morning when stream temperature is least. In cold regions the run-off production is strongly affected by the presence of ice with the soil. Nonetheless, soil moisture modifies the thermal properties of the soil: water is characterised by a high thermal inertia and the thermal conductivity of ice is almost four times larger than that of liquid water, and water flow carries a significant amount of sensible heat. These aspects come under one the R of realistic. Hence, the Richards’ equation has been coupled with the energy equation for the unfrozen case. Moreover, the research developed a model to study the heat transfer considering the phase change of water. In both cases robust numerical schemes have been used. There are few models that already coupled the equations. One of these models is GEOtop that was conceived and built in the research group where this work was carried out. Such models have some limitations. One of the main limitations regards their implementations. In fact, these models were built as a monolithic code and this turns in difficulties in maintaining and developing existing codes. In this work the codes have been developed by using Design Patterns. As a result, the codes are easy to maintain, to extend, and to reuse. Considering the CZ, these aspects are of crucial importance. Researchers should have a model that can be extended to include more processes, i.e. increase its complexity and avoiding the code to become too complicated. The models were integrated in the Object Modelling System v3 (OMS3) framework. The system provides various components for precipitation treatment, radiation estimation in complex terrain, evaporation and transpiration that can be connected to each other’s for generating inputs and outputs. Due to the modularity of the system, whilst the components were developed and can be enhanced independently, they can be seamlessly used at run time by connecting them with the OMS3 DSL language based on Groovy. OMS3 provides the basic services and, among them, tools for calibration and implicit parallelization of component runs. In sum, the thesis analyses the relevant literature to date. It presents a detailed description of the physical processes related to the water flow and the energy budget within the soil. Then, it describes the numerical method used to solve and coupled the equations. It also provides the informatics behind WHETGEO 1D (Water HEat Tracers in GEOframe). Finally, the work focuses on the WHETGEO extension for the bidimensional case by showing how the code can be designed to store grid information.

Devido aos riscos de contaminação dos recursos naturais solo e água, ao alto custo, ao tempo e ao esforço humano nas investigações de campo, os modelos matemáticos, aliados às técnicas numéricas e aos avanços computacionais, constituem uma ferramenta importante na previsão do deslocamento de solutos, contribuindo assim, para o controle de alterações ambientais. No Brasil, a modelação de fluxo e transporte de solutos na zona não-saturada é voltada, quase que exclusivamente, aos problemas relacionados às atividades agrícolas. Entretanto, tão importante quanto a problemática dos produtos químicos nas atividades agrícolas é a questão de poluição e contaminação do solo e da água por chorume, gerado pelos resíduos sólidos domiciliares. Neste trabalho, é desenvolvido e validado um modelo computacional unidimensional para simulação de fluxo e transporte de solutos na zona não-saturada do solo. O modelo matemático é dado pela equação diferencial parcial não-linear de Richards, que rege o movimento de água no solo, e a equação diferencial parcial linear de advecção-dispersão, do transporte de solutos, acompanhadas das condições iniciais e de contorno. A equação de Richards é dada em função do potencial matricial da água e a equação de transporte de solutos estima a evolução temporal da concentração de solutos no perfil do solo. Devido à dificuldade de se obter soluções analíticas destas equações, são resolvidas numericamente pelo método de elementos finitos. As referidas equações são resolvidas utilizando-se malhas uniformes inicialmente. Com a finalidade de obter simulações mais eficientes, a um custo computacional reduzido, é empregada a adaptatividade com refinamento h na malha de elementos finitos. A função interpolação polinomial utilizada é de grau 2 ou maior que garante a conservação de massa. Na equação de Richards, a derivada temporal é aproximada por um quociente de diferença finita e é aplicado o esquema de Euler explícito e na equação de advecção-dispersão, é aproximada por um quociente de diferença finita, aplicando-se o esquema de Euler implícito, devido à linearidade da equação. O sistema operacional é o Linux Ubuntu 32 bits, o ambiente de programação é o PZ, escrito em linguagem de programação C++. Na validação do modelo, utilizam-se dados disponíveis na literatura. Os resultados são comparados, utilizando-se malhas uniformes e malhas adaptativas com refinamento h. Usando-se as malhas uniformes para o problema de Richards e de transporte de potássio, o tempo de execução é de 22 minutos e a memória utilizada de 6164 Kb. Com as malhas adaptadas, o tempo de execução é de 3 minutos e 27 segundos, consumindo 5876 Kb de memória. Houve, portanto, uma redução de 84,32% no tempo de execução, usando-se malhas adaptativas. A utilização da função interpolação polinomial de grau 2 ou maior e o refinamento h, permitem uma boa concordância do modelo na comparação com soluções disponíveis na literatura. / Due to the risks of contamination of soil and water resources, the high cost, time and human effort in the field investigations, the mathematical models, combined with numerical techniques and computational advances, are important tools in forecasting the movement of solutes thereby contributing to the control of environmental alteration. In Brazil, modeling of flow and solute transport in the unsaturated zone is focused, almost exclusively, on problems related to agricultural activities. However, as important as the problematical of chemicals products in agricultural activities is the issue of pollution and contamination of soil and water by leachate, generated by municipal solid wastes. In this work, an one-dimensional computational model for simulation of flow and solute transport in the unsaturated soil has been developed and validated. The mathematical model is given by the Richards\'s non-linear partial differential equation, which determines the movement of water in the soil, and the advection-dispersion linear partial differential equation, of the solute transport, together with initial and boundary conditions. The Richards equation is a function of the water pressure head and the solute transport equation estimate the temporal evolution of the solutes concentration in the soil profile. Due to the difficulty of obtaining analytical solutions of these equations, they are solved numerically using the finite element method. The governing equations are solved using initially a uniform mesh. In order to obtain more efficient simulations with low computational cost, adaptativity with h refinement on the finite element mesh is implemented. The interpolation function is of degree two or higher, assuring mass conservation. In Richards\' equation, the temporal derivative is approximated by Euler explicit finite difference. For the advection-dispersion equation, due to the linearity of the equation, an implicit finite difference scheme is used. The code is written in the programming language C++ based on the programming environment PZ using operating system Linux Ubuntu 32 bit. Model results are validated in comparison with data available in the literature. The results are evaluated using uniform meshes and with h refinement adaptive mesh. Using the uniform meshes for the problem of Richards and transport of potassium, the running time is 22 minutes and 6164 Kb of memory is used. With the adapted meshes, the execution time is 3 minutes and 27 seconds, consuming 5,876 Kb of memory. Therefore there was a reduction of 84.32% in execution time, using adaptive meshes. The interpolation function with degree two or higher and the h refinement, with reduction of the computation time, showed a good agreement in comparison with the literature.

Coordenação de Aperfeiçoamento de Pessoal de Nível Superior / A engenharia geotécnica é uma das grandes áreas da engenharia civil que estuda a interação entre as construções realizadas pelo homem ou de fenômenos naturais com o ambiente geológico, que na grande maioria das vezes trata-se de solos parcialmente saturados. Neste sentido, o desempenho de obras como estabilização, contenção de barragens, muros de contenção, fundações e estradas estão condicionados a uma correta predição do fluxo de água no interior dos solos. Porém, como a área das regiões a serem estudas com relação à predição do fluxo de água são comumente da ordem de quilômetros quadrados, as soluções dos modelos matemáticos exigem malhas computacionais de grandes proporções, ocasionando sérias limitações associadas aos requisitos de memória computacional e tempo de processamento. A fim de contornar estas limitações, métodos numéricos eficientes devem ser empregados na solução do problema em análise. Portanto, métodos iterativos para solução de sistemas não lineares e lineares esparsos de grande porte devem ser utilizados neste tipo de aplicação. Em suma, visto a relevância do tema, esta pesquisa aproximou uma solução para a equação diferencial parcial de Richards pelo método dos volumes finitos em duas dimensões, empregando o método de Picard e Newton com maior eficiência computacional. Para tanto, foram utilizadas técnicas iterativas de resolução de sistemas lineares baseados no espaço de Krylov com matrizes pré-condicionadoras com a biblioteca numérica Portable, Extensible Toolkit for Scientific Computation (PETSc). Os resultados indicam que quando se resolve a equação de Richards considerando-se o método de PICARD-KRYLOV, não importando o modelo de avaliação do solo, a melhor combinação para resolução dos sistemas lineares é o método dos gradientes biconjugados estabilizado mais o pré-condicionador SOR. Por outro lado, quando se utiliza as equações de van Genuchten deve ser optar pela combinação do método dos gradientes conjugados em conjunto com pré-condicionador SOR. Quando se adota o método de NEWTON-KRYLOV, o método gradientes biconjugados estabilizado é o mais eficiente na resolução do sistema linear do passo de Newton, com relação ao pré-condicionador deve-se dar preferência ao bloco Jacobi. Por fim, há evidências que apontam que o método PICARD-KRYLOV pode ser mais vantajoso que o método de NEWTON-KRYLOV, quando empregados na resolução da equação diferencial parcial de Richards. / Geotechnical Engineering is the area of Civil Engineering that studies the interaction between constructions carried out by man or natural phenomena with geological environment, which most of times is partially saturated soil. In this sense, work developing as stabilization, dam containing, retaining walls, foundations and highways are conditioned to a right prediction of water flow into the soil. However, considering the water flow, the studied region areas are commonly on the order of square kilometers, mathematical models solutions require computational meshes of large proportions, causing serious limitations linked to computational memory requirements and processing time. In order to overcome these limitations, efficient numerical methods must be used in the solution of the considered problem. Hence iterative methods for solving nonlinear and large sparse linear systems must be used in this type of application. In short, this study approached a solution to the Richard partial differential equation by the two dimensions finite volume method, bringing Picard and Newton method with greater efficiency. Linear system resolution iterative techniques based on Krylov space with pre-conditioners matrix were used. Portable Extensible Toolkit for Scientific Computation (PETSc) numerical library was a tool used during the task. The results indicate when a Richards equation is solved considering thr PICARD-KRYLOV method, no matter the soil evaluation model, the best combination for solving linear systems is the stabilized double gradient method and the SOR preconditioning. On the other hand, when the van Genuchten equations are used the gradients methods with the SOR preconditioning must be chosen. Adopting the NEWTON-KRYLOV method, the stabilized double gradient method is more efficient in soling Newton linear system, in relation to the preconditioning it must be giving preference to the Jacob block. Finally, there are strong indications that the PICARDKRYLOV method can be more effective than the NEWTON-KRYLOV one, when used for solving Richards partial differential equation.

Dans un champ captant, comme celui qui alimente l'agglomération lyonnaise, le fonctionnement de chaque « objet » (bassins d'infiltration, puits, rivières, nappe, zone non saturée) et leurs interactions sont complexes et mal connus. Dans un premier temps, une série d'essais d'infiltration à différentes échelles dans un bassin artificiel de réalimentation a permis de mieux appréhender le fonctionnement de cet ouvrage et de chiffrer les paramètres hydrodynamiques de la zone non saturée. Les résultats des interprétations par méthodes numériques fondées sur l'équation de Richards ont montré que les flux infiltrés dans les bassins dépendent principalement de la conductivité hydraulique à saturation d'une couche située directement sous le sable calibré couvrant le fond du bassin, identifiée comme étant moins perméable que la nappe. Cette couche conditionne l'existence d'une zone non saturée. La réalisation d'essais de pompage dans l'aquifère sur les forages d'exploitation et sur un dispositif spécialement mis en place durant ce travail a permis de déterminer les paramètres hydrodynamiques de la nappe. Une analyse des observations et une modélisation conceptuelle en 2D, puis en 3D ont permis d'identifier les mécanismes prépondérants (stratifications, apports et prélèvements) et de simuler correctement à la fois les flux infiltrés dans un des bassins d'infiltration et la remontée de la nappe. A l'échelle d'un bassin, les flux infiltrés sont variables dans le temps, ils dépendent de l'état de colmatage de la surface d'infiltration mais également de la température de l'eau infiltrée et de l'état hydrique initial du sol sous le bassin. Les analyses de sensibilité réalisées avec les modèles mis en place indiquent que la conductivité hydraulique à saturation de l'aquifère, mais également la proximité des conditions aux limites imposées dans la nappe (les rivières et les puits de pompage) influencent de manière prépondérante la remontée de la nappe. Une modélisation 3D d'un autre secteur du champ captant comprenant deux bassins d'infiltration, deux bras de rivière ainsi que des puits de pompage a été réalisée. La condition à la limite imposée sur les rivières est du troisième type en accord avec l'observation d'un sous-écoulement en nappe. Les échanges nappe/rivières sont calés sur des observations à partir d'une chronique de propagation d'une onde de crue dans la nappe. Des piézomètres en flûte de pan, spécifiquement implantés à proximité d'un bassin, ont permis d'observer des différences de charge hydraulique fortes à différents niveaux de l'aquifère lorsque le bassin d'infiltration est en eau. La modélisation 3D est conforme à ces observations. Elle a confirmé l'importance du rôle d'une hétérogénéité de type argilo-sableuse (de conductivité hydraulique à saturation inférieure aux autres lithologies présentes dans l'aquifère) dans les écoulements (direction et flux). Le modèle développé représente correctement les flux infiltrés via les bassins ainsi que les fluctuations de la piézométrie de la nappe. Il permet de vérifier l'inversion des écoulements par rapport aux infiltrations de la rivière, d'identifier les puits alimentés par les bassins d'infiltration et également de mettre en évidence les flux de nappe passant sous la rivière. / In a well field of the Lyon metropolitan area, designed for drinking water supply, behaviour of each object (infiltration basins, wells, rivers, aquifer, and unsaturated zone) and their interactions are complex and not well-known. As a first step, infiltration tests at different spatial scales in one artificial basin were performed to better understand the basin operation and to estimate the hydrodynamic parameters of the unsaturated zone. Results of interpretation, using numerical methods based on Richards equation, reveal that infiltrated basin fluxes mainly depend on saturated hydraulic conductivity of a layer located just below the calibrated sand layer that cover the basin bottom. Indeed this layer has been estimated to be less permeable than the aquifer, which allows the existence of the unsaturated zone below. Pumping tests in the groundwater have been performed using production wells and a well specially implemented during this thesis work in order to estimate aquifer hydrodynamic parameters. Observations analysis and a conceptual modelling, in 2D and then in 3D, lead to a better understanding of the controlling mechanisms (stratification, input and output) and to simulate both basin infiltration rates and water table rise. Considering the whole basin scale, input fluxes are transient, related to the clogging statement of the infiltration area but also to the temperature of inflow water and the initial statement of the soil just below the basin. Sensibility analyses using the models highlight that the amount of the water table rise is mainly influenced by the aquifer saturated hydraulic conductivity and also by the location of imposed boundaries in the aquifer (rivers and pumping wells). The model properly accounts basin inflow fluxes and water table fluctuations. The model is able to verify if flows are reversed in relation to river exchanges, if wells are fed by infiltration basins and it highlights aquifer flows below the river. A 3D modelling has been realised in another area of the well field, comprising two infiltration basins, two river arms and pumping wells. In agreement with underflow in the aquifer, rivers are imposed in the model as third kind boundary conditions. Aquifer and river exchanges are calibrated with observed data of one aquifer flood-wave propagation. Significant differences of hydraulic heads have been observed at different depths of the aquifer using panpipes piezometers, specifically implemented, close to one infiltration basin. Theses differences are closely related to basin operation. These observations are properly calculated by the 3D model. Using the model, the effect of one sandy-clay heterogeneous layer (whose saturated hydraulic conductivity is lower than the ones of other aquifer lithologies) on aquifer flows (direction and flux) is notable. The model properly accounts basin inflow fluxes and water table fluctuations. The model is able to verify if flows are reversed in relation to river exchanges, if wells are fed by infiltration basins and it highlights aquifer flow below the river.

Devido aos riscos de contaminação dos recursos naturais solo e água, ao alto custo, ao tempo e ao esforço humano nas investigações de campo, os modelos matemáticos, aliados às técnicas numéricas e aos avanços computacionais, constituem uma ferramenta importante na previsão do deslocamento de solutos, contribuindo assim, para o controle de alterações ambientais. No Brasil, a modelação de fluxo e transporte de solutos na zona não-saturada é voltada, quase que exclusivamente, aos problemas relacionados às atividades agrícolas. Entretanto, tão importante quanto a problemática dos produtos químicos nas atividades agrícolas é a questão de poluição e contaminação do solo e da água por chorume, gerado pelos resíduos sólidos domiciliares. Neste trabalho, é desenvolvido e validado um modelo computacional unidimensional para simulação de fluxo e transporte de solutos na zona não-saturada do solo. O modelo matemático é dado pela equação diferencial parcial não-linear de Richards, que rege o movimento de água no solo, e a equação diferencial parcial linear de advecção-dispersão, do transporte de solutos, acompanhadas das condições iniciais e de contorno. A equação de Richards é dada em função do potencial matricial da água e a equação de transporte de solutos estima a evolução temporal da concentração de solutos no perfil do solo. Devido à dificuldade de se obter soluções analíticas destas equações, são resolvidas numericamente pelo método de elementos finitos. As referidas equações são resolvidas utilizando-se malhas uniformes inicialmente. Com a finalidade de obter simulações mais eficientes, a um custo computacional reduzido, é empregada a adaptatividade com refinamento h na malha de elementos finitos. A função interpolação polinomial utilizada é de grau 2 ou maior que garante a conservação de massa. Na equação de Richards, a derivada temporal é aproximada por um quociente de diferença finita e é aplicado o esquema de Euler explícito e na equação de advecção-dispersão, é aproximada por um quociente de diferença finita, aplicando-se o esquema de Euler implícito, devido à linearidade da equação. O sistema operacional é o Linux Ubuntu 32 bits, o ambiente de programação é o PZ, escrito em linguagem de programação C++. Na validação do modelo, utilizam-se dados disponíveis na literatura. Os resultados são comparados, utilizando-se malhas uniformes e malhas adaptativas com refinamento h. Usando-se as malhas uniformes para o problema de Richards e de transporte de potássio, o tempo de execução é de 22 minutos e a memória utilizada de 6164 Kb. Com as malhas adaptadas, o tempo de execução é de 3 minutos e 27 segundos, consumindo 5876 Kb de memória. Houve, portanto, uma redução de 84,32% no tempo de execução, usando-se malhas adaptativas. A utilização da função interpolação polinomial de grau 2 ou maior e o refinamento h, permitem uma boa concordância do modelo na comparação com soluções disponíveis na literatura. / Due to the risks of contamination of soil and water resources, the high cost, time and human effort in the field investigations, the mathematical models, combined with numerical techniques and computational advances, are important tools in forecasting the movement of solutes thereby contributing to the control of environmental alteration. In Brazil, modeling of flow and solute transport in the unsaturated zone is focused, almost exclusively, on problems related to agricultural activities. However, as important as the problematical of chemicals products in agricultural activities is the issue of pollution and contamination of soil and water by leachate, generated by municipal solid wastes. In this work, an one-dimensional computational model for simulation of flow and solute transport in the unsaturated soil has been developed and validated. The mathematical model is given by the Richards\'s non-linear partial differential equation, which determines the movement of water in the soil, and the advection-dispersion linear partial differential equation, of the solute transport, together with initial and boundary conditions. The Richards equation is a function of the water pressure head and the solute transport equation estimate the temporal evolution of the solutes concentration in the soil profile. Due to the difficulty of obtaining analytical solutions of these equations, they are solved numerically using the finite element method. The governing equations are solved using initially a uniform mesh. In order to obtain more efficient simulations with low computational cost, adaptativity with h refinement on the finite element mesh is implemented. The interpolation function is of degree two or higher, assuring mass conservation. In Richards\' equation, the temporal derivative is approximated by Euler explicit finite difference. For the advection-dispersion equation, due to the linearity of the equation, an implicit finite difference scheme is used. The code is written in the programming language C++ based on the programming environment PZ using operating system Linux Ubuntu 32 bit. Model results are validated in comparison with data available in the literature. The results are evaluated using uniform meshes and with h refinement adaptive mesh. Using the uniform meshes for the problem of Richards and transport of potassium, the running time is 22 minutes and 6164 Kb of memory is used. With the adapted meshes, the execution time is 3 minutes and 27 seconds, consuming 5,876 Kb of memory. Therefore there was a reduction of 84.32% in execution time, using adaptive meshes. The interpolation function with degree two or higher and the h refinement, with reduction of the computation time, showed a good agreement in comparison with the literature.