Theoretical and numerical tools for studying the Critical Zone from plot to catchments

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.