Implementing the Water, HEat and Transport model in GEOframe (WHETGEO1D v.1.0): algorithms, informatics, design patterns, open science features, and 1D deployment
 ^{1}Department of Civil, Environmental and Mechanical Engineering, University of Trento, Via Mesiano 77, 38123 Trento, Italy
 ^{2}Center Agriculture Food Environment, University of Trento, Via Mesiano 77, 38123 Trento, Italy
 ^{1}Department of Civil, Environmental and Mechanical Engineering, University of Trento, Via Mesiano 77, 38123 Trento, Italy
 ^{2}Center Agriculture Food Environment, University of Trento, Via Mesiano 77, 38123 Trento, Italy
Correspondence: Niccolò Tubini (niccolo.tubini@unitn.it)
Hide author detailsCorrespondence: Niccolò Tubini (niccolo.tubini@unitn.it)
This paper presents WHETGEO and its 1D deployment: a new physically based model simulating the water and energy budgets in a soil column. The purpose of this contribution is twofold. First, we discuss the mathematical and numerical issues involved in solving the Richardson–Richards equation, conventionally known as the Richards equation, and the heat equation in heterogeneous soils. In particular, for the Richardson–Richards equation (R^{2}) we take advantage of the nested Newton–Casulli–Zanolli (NCZ) algorithm that ensures the convergence of the numerical solution in any condition. Second, starting from numerical and modelling needs, we present the design of software that is intended to be the first building block of a new customizable landsurface model that is integrated with processbased hydrology. WHETGEO is developed as an opensource code, adopting the objectoriented paradigm and a generic programming approach in order to improve its usability and expandability. WHETGEO is fully integrated into the GEOframe/OMS3 system, allowing the use of the many ancillary tools it provides. Finally, the paper presents the 1D deployment of WHETGEO, WHETGEO1D, which has been tested against the available analytical solutions presented in the Appendix.
 Article
(8660 KB) 
Supplement
(547 KB)  BibTeX
 EndNote
The Earth's critical zone (CZ) is defined as the heterogeneous nearsurface environment in which complex interactions involving rock, soil, water, air, and living organisms regulate the natural habitat and determine the availability of lifesustaining resources (National Research Council, 2001). Clear interest in studying the CZ is spurred on by everincreasing pressure due to growth in the human population and climatic changes. Central to simulating the processes in the CZ is the study of soil moisture dynamics (Clark et al., 2015a; Tubini et al., 2021b). In the following we suggest that studying the CZ requires tools that are not yet readily available to researchers; then we propose one of our own. These tools should be flexible enough to allow the quick embedding of advancements in science.
1.1 Setting up the water budget
It is generally accepted in the field that the motion of water in soil is regulated by a Darcian–Stokesian flow, meaning that any force is immediately dissipated and water under a gradient of generalized forces acquires a velocity (the Darcy velocity) which is linearly proportional to the gradient of the generalized force. This law is known as the Darcy–Buckingham law and reads
where the forces acting are gravity (z; L) and the matric potential ψ [L] and where J [L T^{−1}] is the Darcian flux, K [L T^{−1}] is the hydraulic conductivity, θ [L^{3} L^{−3}] is the dimensionless volumetric water content, ∇ [L^{−1}] is the gradient operator, and z [L] is the vertical coordinate positiveupward. The assumptions under which such a law is derived from Newton's law are presented in Whitaker (1986) and Di Nucci (2014). The hydraulic conductivity depends on soil type (texture and structure) and water content, while the thermodynamic forces must be understood as gradients of the chemical potential of water, which, in turn, can be split in matric potential, osmotic potential, and other potentials (Nobel, 1999, chap. 2). However, in Eq. (1) we consider only the action of the matric potential. On the basis of the law of motion in Eq. (1), the mass conservation reads
where ∇⋅ [L^{−1}] is the divergence operator. Equation (2) is usually known as the Richards equation (Richards, 1931) but was previously formulated by Richardson (1922). Therefore, in the following we call it the R^{2} equation to remind the reader of this double origin. There are very informative reviews that cover its general, historical, and numerical aspects, such as Paniconi and Putti (2015), Farthing and Ogden (2017), and Zha et al. (2019). Therefore, it is not deemed necessary to further summarize the matter here. The R^{2} equation is a function of two variables, θ and ψ, and its resolution requires another relation between these two quantities. This relation is known as soil water retention curves (SWRCs), written in the plural because we have many SWRCs depending on soil characteristics. The reader may be aware that the R^{2} is an exact description of unsaturated flow only if we assume that soil is a bundle of capillaries and that the largest capillaries drain first and are filled last (Mualem, 1976). In fact, in this case a relation can be obtained between the radius of the capillaries and the suction, which was fully derived (Kosugi, 1999). However, there are various reasons to take the capillary bundle concept as a rough approximation of natural soil. Some of the issues are enumerated below.

Firstly, pores in soils are not bundles of welldefined capillaries of a single diameter; in fact, they can have quite random structures, as revealed, for instance, by tomography (Yang et al., 2018).

Secondly, logic and pore scale simulations, as in Tomin and Lunati (2016), for example, indicate that fluids fill the cavities where they fall, and only eventually are they redistributed according to the microstructure of the soil; that is to say that fluids do not move instantaneously from the largest pores to the smallest ones.

A set of relatively large pores can, in certain conditions, preferentially drive the flow of water in a short timescale according to laminar viscous flow driven by gravity before any redistribution happens (Germann and Beven, 1981).

The role of living matter, such as bacteria, animals, fungi, vegetation, and roots, is usually eliminated from the hydrological picture, but it should have a relevant place (Benard et al., 2019).
In addition are the following points.
 5.
Capillary forces are not the only ones acting at the microscale (Lu, 2016). In fact, measured suction values are far below pressures that can be sustained by capillarity alone.
 6.
Temperature affects water viscosity; infiltration is faster at warm temperatures and slower at cold ones (Constantz and Murphy, 1991).
 7.
In highlatitude and highelevation environments, soils may be subject to freezing and thawing processes which affect pore volume and water dynamics (Dall'Amico et al., 2011).
These facts certainly do not threaten the nature of mass conservation in Eq. (2). However, they can certainly alter the statistics which generate the closure equations, i.e. the SWRCs we currently use.

Requirement I – without entering into further details, we can observe that the aforementioned issues have consequences that would require new software to include the possibility of adopting new parameterizations of SWRCs and hydraulic conductivity quickly, easily, and neatly.
1.2 The three or four worlds
The flow of water obeys the general laws of physics for conservation of mass and momentum, but, since the seminal works of Freeze and Harlan (1969), the scientific community has split up (Furman, 2008) into three groups: groundwater people, vadose zone scientists, and surface water hydrologists. This compartmentalization of the scientific community was fostered to deepen the knowledge within single branches, with the interactions between the different parts having been governed in models by assigning boundary conditions (Furman, 2008). However, these boundary conditions are intrinsically inadequate and inappropriate in representing the physics of interactions between different domains whose interactions depend strictly on the state of the system. When these conditions are prescribed a priori (Furman, 2008), the proper dynamics of the CZ fluxes cannot be obtained. There is a need to overcome this situation.

Requirement II – the boundary conditions hardwired into algorithm implementation should be removed in favour of a simultaneous treatment of the three compartments (surface waters, vadose zone, and groundwater).
Fortunately, Šimůnek et al. (2012) found the way to smoothly extend the Richards equation into the groundwater equation. This and other similar approaches are now used in various codes, such as Hydrus, ParFlow (Ashby and Falgout, 1996; Jones and Woodward, 2001; Kollet and Maxwell, 2006), CATHY (Paniconi and Wood, 1993; Paniconi and Putti, 1994), and GEOtop 2.0 (Rigon et al., 2006a; Endrizzi et al., 2014). To extend the R^{2} equation into the saturated domain it is necessary to include the contribution of groundwater storativity due to matrix and fluid compressibility. The common approach is to write the R^{2} equation as
where S_{s} [L^{−1}] is the specific storage coefficient, defined as
with ρ [M L^{−3}] being the water density, g [L T^{−2}] gravitational acceleration, n [L^{3} L^{−3}] the soil porosity, β [L T^{2} M^{−1}] the liquid compressibility, and α [L T^{2} M^{−1}] the matrix compressibility. On the lefthand side of Eq. (3), the first term accounts for changes in liquid saturation, while the second term accounts for the compression or expansion of the porous medium and the water. The lefthandside term in Eq. (3) can be rewritten as
where c [L^{−1}] is the water retention capacity. Comparing the two terms in brackets, we can see that for ψ<0, then $c\gg {S}_{\mathrm{s}}\frac{\mathit{\theta}}{{\mathit{\theta}}_{\mathrm{s}}}$; this means that under unsaturated conditions, the contribution of the specific storage is negligible. Whereas when the soil is saturated and ψ>0, then c=0 and therefore what counts is the specific storage. Because of this, it is possible to account for groundwater specific storage simply by modifying the SWRC as
Furthermore, switching from Richards to shallow water was made possible in the equation writing thanks to, for example, Casulli (2017) and Gugole et al. (2018). Therefore, switching to a fully integrated, simultaneous treatment of the three domains is now possible.
1.3 The necessary coupling with the energy budget
As remarked in point 6 above, temperature affects water viscosity, which effectively doubles in passing from 5 to 20 ^{∘}C (Eisenberg et al., 2005), with a positive feedback on the infiltration process. This has been clearly observed in natural systems (Ronan et al., 1998; Eisenberg et al., 2005; Engeler et al., 2011) wherein infiltration rates follow diurnal and seasonal temperature cycles. In fact, according to Muskat and Meres (1936), the unsaturated hydraulic conductivity can be expressed as
where κ_{r}(θ) [−] is the relative permeability, κ [L^{2}] is the intrinsic permeability, ρ [L^{3} M^{−1}] is the liquid density, g is the acceleration of gravity, and ν [L^{2} T^{−1}] is the kinematic viscosity of the liquid. Thus, for constant θ, variations in K(θ) due to temperature can be accounted for as (Constantz and Murphy, 1991)
In Eq. (8), T_{1} is a reference temperature, while T_{2} is the soil water temperature. Temperature is also responsible for the phase change of water (point 7), and because of pore ice, as well as excess ice, infiltration rates and subsurface flows are significantly modified (Walvoord et al., 2012).

Requirement III – to account for thermodynamic effects, temperature should at least be present in the R^{2} equation as a parameter, as in Eq. (8). However, for a more accurate approximation of the water dynamics, the option to solve the water and energy budgets simultaneously must be present.
Soil thermal properties are important physical parameters in modelling landsurface processes (Dai et al., 2019) since they control the partitioning of energy at the soil surface and its redistribution within the soil (Ochsner et al., 2001). For a multiphase material like soil, their definition is always problematic since they depend on the physical properties of each phase and their variations (Dong et al., 2015; Dai et al., 2019; Nicolsky and Romanovsky, 2018). In the literature different models have been proposed with such a scope, and further studies on it are recommended (Dai et al., 2019): nonetheless, when considering the phase change of water, the estimation of the unfrozen and frozen water fraction is still an unresolved issue for which different models, usually referred to as SFCCs or soil freezing characteristics curves, have been proposed (Kurylyk and Watanabe, 2013). Thus, it is clear that the aspects related to the estimation of soil thermal properties fall fully within Requirement I too. Moreover, there are other reasons for which the equations of the water and energy budgets should be solved in a coupled manner (Rigon et al., 2006a): for instance, this makes it possible to include an appropriate treatment of evaporation and transpiration processes (Bonan, 2019; Bisht and Riley, 2019), as well as of the heat advection (Frampton et al., 2013; Walvoord and Kurylyk, 2016; Wierenga et al., 1970; Ronan et al., 1998; Engeler et al., 2011; Zhang et al., 2019).
Finally, there is a great urge to model solute transport according to the water movements. The range of applications for solute, tracers, and pollutants spans from agriculture to industry to research itself. In fact, in recent years there has been a tumultuous increase in the number of studies using tracers to assess the various pathways of water (Hrachowitz et al., 2016). However, so far these studies have mostly used lumped models with limited capability to investigate water age selection processes, with are processes that have become very important in the most recent literature, e.g. Penna et al. (2018). Using more complex modelling can benefit both the investigation of the processes and the construction of more refined water budget closures. Even though in this paper we do not detail the work on tracers, they must be kept in mind in software design so that the modules can be implemented eventually.
1.4 Heat transport
Under the conditions defined above, the governing equation for the transport of energy in variably saturated porous media is given by the following energy conservation equation:
where h is the specific enthalpy of the medium [L^{2} T^{−2}], λ [$\mathrm{M}\phantom{\rule{0.125em}{0ex}}\mathrm{L}\phantom{\rule{0.125em}{0ex}}{\mathrm{\Theta}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{T}}^{\mathrm{3}}$] is the thermal conductivity of the soil, T [Θ] is the temperature, ρ_{w} [M L^{−3}] is the water density, c_{w} [${\mathrm{L}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{\Theta}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{T}}^{\mathrm{2}}$] is the specific heat capacity of water, T_{ref} [Θ] is a reference temperature used to define the enthalpy, and J_{w} is the water flux [L T^{−1}]. The first term on the righthand side is the heat conduction flux, described by the Fourier's law, and the second term is the sensible heat advection of liquid water. The specific enthalpy of a control volume of soil V_{c} [L^{3}] can be calculated as the sum of the enthalpy of the soil particles and liquid water:
where ρ_{sp} and ρ_{w} are the densities of the soil particles and the water, and c_{sp} and c_{w} are the specific heat capacities of the soil particles and the water. Equation (9) is the socalled conservative form.
1.5 Informatics
Amateurism in software (model) building may be reflected in an incorrect analysis of the hydrological processes, threatening the development of science from the offset. Also, in recent years a new paradigm in evaluating scientific productivity, and therefore software productivity, has gained importance (see e.g. project IDEAS, 2021), with the aim of reducing the effort, time, and cost of software development, maintenance, and support. More in general, the scientific community has recognized the need to transition to open, accessible, reusable, and reproducible research, which requires the adoption of new software engineering practices. Our intention, therefore, is to build robust and reliable tools that are without unknown side effects and designed to be easily inspected by third parties.
As discussed in Serafin (2019), computing programmes have usually been developed as monolithic code; this has serious drawbacks for maintainability and development, as has been shown by our own experience with the GEOtop model (Rigon et al., 2006a; Endrizzi et al., 2014) and by experiences with other modelling frameworks (Clark et al., 2015b, 2021; Bisht and Riley, 2019). The main limit of the monolithic approach can be identified in its absence of separation of concerns (Martin, 2009; Newman, 2015). This results in huge, unreadable code bases that mix several different scientific and/or mathematical concepts (David et al., 2013), which are characteristics that impede easy reviews of models. Being aware of these past experiences, the WHETGEO1D code has therefore been developed by adopting objectoriented programming (OOP) and a generic programming strategy whose contents are explained below. In this manner it is possible to split up the code into smaller structures, with each one implementing a specific scientific and/or mathematical concept and also to decouple the algorithms from concrete data representation with gains in code maintainability and algorithm inspection.
Moreover, WHETGEO1D has been integrated into the Object Modelling System v3 (OMS3) framework (David et al., 2013); see Appendix B. OMS3 is a componentbased environmental modelling framework that allows the developers to create a component for each single modelling concept, thus having a second level of separation of concerns. Because of the modularity of OMS3, the WHETGEO1D components can be used seamlessly at runtime by connecting them with the OMS3 DSL language based on Groovy (https://groovylang.org, last access: 21 December 2021). OMS3 provides the basic services and, among these, tools for calibration and implicit parallelization of component runs.
WHETGEO is also part of GEOframe (Formetta et al., 2014a; Bancheri, 2017; Bancheri et al., 2020), a system of OMS3 components; see Appendix A. GEOframe provides various components for precipitation treatment, radiation estimation in complex terrain, and evaporation and transpiration. Conversely, with WHETGEO1D the GEOframe system is enriched with new components now able to cope with different spatial resolutions and configurations compared to what was offered by the existing tools (Bancheri et al., 2020).
1.6 Organization and scope
This paper describes the implementation and content of the WHETGEO1D (Water, HEat and Transport in GEOframe) software, in observance of Requirements I to III and aware of the hydrologic facts described in points 1 to 7 of Sect. 1.1. Actually, we do not seek to address all the issues listed but focus instead on getting the proper software design that their implementation requires. The first objective is to not have to break and rewrite existing codes and then be able to code incrementally, in step with the progressing research in hydrological processes.
Further requirements derive from numerics and mathematics. In fact, translating the R^{2} and the heat equations into discretized numerical codes implies overcoming some challenges that are well known in the hydrology community, as illustrated in the next sections. Our contribution aims to make available to the public for the first time integration algorithms that ensure the conservation of mass and energy for any time step size and for a great variety of conditions. It is interesting to note that some of these algorithms were presented almost a decade ago in the math community, remaining a little concealed, however. Although WHETGEO is meant to deal with the transport of solutes and tracers, this topic will be covered in a dedicated paper since its physics and numerics are slightly different from those of the water and energy budgets.
Moreover, WHETGEO1D represents the first brick of a new expandable landsurface model (Bisht and Riley, 2019): the forthcoming development of WHETGEO1D and its inclusion in the GEOSPACE (Soil Plant Atmosphere Continuum Estimator) model that aims to simulate the soil–plant continuum (D'Amato et al., 2022).
In the following section we discuss the mathematical issues and the numerics of our implementation. In the subsequent one, we describe the informatics and software engineering issues, as well as how we solved them. Finally, we discuss WHETGEO1D implementation by testing it against some analytical solutions available for simplified parameterizations and boundary conditions; for details of these see Appendix C.
Translating these equations into numerical discretized codes implies overcoming some challenges, as we illustrate in the next sections by exploring, at first, the issues posed by each of the equations.
2.1 General issues of the R^{2} equation
Equation (2) is said to be written in “mixed form” because it is expressed in terms of θ and ψ (and uses the SWRC to connect the two variables).
The “ψbased form” is derived from Eq. (2) by applying the chain rule for derivatives:
where
with dimension [L^{−1}] as the specific moisture capacity, also called hydraulic capacity. Even though Eqs. (2) and (11) are analytically equivalent under the assumption that the water content is a differentiable variable, this is not generally true in the discrete domain where the derivative chain rule is not always valid (Farthing and Ogden, 2017). Because of this the ψbased form may suffer from large balance errors in the presence of big nonlinearities and strong gradients, as discussed in Casulli and Zanolli (2010), Farthing and Ogden (2017), Zha et al. (2019), and literature therein. The specific moisture capacity c, which appears in the storage term, itself depends on ψ and so is not constant over a discrete time interval during which ψ changes value. Let us discretize the time derivative in Eq. (11) by using the backward Euler scheme and obtain
where ${\stackrel{\mathrm{\u0303}}{c}}_{i}$ is the discrete operator of the soil moisture capacity, c(ψ). In order to preserve the chain rule of derivatives at the discrete level, i.e. the equality $c\partial \mathit{\psi}/\partial t=\partial \mathit{\theta}\left(\mathit{\psi}\right)/\partial t$, ${\stackrel{\mathrm{\u0303}}{c}}_{i}$ has to satisfy the requirement (Roe, 1981)
As can be seen from the above equation, the right definition of ${\stackrel{\mathrm{\u0303}}{c}}_{i}$ depends on the solution itself. To overcome this problem, in the literature different techniques have been presented to improve the evaluation of ${\stackrel{\mathrm{\u0303}}{c}}_{i}$, but none ensures mass conservation (Farthing and Ogden, 2017).
There is a third form of Eq. (2), the socalled “θbased form”, that reads as
where all is made explicit by inverting the SWRC, and $D\left(\mathit{\theta}\right)=K\left(\mathit{\theta}\right)\partial \mathit{\psi}/\partial \mathit{\theta}$ is the soil water diffusivity [L^{2} T^{−1}]. The first term on the righthand side represents the water flow due to capillary forces, while the second term is the contribution due to gravity (Farthing and Ogden, 2017). The θbased form is massconserving, and it can be solved perfectly by mass conservative methods (Casulli and Zanolli, 2010). However, its applicability is limited to the unsaturated zone since water content varies between θ_{r} and θ_{s}, whereas water suction is not bounded. This formulation is intrinsically not suited to fulfilling our Requirement III. Moreover, water content is discontinuous across layered soil since the SWRCs are soilspecific, whereas water suction is continuous even in inhomogeneous soils (Farthing and Ogden, 2017; Bonan, 2019).
In WHETGEO, we directly use the conservative form of the R^{2} equation, Eq. (2), which seems the easiest way to deal with both the mass conservation issues and the extension of the equation to the saturated case.
2.1.1 The discretization of the R^{2} equation
The implicit finitevolume discretization of Eq. (2) reads as
where Δt is the time step size,
is an optional source–sink term in volume, and θ_{i}(ψ) is the ith water volume given by
Equation (16) can be written in matrix form as
where ψ={ψ_{i}} is the tuple of unknowns, θ(ψ)=θ_{i}(ψ_{i}) is a tuple function representing the discrete water volume, T is the flux matrix, and b is the righthandside vector of Eq. (16), which is properly augmented by the known Dirichlet boundary condition when necessary. For a given initial condition ${\mathit{\psi}}_{i}^{\mathrm{0}}$, at any time step $n=\mathrm{1},\mathrm{2},\mathrm{\dots}$, Eq. (16) constitutes a nonlinear system for ${\mathit{\psi}}_{i}^{n+\mathrm{1}}$, with the nonlinearity affecting only the diagonal of the system and being represented by the water volume ${\mathit{\theta}}_{i}\left({\mathit{\psi}}_{i}^{n+\mathrm{1}}\right)$. This set of equations is a consistent and conservative discretization of Eq. (2). Therefore, regardless of the chosen spatial and temporal resolution, ${\mathit{\psi}}_{i}^{n+\mathrm{1}}$ is a conservative approximation of the new water suction.
2.1.2 Surface boundary condition
The definition of the type of surface boundary condition (Neumann vs. Dirichlet) for the R^{2} equation is a nontrivial task since it can depend on the state of the system. In the literature several approaches are used (Furman, 2008). These approaches are mainly based on a switch of the type of the boundary condition from a prescribed head to prescribed flux and vice versa. This switching often causes numerical difficulties that need to be addressed (Furman, 2008).
To overcome this problem we have included an additional computational node at the soil surface. As will be made clear in the following, for it we prescribe an “equation state” like that presented in Casulli (2009):
where H [L] is the water depth, which also represents the pressure if the ponding is assumed to happen in hydrostatic conditions. By doing so, it is possible to prescribe as the surface boundary condition the rainfall intensity (Neumann type) without resorting to any switching techniques to reproduce infiltration excess or saturation excess processes. In this case the system in Eq. (19) must be modified to account for the additional computational node describing the state of the soil surface.
Here, J^{n} is the rainfall intensity and represents the Neumann boundary condition used to drive the system at the soil surface. For any time step, Eq. (21) can be written in matrix form, similar to (Eq. 19), as
where ψ={ψ_{i}} is the tuple of unknowns, V(ψ)=(θ_{i}(ψ_{i})) for $i=\mathrm{1},\mathrm{2},\mathrm{\dots},N\mathrm{1}$, V_{N}(ψ)=H(ψ) is a tuple function representing the discrete water volume, T is the flux matrix, and b is the righthandside vector of Eq. (16), which is properly augmented by the known Dirichlet boundary condition when necessary. For a given initial condition ${\mathit{\psi}}_{i}^{\mathrm{0}}$, at any time step $n=\mathrm{1},\mathrm{2},\mathrm{\dots}$ Eq. (22) constitutes a nonlinear system for ${\mathit{\psi}}_{i}^{n+\mathrm{1}}$, with the nonlinearity affecting only the diagonal of the system and being represented by the water volume ${V}_{i}\left({\mathit{\psi}}_{i}^{n+\mathrm{1}}\right)$. Therefore, regardless of the chosen spatial and temporal resolution, ${\mathit{\psi}}_{i}^{n+\mathrm{1}}$ is a conservative approximation of the new water suction.
2.2 Heat transport numerics issues
Equation (9) is said to be written in conservative form and expresses an important property, which is the conservation of the scalar quantity, in this case the specific enthalpy. It is interesting to note that by making use of the mass conservation equation (Eq. 2), Eq. (9) can be written in an analytically equivalent form, the socalled nonconservative form (Sophocleous, 1979; Šimůnek et al., 2005):
Equation (23) expresses another important property, which is the maximum principle (Casulli and Zanolli, 2005); i.e. the analytical solution is always bounded, above and below, by the maximum and minimum of its initial and boundary values, as shown in Greenspan and Casulli (1988, chap. 7.3).
Although Eqs. (9) and (23) are analytically equivalent, once they are discretized the corresponding numerical solution will, in general, either be conservative or satisfy a discrete max–min property (Casulli and Zanolli, 2005), but not both as would be required.
As in the case of the water flow, thermal budget issues can be subdivided into three aspects: the discretization of the equation, the inclusion of the appropriate boundary conditions, and the implementation of some closure equation for the thermal capacity and conductivity.
2.2.1 The discretization of the heat equation
The key feature (Casulli and Zanolli, 2005) to obtaining a numerical solution for the heat transport equation that is both conservative and possesses the max–min property is to solve the conservative form of the heat equation by using the velocity field obtained in solving the continuity equation (Eq. 2). By making use of the upwind scheme for the advection part and the centred difference scheme for the diffusion part we have
where
When the heat equation does not consider water phase changes, it is decoupled from the R^{2} equation and the finitevolume discretization leads to a linear algebraic system of equations. However, once freezing and thawing processes are considered, the heat equation is fully coupled with the R^{2} equation, as in Dall'Amico et al. (2011) for instance, and the enthalpy function becomes nonlinear. At this point, since the enthalpy function is nonlinear the NCZ algorithm is required to linearize it, as shown in Tubini et al. (2021b). So far, we have not considered the problem of water flow in freezing soils; however, being aware of this issue is important for the future developments and code design.
2.2.2 Driving the heat equation with the surface energy budget
At the soil surface the heat equation is driven by the surface energy balance. The heat flux exchanged between the soil and the atmosphere, the surface heat flux, ℱ [M T^{−3}], is given as
where S_{in} is the incoming shortwave radiation, S_{out} is the outgoing shortwave radiation, L_{in} is the incoming longwave radiation, L_{out} is the outgoing longwave radiation, and H and LE are respectively the turbulent fluxes of sensible heat and latent heat. Fluxes are positive when directed toward the soil surface and all have the dimension of an energy per unit area per unit time [M T^{−3}].
Similarly to the definition of the surface boundary condition for the water flow, the surface boundary condition for the energy equation is also systemdependent. In fact, in Eq. (26) the only fluxes that do not depend on the soil temperature and/or moisture are the incoming shortwave and longwave radiation fluxes, S_{in} and L_{in}. The outgoing shortwave radiation flux is usually parameterized as
where the surface albedo α [−] can be assumed to vary with the soil moisture content (Saito et al., 2006) and radiation wavelength. The outgoing longwave surface radiation is
where T_{s} [Θ] the temperature of the topmost layer of soil, ϵ is the soil emissivity, and σ is the Stefan–Boltzmann constant. The sensible heat flux H is taken as
where ρ_{a} is the air density [M L^{−3}], and c_{a} is the thermal capacity of air per unit mass [${\mathrm{L}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{T}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{\Theta}}^{\mathrm{1}}$]. Regarding the aerodynamic resistances r_{H} [T L^{−1}], it should be noted that it can be evaluated with different degrees of approximation and may require a specific modelling solution. For instance, the aerodynamic resistance r_{H} can be evaluated with models ranging from semiempirical models to the Monin–Obukhov similarity (Liu et al., 2007) or even by solving the turbulent dynamics with direct methods (Raupach and Thom, 1981; Mcdonough, 2004).
The latent heat flux is taken here as given by a formula of the type
where l [M L^{2} T^{−2}] is the specific latent heat of vaporization of water, E_{P} is the potential evapotranspiration, and r_{H} and r_{v} [T L^{−1}] are respectively the aerodynamic resistance and the soil surface resistance to water vapour flow. The latent heat flux it is the sum of two distinct processes: evaporation and transpiration. Compared to the other fluxes, latent heat flux presents further complications because evaporation is both an energy and a waterlimited process, and transpiration also depends on the physiology of trees (as well as root distribution and growth and leaf cover). The latent heat flux is associated with the water flux that must be accounted for in the R^{2} equation. Here we present a simplified treatment of the latent heat flux as an external driving force and/or prescribed boundary condition. A more exhaustive and physically based treatment of the latent heat flux, and the related water flux, is addressed in the ongoing development of the GEOSPACE model (D'Amato, 2021).
Including the surface energy budget boundary condition requires the computation of additional quantities such as the incoming radiation fluxes, the shortwave radiation and the longwave radiation, and the potential evapotranspiration flux. These quantities can be easily computed within the GEOframe system in which WHETGEO1D is embedded. The proper estimation of the incoming radiation fluxes is far from being a simple task, and it is often oversimplified in hydrological problems. Our approach is to use the tools already developed inside the system GEOframe, which were tested independently and accurately (Formetta et al., 2013, 2014a). Similarly, the evapotranspiration can be computed with other GEOframe components (Bottazzi, 2020; Bottazzi et al., 2021).
2.3 Algorithms
By using a numerical method, here the finitevolume method, a partial differential equation is transformed into a system of nonlinear algebraic equations, as has already been shown. The system has to be solved with iterative methods, and, at their core, these reduce the problem to using a linear system solver. The solver can be of various types, according to the dimension of the problem. For instance, in 1D, the final system of finitevolume problems we present is tridiagonal and can be conveniently solved with the Thomas algorithm (Quarteroni et al., 2010), which is a fast direct method. In 2D or 3D, the final matrix is not tridiagonal, and a different solution method must be used: for instance, the conjugate gradient (Shewchuk, 1994). These algorithms are well known and do not need to be explained here.
However, the reduction of a nonlinear system to a linear one is not trivial. We illustrate the issues by taking the R^{2} equation as an example. As discussed in depth in Zha et al. (2019), Farthing and Ogden (2017), and references therein, the linearization of the R^{2} equation is challenging. Following the work of Celia et al. (1990), a lot of advancements have been made in this direction: Hydrus, CATHY, and ParFlow use variants of the Newton and Picard iteration methods (Zha et al., 2019; Paniconi and Putti, 1994), while GEOtop 2.0 implements a suitable globally convergent Newton method (Kelley, 2003). Although current algorithms are relatively stable, they may fail to converge or require a considerable computational cost (Zha et al., 2019). This has significant impacts on both the reliability of the solution, which can have mass balance errors, and the computational cost to produce it (Farthing and Ogden, 2017; Zha et al., 2019). Since Casulli and Zanolli (2010) and Brugnano and Casulli (2008), a new method was found, called nested Newton by the authors and NCZ in the following, that guarantees convergence in any situation, even with the use of large time steps and grid sizes.
As clearly pointed out by Casulli and Zanolli (2010), what makes the linearization of the R^{2} equation difficult is the nonmonotonic behaviour of the soil moisture capacity. A mathematical proof of convergence for NCZ exists (Brugnano and Casulli, 2008, 2009; Casulli and Zanolli, 2010, 2012), which is not repeated here. However, we take the time to illustrate this new algorithm with care.
Let us start again from the nonlinear system (Casulli and Zanolli, 2012):
where ψ=(ψ_{i}) is the tuple of unknowns, V(ψ)=(V_{i}(ψ_{i})) is a nonnegative vectorial function, and the V_{i}(ψ_{i}) terms are defined for all ψ_{i}∈ℝ and can be expressed as
For all $i=\mathrm{1},\mathrm{2},\mathrm{\dots},N$, the following assumptions are made for the functions a_{i}(ψ) (we are quite literally following Casulli and Zanolli, 2010 here):
 A1:
a_{i}(ψ) is defined for all ψ∈ℝ and is a nonnegative function with bounded variations;
 A2:
${\mathit{\psi}}_{i}^{*}\in \mathbb{R}$ exists such that a_{i}(ψ) is strictly positive and nondecreasing in $(\mathrm{\infty},{\mathit{\psi}}_{i}^{*})$ and nonincreasing in $({\mathit{\psi}}_{i}^{*},+\mathrm{\infty})$.
T in Eq. (31) is the socalled matrix flux, and it is a symmetric and (at least) positive semidefinite matrix satisfying one of the following properties:
 T1:
T is a Stieltjes matrix, i.e. a symmetric M matrix, or
 T2:
T is irreducible, null(T)≡span(v), with v>0 (componentwise), and T+D is a Stieltjes matrix for all diagonal matrices D O, with O denoting the null matrix.
Finally, b is the vector of the known terms. When T satisfies property T2, then for Eq. (31) to be physically and mathematically compatible, the following assumption about b is required:
where ${\mathit{V}}^{\text{Max}}={\int}_{\mathrm{\infty}}^{+\mathrm{\infty}}{a}_{i}\left(\mathit{\xi}\right)\mathrm{d}\mathit{\xi}$.
Having assumed that the a_{i}(ψ) terms are nonnegative functions of bounded variations, they are differentiable almost everywhere, admit only discontinuities of the first kind, and can be expressed as the difference of two nonnegative, bounded, and nondecreasing functions, say p_{i}(ψ) and q_{i}(ψ), such that
for all ψ∈ℝ. When a(ψ) terms satisfy assumptions A1 and A2, the corresponding decomposition (known as the Jordan decomposition as in Chistyakov, 1997, and presented in Fig. 2) is given by
where ${\mathit{\psi}}_{i}^{*}$ is the position of the maximum of p_{i}. Thereafter, V(ψ) can be expressed as
where the ith component of V_{1}(ψ) and V_{2}(ψ) is defined as
By making use of Eq. (36) the algebraic system in Eq. (31) can be written as
It is necessary here to point out exactly how the nonlinear system in Eq. (31) reads when considering only the R^{2} equation and when the water depth function is used to properly define the surface boundary condition. In the first case, i.e. when Neumann or Dirichlet boundary conditions are used, the vectorial function is defined as V(ψ)=(θ_{i}(ψ_{i})) for $i=\mathrm{1},\mathrm{2},\mathrm{\dots},N$.
Instead, when we consider the water depth function to describe the computational node at the soil surface, the vectorial function is defined as V(ψ)=(θ_{i}(ψ_{i})) for $i=\mathrm{1},\mathrm{2},\mathrm{\dots},N\mathrm{1}$ and V_{N}(ψ)=H(ψ). Therefore, the nonlinear system in Eq. (38) is valid to describe both the subsurface and surface waters when the symbols are appropriately understood.
This aspect, the use of two different equation states, and the fact that the NCZ algorithm can be successfully reused to solve other problems (Casulli and Zanolli, 2012; Tubini et al., 2021b), requires a careful design of its implementation, as discussed in the following sections.
The concepts and requirements previously illustrated must be cast into software whose usability, expandability, and inspectability are demanded by good software design, which adds further requirements. The software design requirements and the software deployment are explained below.
3.1 Design requirements
One of the major difficulties encountered by a research group concerns the development and reuse of scientific software (Berti, 2000) and the writing of structurally clean code, i.e. a code that is easily readable and understandable, with objects that have a specified and possibly unique responsibility (Martin, 2009).
An objectoriented programming approach, with the adoption of standard design patterns (DPs) (Gamma et al., 1995; Freeman et al., 2004) and the creation of new ones, has been adopted for the internal class design and hierarchy.
The design principles followed by the WHETGEO1D software can be summarized as follows.
 a.
The software should be opensource to allow inspection and improvements by third parties.
 b.
For the same reason it should be organized into parts, each with a clear functional meaning and possibly a single responsibility.
 c.
The software can be extended with minimal effort and without modifications (according to the “open to extensions, closed to modifications” principle; Freeman et al., 2004). In particular, the parts to be modified are those that, according to the discussion in the previous sections, could be changed to try new closures, i.e. the SWRCs and the hydraulic conductivities in the case of the R^{2} equation and the thermal capacity and thermal conductivity in the case of the energy budget. Adding a new SWRC type or a new conductivity function should be easy.
 d.
The largest set of boundary conditions should be smoothly manageable.
 e.
The implementation of equations should be abstract, according to the principle of “programming to interfaces and not to concrete classes”, which is the core of contemporary OOP (Gamma et al., 1995). The different equations describing the processes should be implemented within the set of classes by implementing a common interface.
 f.
The implementation of algorithms should not depend on the data formats of inputs and outputs.
Another requirement has to do with the user experience. In fact, solvers of partial differential equations (PDEs) (Menard et al., 2020) tend to be complex to understand and run when features are added. In particular, the number of inputs grows exponentially when features are added, and the user has to overcome a steep learning curve before being able to use these software packages to appreciate all the cases implemented and their physics.
 g.
To simplify this situation, WHETGEO1D has to be implemented is such a way that any of the alternative implementations must come only with their own parameters and variables, as well as appear to the user to be as simple as possible, though not too simple.
There is finally a last requirement to consider.
 h.
For computational and research purposes, there will be one, two, and threedimensional (1D, 2D, 3D) implementations of the aforementioned equations. Therefore, as much of the code as possible should be shared across these. In particular, the NCZ and Newton algorithms should be shareable across the various applications.
This requirement implies that the geometry of the domain, as well as the topology, should be specified in an abstract manner to cope with the specifics of each dimensionality.
The rest of this section is organized to respond to points (a) to (h). Point (a) is actually responded to in the next section describing where the software can be downloaded and with which open license. Points (b) and (c) are accomplished by studying an appropriate design of classes and the use of design patterns (Gamma et al., 1995; Freeman et al., 2004). For (d), generic programming is used and specific classes are implemented. Point (e) is resolved by deploying a set of classes that implement a common interface with an extension that allows it to obtain the required functionalities.
To respond to the issues raised in (f) and (g), WHETGEO1D is implemented as various components of OMS3, as shown below, each one with its own inputs. Therefore, the number of flags to check and the number of unused inputs are reduced to the minimum required by the solvers and parameterizations that the user chooses. If users want to solve the R^{2} equation alone, for instance, they can pick out the appropriate component and do not need to know about the inputs and details of the energy budget. The separation into components has two other advantages. First, it eases the testing of a single process against available analytical solutions and against other model results (Bisht and Riley, 2019). Second, it improves the model structure, facilitating the representation of new processes (Clark et al., 2021).
Point (h) is solved by deploying new components for the 1D, 2D, and 3D cases. In the following section we mainly deal with points (b) to (e). Before discussing details of some classes, a few general choices have to be reported. Data of any type are stored internally in vectors of doubles, which in turn encapsulated in appropriate Java objects. OOP good practice would suggest that an object should be immutable (Bloch, 2001), but we decided that the main classes have to be mutable and allocated once forever as singletons (Gamma et al., 1995; Freeman et al., 2004). This potentially exposes the software to side effects but frees it to allocate new objects at any time step and decreases the computational burden and memory occupancy generated by deallocating unused obsolete objects at runtime. This approach may be considered a specific design pattern for partial differential equation solvers.
3.2 The software organization
The more visible effect of our choices is that we have built various OMS3 components.

whetgeo1d1.0beta

netcdf1.0beta

closureequation1.0beta

buffer1.0beta

numerical1.0beta
Internally, the classes are assembled by using some interfaces and abstract classes, since WHETGEO1D is coded using the Java language.
In order to improve the reusability of the Java code we adopted a generic programming approach (Berti, 2000) that consists of decoupling of algorithm implementations from the concrete data representation while preserving efficiency. The generic approach has been balanced with domainspecific ones that can improve the computational efficiency of the software, as is the case of the previously mentioned Thomas algorithm used in 1D implementations.
Another requirement regards the division of software classes into three main groups, as the lack of a proper separation between the parameterization of physical processes and their numerical solutions has been recognized as one of the weak points of existing landsurface models (Clark et al., 2015b, 2021). One group describes the mathematical–physical problem, the second one implements the numerical solution (Berti, 2000), and the third one contains the growing group of concrete classes. The first group contains the SWRC, as well as hydraulic and thermal conductivity, and it forms a standalone library since its content can eventually be reused in the 2D and 3D version of WHETGEO1D. Similarly, we grouped the classes that solve linear and nonlinear algebraic systems, containing the Thomas conjugate gradient and the various Netwon types of algorithms, into a second standalone library. The third group of classes gathers the concrete implementations and the variety of OMS3 components that are deployed.
The classes used, and their repository for thirdparty inspection, are illustrated in the 00_Notebooks referred to in Sect. 4 and in the Supplement. However, there are three pivotal groups of classes that we want to mention here: these contain the description of the geometry of the integration domain, the closure equation, and the state equation.
3.2.1 Computational domain, i.e. the Geometry
class
One of the key aspects to have in a generic solver regards the management of the grid and, in particular, the definition of its topology, or how grid elements are connected to each other. In the 1D case the description of the topology is quite simple since it can be implicitly contained in the vector representation: each element of the vector corresponds to a control volume of the grid, and it is only connected with the elements preceding and following it. It is worth noting that this approach is peculiar to 1D problems and cannot be adopted for the 2D and 3D domains, where, especially when unstructured grids are used, the grid topology requires a smart implementation of the incidence and adjacency matrices.
For each control volume it is necessary to store its geometrical quantities, their position and dimension, its variables, its parameter set, and the form of the equation to be solved there, referred to in the following as “equation state”. The appropriate arrangement of information, together with the internal design of the classes, allows us to create a generic finitevolume solver.
3.2.2 Closure equations, i.e. the ClosureEquation
abstract class
The ClosureEquation
class is shown in the Unified Modeling Language (UML) diagram in
Fig. 3. As explained in the Introduction, one of the core
concepts of modelling water and heat transport in soils is the SWRC. Soil is a
multiphase material, and thus knowledge of its composition is of crucial
importance in defining its unsaturated hydraulic conductivity as well as its thermal
properties, specific internal energy, and thermal conductivity.
An abstract class, ClosureEquation
, is defined to contain only abstract
methods that would be overwritten by the concrete classes implementing it. The
ClosureEquation
class essentially defines a new data type. A closer
inspection of Fig. 3 reveals that the
ClosureEquation
is composed by aggregation with the
Parameter
class, which contains all the physical parameters of the
model. Moreover, the Parameter
class is implemented by using the
Singleton pattern (Freeman et al., 2004). This design pattern is functional
to have only one instance of the Parameter
class shared by all the
ClosureEquation
objects and to provide a global point of access to
the Parameter
instance.
The Simple Factory pattern SoilWaterRetentionCurvesFactory
accomplishes the task of implementing the concrete classes. By preferring
polymorphism to inheritance and using the Simple Factory pattern
(Gamma et al., 1995; Freeman et al., 2004), the developers can easily include
and extend existing code or new formulations or parameterizations of
SWRC. Also, the Simple Factory fulfils the dependency inversion principle
(Eckel, 2003), and thus new extensions cannot affect the functioning
of existing code. The same closure equation, for instance a particular SWRC,
can be used to compute the soil water volume when solving the R^{2} equation
and the specific enthalpy of the soil when solving the heat equation.
3.2.3 Equation state, i.e. the EquationState
class
The EquationState
class in Fig. 4 contains
the implementation of the discretized form of the equation state of the PDE
under scrutiny. It contains a reference to the ClosureEquation
object, the Geometry
object, and the ProblemVariables
object.
Notably, the solver of a PDE problem can refer to the abstract class and to
its abstract object without implying the specific concrete equation to be
solved or its concrete parameterizations. Moreover, the compositionality of
the EquationState
allows the creation of new solvers from existing
closures without the need to add new subclasses. As shown in the UML in
Fig. 4, the EquationState
class defines methods used to linearize the
PDE when it is nonlinear. For instance, it computes the first and second
derivatives and the functions necessary to define the Jordan decomposition as
required by the NCZ algorithm. Specifically, these methods are p
,
q
, pIntegral
, qIntegral
, and computeXStar
.
In our code design the ClosureEquation
class is limited to computing
a physical parameterization, whereas the EquationState
class is used
to discretize the equation state of the PDE and whenever required to properly
linearize it. Any new concrete EquationState
subclass can either have
the same physics of another with a different solver or different physics
with the same solver.
3.3 Generic programming at work
As explained in Sect. 2.1.2, the definition of the surface boundary condition for the R^{2} equation can require the introduction of an additional computation node at the soil surface to simulate the water depth. This means that we have two different equation states: one for the soil water content and one for the water depth. On the other hand, when considering the Neumann or Dirichlet boundary condition we have only one equation state for the soil water content. The nonlinear solver of the NCZ algorithm must work seamlessly with any type of boundary condition used and with any type of equation involved in the problem.
Let us consider, for example, the case of the R^{2} equation with water depth
as the surface boundary condition. Figure 5 reports
the pseudocode for computing the nonlinear functions in the NCZ algorithm
using a procedural approach. When traversing the computational grid we need to
resort to an ifelse
statement to compute the nonlinear function of
each control volume: the control volume N represents the control volume for
the surface water, and in this case we need to use the water depth equation
state, H(ψ) in Fig. 5. The remaining control
volumes represent the soil moisture, and for them we need to use the SWRC
equation state, θ(ψ) in Fig. 5. The main
limit of this approach is that the computation of the nonlinear functions,
H(ψ) and θ(ψ), is hardwired into the code of the NCZ
algorithm. This presents a shortcoming for the reusability of the code since,
as the boundary condition changes and the physical problem changes, it is
necessary to modify the for
loop and the name of the objects
computing the nonlinear functions.
Adopting the OOP and generic programming approach
(Fig. 6), it is possible to implement the NCZ algorithm
in such a way that enhances its reusability. The key feature is the decoupling
of the computational grid from the algorithm (data)
(Berti, 2000). This is achieved through two elements. The first
consists of creating a container of the objects that deal with the equation
states of the problem, equationState
, eS in
Fig. 6. Second, we use a label
equationStateID
, eSID in Fig. 6, to specify
the behaviour of each control volume. So, the behaviour of each control volume
is determined by this label and not by the position of the element in the
grid. Specifically, when we traverse the grid we use the
equationStateID
to determine which object inside the container
equationState
to use.
The NCZ algorithm has been implemented in the NestedNewtonThomas
class. The NestedNewtonThomas
contains a reference to the
Thomas
object, whose task is to solve a linear system, and to a list
of EquationState
objects (Fig. 7). Considering the
ubiquity of nonlinear problems in hydrology and the robustness of the NCZ
algorithm, the NCZ algorithm has been encapsulated in a standalone library.
The example has been illustrated in 1D, but it becomes even more effective when working in 2D or 3D, especially with an unstructured grid.
While most of what has been written so far is of general application, the deployment shown here is 1D. Information on WHETGEO1D for users and developers is provided in the Supplement, where there is a Jupyter Notebook that contains the guidelines for executing the codes for any of the components. Its name starts with “00_”, and we call it “Notebook Zero” of the components. The latest executable code can be downloaded from

https://github.com/geoframecomponents/WHETGEO1D (last access: 21 December 2021)
and can be compiled by following the instructions therein. The version of the OMS3 compiled project can be found here: https://github.com/GEOframeOMSProjects/OMS_Project_WHETGEO1D (last access: 21 December 2021). The code can be executed in the OMS3 console, which can be downloaded and installed according to the instructions given at

https://geoframe.blogspot.com/2020/01/thewinterschoolongeoframesystemis.html (last access: 21 December 2021).
Some brief information about GEOframe can be found in Appendix B, and more comprehensive information is at

https://abouthydrology.blogspot.com/2015/03/jgrassnewageessentials.html (last access: 21 December 2021) and

https://geoframe.blogspot.com/2015/03/ (last access: 21 December 2021).
To run the tests, please follow the instructions on the GitHub repository of the GEOframe components. If a user wants to compile the code themselves, they can use the appropriate Gradle script that guarantees independence from any IDE.
For further information about input and output formats for WHETGEO1D, please see the Notebook 00_WHETGEO1D_Richards.ipynb in the folder documentation of the Zenodo distribution.
4.1 Workflow for users
Examples of uses of WHETGEO1D can be found in the form of Python Notebooks in the directory Notebooks/Jupyter. Documentation can be found in the form of Python Notebooks in the directory documentation. Simulations with WHETGEO1D are run as OMS3 simulations. Therefore, the first operation to accomplish is to prepare the appropriate .sim files. For new users, many simulation files are available in the directory simulation of the Zenodo distribution. As shown in Fig. 8, in the modelling solutions that involve WHETGEO1D, there is always a “main” component that is in charge of running the core code for solving the PDE. The inputs and the outputs are treated by other OMS3 components. They are tied together by a domainspecific language (DSL) based on Groovy. This allows for great flexibility in using various input and output formats.
4.2 Inputs and outputs
Input data can be broadly classified into time series, computational grid data, and simulation parameters. Time series are used to specify the boundary conditions of the problem. Time series are contained in .csv files with a specific format that is OMS3compliant (David, 2021). With computational grid data we refer to the domain discretization, initial condition, and soil parameters. All these data are stored in a NetCDF file. Time series and computational grid data are elaborated with dedicated Python modules distributed under the geoframepy package (Tubini and Rigon, 2021a). The simulation parameters, such as the start date and end date of the simulation, time step size, and file paths, are specified by the user in the OMS3 .sim file.
For the design of the output workflow we took advantage of the OMS3 system
that allows the user to connect standalone
components. Figure 8 shows the output workflow for
saving output data with Main
, Buffer
, and
netCDF writer
as the standalone OMS3 components. Main
stands for the generic component having the responsibility of solving the
PDE. Buffer
has the responsibility of temporarily storing output
data, and Writer
handles the saving of data to the disk. The
Buffer
component has the sole purpose of storing data, and this has
two important advantages. The first is that it limits the number of accesses
to the disk to save output, i.e. reducing the computational time. The second
is that it introduces a layer separating the Main
component from the
netCDF writer
. This increases the flexibility of the modelling
solution, as future developers can adopt different file formats or develop
different writer components that, instead of saving all the outputs, can save
discrete outputs or aggregated outputs. The advantage is that developers need
only know the legacy of the Buffer
component and customize both the
output file format and memory optimization strategy, such as chuncking,
according to their need. Currently all outputs are stored in a NetCDF3 format
(Unidata, 2021). NetCDF is a selfdescribing portable data format developed
and maintained by UCAR Unidata. NetCDF is commonly used by the geoscience
community, and there is an evergrowing number of tools for processing and
visualization.
The choice of including the Buffer
component in the workflow of the
modelling solution is motivated by earlier experiences of the GEOtop community
with the GEOtop model and by more recent experience with the FreeThaw1D model
(Tubini et al., 2021b), according to which long spinups of the model (∼1500 years)
were required with consequently large output files (∼ GB). Furthermore, in
anticipation of the 2D and 3D developments, the NetCDF3 format is probably
not the most appropriate and could be abandoned in favour of more performing
file formats (Unidata, 2021a, b), such as NetCDF4 or HDF5.
WHETGEO1D can be integrated with the builtin calibration component
LUCA
(Hay et al., 2006; Formetta et al., 2014a) and the
Verification
component, as shown in
Fig. 9. The former is used to calibrate
optimal parameters and the latter to compute the indices of goodness of simulated
data versus measured data. Besides the LUCA
component and the
Verification
component, it is necessary to add two more components,
specifically the Buffer calibration parameters
and the
Measurement points data
. The Buffer
calibration parameters are
needed to interface the Main
component with the LUCA
component. In fact, in the WHETGEO1D Main
component, physical
parameters are stored as vectors, whilst LUCA
handles calibration
parameters as scalars (single value). The Buffer calibration
component receives the optimal parameter set from the LUCA
component
and returns them packed in appropriate vectors. The Verification
component receives as input two OMS3compliant time series: one for measured
data and one for simulated data. In this case it is necessary to extract from
the simulation output only the simulated data at the measurement points of the
variable used to calibrate the model. These data are then saved as OMS3 time
series. It interesting that the integration of WHETGEO1D with the OMS3
builtin calibration components is achieved by adding two new standalone
components without modifying the source code of the existing components,
i.e. the Main
component, the Buffer
component, and the
netCDF writer
.
4.3 Workflow for developers
Here, as an example, we present how to add the Brooks–Corey (Brooks and Corey, 1964) model as an extension of the code base. The constitutive relationships are given by
where θ_{s} and θ_{r} are the saturated and residual values of the volumetric water content, respectively, ψ_{d} is the airentry water suction value, n is the pore size distribution index, and K_{s} is the saturated hydraulic conductivity at saturation.
The standard approach to adding a new SWRC parameterization, here the Brooks
and Corey model, requires the definition of a new class that extends the abstract
class ClosureEquation
. This new class, SWRCBrooksCorey
,
provides the implementation of the abstract methods defined in the super class
ClosureEquation
and inherits the association with the Parameters
class. Specifically, the SWRCBrooksCorey
class overrides the
following methods:

f
calculates the water content for a given water suction value and set of parameters (Eq. 39); 
df
calculates the first derivative of Eq. (39); 
ddf
calculates the second derivative of Eq. (39).
In order to use the Brooks–Corey model in the Richards equation it is
necessary to define a new class, SoilWaterVolumeBrooksCorey
, that
extends the abstract class StateEquation
. Specifically, the
SoilWaterVolumeBrooksCorey
class overrides the following methods.

equationState
calculates the water volume using the Brook–Corey model for a given water suction value and set of parameters (Eq. 18). 
dEquationState
calculates the first derivative of theequationState
function, in case this the moisture capacity function. This method is used within the linearization algorithm. 
ddEquationState
calculates the second derivative of theequationState
function. This method is relevant for models for which the ψ^{*} cannot be computed analytically but requires the application of a rootfinding method such as the bisection method. An example is the soil internal energy function when considering the phase change of water (Tubini et al., 2021b). 
p
calculates the p function of the Jordan decomposition. 
pIntegral
calculates the V_{1} function of the Jordan decomposition. 
computeXStar
calculates the ψ^{*} value to properly define the functions p and V_{1}. 
initialGuess
calculates the initial guess for the linearization algorithm.
In this paper we discussed the issues raised by implementing a new expandable system to model the Earth's CZ, WHETGEO, the core of which is a reliable and robust integrator of the Richards–Richardson equation and the associated energy budget.
We would like to stress the following.

WHETGEO makes available the NCZ algorithm, which has a priori convergence ensured for any choice of time step and for a great variety of boundary and initial conditions.

The coupling between the R^{2} equation and the heat advection diffusion equation is done using a recent algorithm (Casulli and Zanolli, 2005) that, to our knowledge, has not been applied before for the same scope.

The application of our methods could benefit a large range of Earth science models, which at present require complicated workarounds (Regenass et al., 2021).

The adoption of the mixed form of the equation allows the joint modelling of the saturated and unsaturated soil and groundwater conditions.

Moreover, the coupling strategy with ponding water makes WHETGEO naturally suitable to investigate the partitioning of infiltration and surface runoff without any of the issues documented for other software (e.g. Šimůnek et al., 2005).
The implementation of WHETGEO has been shown to solve three software requirements and the (a) to (h) design specifications. Each of these was analysed and thus influenced the choice of algorithms and code implementation. Of the hydrological issues presented in seven observations, the issues 1 to 4 were also solved, and the code provided some answers for observations 5 to 7 because of the following.

The code allows for the easy changing of soil water retention curves; therefore, for instance, it would be easy to implement new ones taking more accurate account of the different processes in the various ranges of suction. A description of how to add a new soil water retention curve is presented in Sect. 4.3.

The effect of temperature on viscosity has been included, and, in the Supplement, we have shown that temperature changes in soil have relevant effects on infiltration and runoff production.

Freezing and thawing processes coupled with the R^{2} equation will be inserted. As an intermediate step, the heat conduction with phase change has been implemented in the energy budget solver.
The flexibility of the code was obtained with the creation of the 𝚂𝚝𝚊𝚝𝚎𝙴𝚚𝚞𝚊𝚝𝚒𝚘𝚗 and the 𝙲𝚕𝚘𝚜𝚞𝚛𝚎𝙴𝚚𝚞𝚊𝚝𝚒𝚘𝚗 classes, as illustrated in Sect. 3. Their organization can be considered a new programming pattern that can be reused by programmers to implement the solution of partial differential equations independently of our use of the Java language. Furthermore, we also chose to implement the objects containing data by adopting the Singleton pattern. This choice violates the safety prescription that objects should preferably be immutable, but, by avoiding the continuous allocation and deallocation (through the Java garbage collector) of new objects, we aimed at saving computational time. Also, this scheme can be considered a design choice that can be replicated by anyone interested in scientific computing.
The first deployment of the concepts was the 1D standalone water budget as well as the coupled water and energy budget WHETGEO versions. The water budget was tested against analytical solutions presented in Srivastava and Yeh (1991) and Vanderborght et al. (2005). Some behavioural simulations were also performed to show some features of the code, such as the ability to deal with switching boundary conditions. Also, we showed that WHETGEO can be easily extended to simulate the thermal regime in frozen soils, as described in Tubini et al. (2021b), by merely adding the SFC model presented by Dall'Amico et al. (2011).
The model is an entirely new code that was accurately designed to allow easy inspection and to support code literacy and class reuse. It is opensource and built with opensource tools. It and its documentation fulfil the requirements of open science (Hall et al., 2021). In a science for which a lot of daily research activities are based on computer simulations, this is a requirement that most existing models do not fulfil but that the authors believe is profitable for the progress of science. Also, WHETGEO is built with a chain of opensource tools, making it available to all researchers who want to peruse it without constraint. Its documentation was produced using only opensource tools and is available together with the distribution of the main code.
WHETGEO1D was implemented as a Java component within the GEOframe, an opensource, componentbased hydrological modelling system. Within GEOframe, each part of the hydrological cycle is implemented in a selfcontained building block, an OMS3 component (David et al., 2013). Components can be joined together to obtain multiple modelling solutions that can accomplish simple to very complicated tasks. GEOframe has shown great flexibility and robustness in several applications (Bancheri et al., 2020; Abera et al., 2017a, b). There are more than 50 components available that can be grouped into the following categories:

geomorphic and DEM analyses;

spatial extrapolation or interpolation of meteorological variables;

estimation of the radiation budget;

estimation of evapotranspiration;

estimation of runoff production with integral distributed models;

channel routing;

travel time analysis;

calibration algorithms.
Using the components for geomorphic and DEM analyses (Rigon et al., 2006b), the basin can be discretized into hydrological response units (HRUs), i.e. hydrologically similar parts, such as a catchment or a hillslope or one of its parts. The meteorological forcing data can be spatially interpolated using a geostatistical approach, such as the kriging technique (Bancheri et al., 2018). Both shortwave and longwave radiation components are available for the estimation of the radiation budget (Formetta et al., 2013, 2016). Evapotranspiration can be estimated using three different formulations: the FAO evapotranspiration model (Allen et al., 1998), the Priestley–Taylor model (Priestley and Taylor, 1972), and the Prospero model (Bottazzi, 2020; Bottazzi et al., 2021). Snow melting and the snow water equivalent can also be simulated with three models, as described in Formetta et al. (2014b). Runoff production is performed by using the Embedded Reservoir Model (ERM) or a combination of its reservoirs (Bancheri et al., 2020). The discharge generated at each hillslope is routed to the outlet using the Muskingum–Cunge method (Bancheri et al., 2020). Travel time analysis of a generic pollutant within the catchment can be done using the approach proposed in Rigon et al. (2016b, a). Model parameters can be calibrated using two algorithms and several objective functions: Let us calibrate (LUCA) (Hay et al., 2006) and particle swarm optimization (PSO) (Kennedy and Eberhart, 1995). A graphbased structure, called NET3 (Serafin, 2019), is employed for the management of process simulations. NET3 is designed using a river network and graph structure analogy, whereby each HRU is a node of the graph, and the channel links are the connections between the nodes. In any NET3 node, a different modelling solution can be implemented and nodes (HRUs or channels) can be connected or disconnected at runtime through scripting. GEOframe is opensource and helps the reproducibility and replicability of research (Bancheri, 2017). Developers and users can easily collaborate, share documentation, and archive examples and data within the GEOframe community.
The Object Modelling System v.3 (OMS3) is a componentbased environmental modelling framework that provides a consistent and efficient way to (1) create science simulation components; (2) develop, parameterize, and evaluate environmental models, as well as modify and adjust them as science advances; and (3) repurpose environmental models for emerging customer requirements (David et al., 2013).
In OMS3 the term component refers to selfcontained, separate software units that implement independent functions in a contextindependent manner (David et al., 2013). This means that developers and researchers can build their model as a composition of standalone components, moving away from the monolithic approach. The entire GEOframe system, and therefore WHETGEO1D, is built upon the OMS3 framework.
Compared to other environmental modelling frameworks (EMFs), OMS3 is characterized by being a noninvasive and lightweight framework (Lloyd et al., 2011). That is to say that the model code is not tightly coupled with the underlying framework – OMS3; i.e. the environmental modeller does not need a deep knowledge of the API, and the modelling components can still function and continue to evolve outside the framework (David et al., 2013). In fact, OMS3 relies on specific annotations to provide metadata for Java code. These annotations describe elements such as classes, fields, and methods and are used by the framework to interpret the component as a building block of the modelling solution (MS), hence controlling its connectivity and data flow (David et al., 2013). It is worth noting that, being metadata, these annotations do not directly affect the execution of the source code outside the OMS3 noninvasive and lightweight framework.
Besides the technical aspects, the adoption of a software framework has a positive effect on “nonfunctional” quality attributes, such as maintainability, portability, reusability, and understandability (David et al., 2013). The componentbased approach allows the developer to break down the problem into smaller parts, each one tackled by a specific component. Hence, the components are joined together to build the desired modelling solution (point b). This facilitates the construction of new MSs, thanks to the plugin system of model components (David et al., 2013; Peckham et al., 2013; Serafin, 2019). Thanks to the modularity, the updating of a component with the most recent scientific advances is facilitated and has no side effects on the other components. The other advantage regards the longterm development of the code. From past experiences, one of the main limits to model development and maintenance was related to the lack of a proper software architectural design (Rizzoli et al., 2006; David et al., 2013; Formetta et al., 2014a; Bancheri, 2017; Serafin, 2019). Moreover, it is interesting to note that the componentbased approach encourages collective model development (Serafin, 2019) and also eases the attribution of authorship since any component is a standalone chunk of code and can be authored separately.
Also, the adoption of an environmental modelling framework promotes the concept of reproducible research, easing thirdparty inspection and providing consistent and verifiable model results (Formetta et al., 2013; Bancheri, 2017; Serafin, 2019).
Another advantage of using OMS3 is represented by the opportunity to keep the code development transparent to the user.
In this section we test the solver of the R^{2} equation against the analytical solutions presented by Srivastava and Yeh (1991) and by Vanderborght et al. (2005). Then we present and discuss two “behavioural” test cases to try out WHETGEO1D in simulating both the infiltration excess and the saturation excess process.
C1 Analytical solution Srivastava and Yeh (1991)
Srivastava and Yeh (1991) derived an analytical solution describing the onedimensional transient infiltration in a homogeneous and layered soil. The hydraulic properties of the soil are described by the following constitutive relations:
where K_{s} is the saturated hydraulic conductivity, θ_{r} is the residual water content, θ_{s} is the saturated water content, and α is the soil pore size distribution parameter representing the desaturation rate of the SWRC. The lower boundary condition is represented by the water table, ψ=0, while the upper boundary condition is subjected to a constant flux, q. The initial condition corresponds to the steadystate profile due to a prescribed initial flux at the soil surface and prescribed pressure at the lower boundary. The analytical solution is derived by linearizing the Richards equation and using Laplace transformation. Details on the analytical solution can be found in Srivastava and Yeh (1991).
C1.1 Homogeneous soil
We consider a onedimensional homogeneous soil layer of 1 m depth (TP1). The saturated hydraulic conductivity value is assumed to be 1.0 cm h^{−1}, with θ_{s}=0.45 m^{3} m^{−3}, θ_{r}=0.2 m^{3} m^{−3}, and α=0.01 cm^{−1}. The initial condition is determined by imposing as a lower boundary condition ψ=0 m and a constant water flux at the soil surface of q_{A}=0.1 cm h^{−1}. For times greater than 0 the water flux at the soil surface is q_{B}=0.9 cm h^{−1}. The domain is discretized with a uniform grid space Δz=0.001 m and the time step is Δt=60 s. The model accuracy is enhanced by allowing two Picard iterations per time step. Figure C1 shows a comparison between the numerical and the analytical solutions.
C1.2 Layered soil
In this numerical problem (TP2) we consider onedimensional vertical infiltration toward the water table through a layered soil. The initial condition is determined by imposing as a lower boundary condition ψ=0 m and a constant water flux at the soil surface of q_{A}=0.1 cm h^{−1}. For times greater than 0 the water flux at the soil surface is q_{B}=0.9 cm h^{−1}. The domain is discretized with a uniform grid space Δz=0.001 m and the time step is Δt=60 s. The model accuracy is enhanced by allowing two Picard iterations per time step. The hydraulic conductivity at the interface is computed as the harmonic mean of the neighbours (Romano et al., 1998). Comparison between the numerical and the analytical solution for water suction is shown in Fig. C3.
C2 Analytical solution of Vanderborght et al. (2005)
The next test case was defined by Vanderborght et al. (2005) to evaluate the steadystate flux in layered soil profiles. For this numerical problem (TP3) we consider a soil column of 2 m depth with one soil type for depth 0–0.5 m overlying another soil type for depth 0.5–2 m, specifically for loam over sand, sand over loam, and clay over sand. The soil parameters are defined in Table C1. The initial condition for water suction is a uniform profile with $\mathit{\psi}=\mathrm{20}$ m, the surface boundary condition is a constant flux of $q=\mathrm{5.79}\times {\mathrm{10}}^{\mathrm{8}}$ m s^{−1}, and at the bottom we impose a free drainage boundary condition. The domain is discretized with a uniform grid space Δz=0.01 m and the time step is Δt=3600 s. In order to reach the steadystate condition the simulation lasts 2 years. Comparison between the numerical and the analytical solution is shown in Fig. C5.
C3 Surface boundary condition
The definition of the surface boundary condition is a nontrivial task since it is a systemdependent boundary condition. The infiltration rate through the soil surface depends on precipitation, rainfall intensity J, and on the moisture condition of the soil. Because of this, the surface boundary condition may change from the Dirichlet type – prescribed water suction – to the Neumann type – prescribed flux – and vice versa. The works by Horton (1933) and Dunne and Black (1970) establish the conceptual framework to explain the runoff generation.
The infiltration excess or Horton runoff occurs when the rainfall intensity is larger than the infiltration capacity of the soil:
Infiltration excess is most commonly observed with shortduration, intense rainfall.
The saturation excess or Dunnian runoff occurs when the soil is saturated and additional water exfiltrates at the soil surface. Saturation excess generally occurs with longduration, moderate rainfall or with a series of successive precipitation events. In this case the soil depth or the presence of shallow fragipan is the determining factor for saturation excess. Another possible cause is the rise of the water table up to the soil surface.
C3.1 Infiltration excess
In this numerical experiment we consider a homogeneous soil of 3 m depth. Soil hydraulic properties are described with the Van Genuchten model (Table C2).
The initial condition is assumed to be hydrostatic with ψ=0 m at the bottom. The surface boundary condition is synthetic rainfall, as in Fig. C6a, lasting 15 min with constant intensity of 0.028 mm s^{−1}. At the bottom we prescribed a Dirichlet boundary condition with constant ψ=0 m so the transient is driven only by the surface boundary condition. In Fig. C6a, the time is indicated when it would be necessary to switch from the Neumanntype to the Dirichlettype boundary condition.
Figure C7 shows a comparison of water ponding at the soil surface considering two different initial conditions of the soil: wet and dry. For the wet case, the initial condition is hydrostatic with ψ=0 m at the bottom. For the dry case, the initial condition is hydrostatic with $\mathit{\psi}=\mathrm{100}$ m at the bottom. In the wet initial condition the hydraulic conductivity is higher than for the dry initial condition; however, in the dry case the capillary gradient is larger and because of this the soil infiltration capacity is higher, as in Fig. C7a. With regards to the water ponding, the maximum value is almost the same in both cases, 1 mm higher in the wet case, but the time evolution is different: in the wet case the water only infiltrates completely 13 h later than the dry case. This delay may seem counterintuitive since wetter conditions are associated with higher values of hydraulic conductivity (Fig. C8), but in the wet soil the capillary gradients are smaller than in the dry soil (Fig. C9).
C3.2 Saturation excess
In this section we present two numerical experiments to simulate the saturation excess process. Saturation excess is more critical, in terms of simulation stability, than the infiltration excess (PC Progress Discussion Forums, 2021). We consider two cases: one in which the water table reaches the soil surface and another in which the total rainfall amount is larger than the maximum waterholding capacity but the rainfall intensity is less than the maximum infiltration rate.
Firstly, we consider a layered soil of 3 m depth. The thicknesses of the loamy layer and clay layer are 0.5 and 2.5 m, respectively. The soil hydraulic properties are described with the Van Genuchten model (Table C3). The initial condition is assumed to be hydrostatic with ψ=0 m at the bottom. At the surface boundary condition we prescribe no rainfall, while at the bottom a variable Dirichlet boundary condition is prescribed (Fig. C10a). The transient is driven by variation of the water table. In Fig. C10a the time is indicated when it would be necessary to switch the surface boundary condition from the Neumanntype to the Dirichlettype boundary condition and vice versa.
Secondly, we consider a layered soil of 3 m depth. The thicknesses of the loamy sand layer and the clay layer are 0.3 and 2.7 m, respectively. The soil hydraulic properties are described with the Van Genuchten model (Table C4). The initial condition is assumed to be hydrostatic with $\mathit{\psi}=\mathrm{2}$ m at the bottom. The surface boundary condition is synthetic rainfall (Fig. C11a), and at the bottom we prescribed a Dirichlet boundary condition with constant $\mathit{\psi}=\mathrm{2}$ m so the transient is driven only by the surface boundary condition. Figure C11b shows the time evolution of the degree of saturation within the soil. Initially water infiltrates into the soil but then the clay layer, which is characterized by a lower conductivity than the loam–sand layer, limits the deep infiltration causing the saturation of the loam–sand layer from below.
Repeating the above numerical experiment with a thicker loam–sandy layer (Fig. C12), there is no water ponding at the soil surface. In this case all the rainfall can infiltrate into the loam–sandy layer thanks to the increased water storage capacity.
Here we present a behavioural test case on the pure heat conduction considering the surface energy balance. The purpose of this simulation is twofold. First, we show an application of WHETGEO1D that exploits existing GEOframe components to model the external components of the surface energy budget, specifically the incoming shortwave radiation (Formetta et al., 2013), the incoming longwave radiation (Formetta et al., 2016), and the latent heat flux. Second, we show how to easily include the phase change of water by adding a new closure equation that describes the SFC model presented by Dall'Amico et al. (2011)
The soil column is 30 m deep, and the initial condition is a constant temperature profile T=12 ^{∘}C. Figure D1a shows the components of the surface energy fluxes and the thermal regime of the uppermost 2 m of the soil column. As can be seen in panel (b), the soil temperature falls below 0 ^{∘}C; therefore, it is not reasonable to neglect freezing and thawing processes.
By including the closure equation of the SFC presented by Dall'Amico et al. (2011), it is possible to consider the phase change of water. Moreover, compared to the FreeThaw1D model (Tubini et al., 2021b), in WHETGEO1D it is possible to drive the simulation of the soil thermal regime by using the surface energy budget. This aspect represents a novelty with respect to FreeThaw1D and was obtained with minimum effort thanks to the code design we adopted. The soil column is of 30 m depth, and the initial condition is a constant temperature profile T=12 ^{∘}C. The parameters of the SFC are presented in Table D1. Figure D1 shows in panel (a) the components of the surface energy fluxes and the thermal regime of the uppermost 2 m of the soil column. As can be seen in panel (b), the soil temperature falls below 0 ^{∘}C, and thus neglecting freezing and thawing processes results in a strong approximation.
Figure D3 shows a comparison between the two simulations.
The source code is written in Java using the objectoriented programming paradigm. The source code can be found at https://github.com/geoframecomponents/WHETGEO1D (last access: 21 December 2021) and a frozen version at https://doi.org/10.5281/zenodo.5112727 (Tubini and Rigon, 2021c). The OMS3 project can be found at https://github.com/GEOframeOMSProjects/OMS_Project_WHETGEO1D (last access: 21 December 2021), and a frozen version is available on Zenodo at https://doi.org/10.5281/zenodo.4749319 (Tubini and Rigon, 2021b). The source code of external dependencies is provided in the README of the GitHub page at https://github.com/geoframecomponents/WHETGEO1D (last access: 21 December 2021). WHETGEO1D is deployed as an opensource code to work alone or within the Object Modelling System version 3 framework (David et al., 2013, https://doi.org/10.1016/j.envsoft.2012.03.006). In the latter case it can be connected at runtime with the many other components developed in the GEOframe system to provide hydrometeorological forcings and other fluxes, such as evapotranspiration. The code must be run within the OMS3 console or using the Dockerized version of OMS3. To set up the environment please follow the steps described in the README file present in the GitHub repository at https://github.com/GEOframeOMSProjects/OMS_Project_WHETGEO1D (last access: 21 December 2021). and in the GEOframe pages at https://geoframe.blogspot.com/2021/05/whetgeo1d.html (last access: 18 May 2021, GEOframe, 2021b). Once you have installed OMS3, please follow the instructions contained in the documentation folder. They contain all the details about simulation inputs and parameters. For those interested in using WHETGEO1D and applying it in their research, it is available as part of the material for the first GEOframe Summer School held in Trento on 4–7 October 2021 at http://geoframe.blogspot.com/2021/06/geoframesummerschool2021gss2021.html (last access: 21 December 2021, GEOframe Summer School, 2021). The material of the school is available here: https://doi.org/10.17605/OSF.IO/Z53C6 (Tubini et al., 2021a). The fair use and publication policy of the GEOframe group is available here: https://osf.io/wgdyq/ (last access: 21 December 2021, GEOframe, 2021a).
The simulations and data presented here can be found at https://doi.org/10.5281/zenodo.4749319 (Tubini and Rigon, 2021b).
The supplement related to this article is available online at: https://doi.org/10.5194/gmd15752022supplement.
RR and NT were responsible for conceptualization. NT and RR developed the methodology. NT and RR were responsible for software engineering. NT performed code writing. RR and NT performed code revision. NT was responsible for data curation and simulations. NT and RR contributed to paper writing, review, and editing. NT and RR were responsible for documentation. RR acquired funding.
The contact author has declared that neither they nor their coauthor has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The authors would like to thank Professor Vincenzo Casulli and Professor Michael Dumbser of the Department of Civil, Environmental and Mechanics Engineering at the University of Trento for their fruitful discussions on the numerical aspects of the work. We thank Concetta D'Amato for the help in preparing the supporting material. We would like to thank Joseph Eamonn Tomasi for his invaluable help with the linguistic revision. We especially thank the two anonymous reviewers for their constructive comments that helped to improve our paper.
This work has been partially supported by a PhD grant by the Department of Civil, Environmental and Mechanics Engineering at the University of Trento and by the Italian MIUR project (PRIN 2017) “WATer mixing in the critical ZONe: observations and predictions under environmental changesWATZON” (project code: 2017SL7ABC).
This paper was edited by MinHui Lo and reviewed by three anonymous referees.
Abera, W., Formetta, G., Borga, M., and Rigon, R.: Estimating the water budget components and their variability in a prealpine basin with JGrassNewAGE, Adv. Water Resour., 104, 37–54, https://doi.org/10.1016/j.advwatres.2017.03.010, 2017a. a
Abera, W., Formetta, G., Brocca, L., and Rigon, R.: Modeling the water budget of the Upper Blue Nile basin using the JGrassNewAge model system and satellite data, Hydrol. Earth Syst. Sci., 21, 3145–3165, https://doi.org/10.5194/hess2131452017, 2017b. a
Allen, R. G., Pereira, L. S., Raes, D., and Smith, M.: Crop evapotranspirationGuidelines for computing crop water requirements, Fao, Rome, FAO Irrigation and drainage paper 56, 300(9), D05109, 1998. a
Ashby, S. F. and Falgout, R. D.: A parallel multigrid preconditioned conjugate gradient algorithm for groundwater flow simulations, Nucl. Sci. Eng., 124, 145–159, https://doi.org/10.13182/NSE96A24230, 1996. a
Bancheri, M.: A flexible approach to the estimation of water budgets and its connection to the travel time theory, PhD thesis, University of Trento, 2017. a, b, c, d
Bancheri, M., Serafin, F., Bottazzi, M., Abera, W., Formetta, G., and Rigon, R.: The design, deployment, and testing of kriging models in GEOframe with SIK0.9.8, Geosci. Model Dev., 11, 2189–2207, https://doi.org/10.5194/gmd1121892018, 2018. a
Bancheri, M., Rigon, R., and Manfreda, S.: The GEOframeNewAge Modelling System Applied in a Data Scarce Environment, Water, 12, 86, https://doi.org/10.3390/w12010086, 2020. a, b, c, d, e
Benard, P., Zarebanadkouki, M., and Carminati, A.: Physics and hydraulics of the rhizosphere network, J. Plant Nutr. Soil Sc., 182, 5–8, https://doi.org/10.1002/jpln.201800042, 2019. a
Berti, G.: Generic software components for Scientific Computing, PhD thesis, Von der Fakultät für Mathematik, Naturwissenschaften und Informatik der Brandenburgischen Technischen Universität Cottbus, 2000. a, b, c, d
Bisht, G. and Riley, W. J.: Development and verification of a numerical library for solving global terrestrial multiphysics problems, J. Adv. Model. Earth Sy., 11, 1516–1542, https://doi.org/10.1029/2018MS001560, 2019. a, b, c, d
Bloch, J.: Effective Java: Programming Language Guide, AddisonWesley Professional, ISBN: 9780201310054, 2001. a
Bonan, G.: Climate change and terrestrial ecosystem modeling, Cambridge University Press, https://doi.org/10.1017/9781107339217, 2019. a, b, c, d, e, f
Bottazzi, M.: Transpiration theory and the Prospero component of GEOframe, PhD thesis, University of Trento, 2020. a, b
Bottazzi, M., Bancheri, M., Mobilia, M., Bertoldi, G., Longobardi, A., and Rigon, R.: Comparing evapotranspiration estimates from the GEOframeProspero model and three empirical models under different climate conditions, Water, 13, 1221–1243, https://doi.org/10.3390/w13091221, 2021. a, b
Brooks, R. H. and Corey, T. A.: Hydraulic properties of porous media Hydrology Paper No. 3, Civil Engineering Department, Colorado State University, Fort Collins, CO, 1964. a
Brugnano, L. and Casulli, V.: Iterative solution of piecewise linear systems, SIAM J. Sci. Comput., 30, 463–472, https://doi.org/10.1137/070681867, 2008. a, b
Brugnano, L. and Casulli, V.: Iterative solution of piecewise linear systems and applications to flows in porous media, SIAM J. Sci. Comput., 31, 1858–1873, https://doi.org/10.1137/08072749X, 2009. a
Casulli, V.: A highresolution wetting and drying algorithm for freesurface hydrodynamics, Int. J. Numer. Meth. Fl., 60, 391–408, https://doi.org/10.1002/fld.1896, 2009. a
Casulli, V.: A coupled surfacesubsurface model for hydrostatic flows under saturated and variably saturated conditions, Int. J. Numer. Meth. Fl., 85, 449–464, https://doi.org/10.1002/fld.4389, 2017. a
Casulli, V. and Zanolli, P.: High resolution methods for multidimensional advection–diffusion problems in freesurface hydrodynamics, Ocean Model., 10, 137–151, https://doi.org/10.1016/j.ocemod.2004.06.007, 2005. a, b, c, d
Casulli, V. and Zanolli, P.: A nested Newtontype algorithm for finite volume methods solving Richards' equation in mixed form, SIAM J. Sci. Comput., 32, 2255–2273, https://doi.org/10.1137/100786320, 2010. a, b, c, d, e, f
Casulli, V. and Zanolli, P.: Iterative solutions of mildly nonlinear systems, J. Comput. Appl. Math., 236, 3937–3947, https://doi.org/10.1016/j.cam.2012.02.042, 2012. a, b, c
Celia, M. A., Bouloutas, E. T., and Zarba, R. L.: A general massconservative numerical solution for the unsaturated flow equation, Water Resour. Res., 26, 1483–1496, https://doi.org/10.1029/WR026i007p01483, 1990. a
Chistyakov, V.: On mappings of bounded variation, J. Dyn. Control Syst., 3, 261, https://doi.org/10.1007/BF02465896, 1997. a
Clark, M. P., Fan, Y., Lawrence, D. M., Adam, J. C., Bolster, D., Gochis, D. J., Hooper, R. P., Kumar, M., Leung, L. R., Mackay, D. S., Maxwell, R. M., Shen, C., Swenson, S. C., and Zeng, X.: Improving the representation of hydrologic processes in Earth System Models, Water Resour. Res., 51, 5929–5956, https://doi.org/10.1002/2015WR017096, 2015a. a
Clark, M. P., Nijssen, B., Lundquist, J. D., Kavetski, D., Rupp, D. E., Woods, R. A., Freer, J. E., Gutmann, E. D., Wood, A. W., Brekke, L. D., Arnold, J. R., Gochis, D. J., and Rasmussen, R. M.: A unified approach for processbased hydrologic modeling: 1. Modeling concept, Water Resour. Res., 51, 2498–2514, https://doi.org/10.1002/2015WR017198, 2015b. a, b
Clark, M. P., Zolfaghari, R., Green, K. R., Trim, S., Knoben, W. J. M., Bennet, A., and Nijssen, B.: The numerical implementation of land models: problem formulation and laugh tests, J. Hydrometeorol., 22, 1627–1648, https://doi.org/10.1175/JHMD200175.1, 2021. a, b, c
Constantz, J. and Murphy, F.: The temperature dependence of ponded infiltration under isothermal conditions, J. Hydrol., 122, 119–128, https://doi.org/10.1016/00221694(91)90175H, 1991. a, b
Dai, Y., Wei, N., Yuan, H., Zhang, S., Shangguan, W., Liu, S., Lu, X., and Xin, Y.: Evaluation of soil thermal conductivity schemes for use in land surface modeling, J. Adv. Model. Earth Sy., 11, 3454–3473, https://doi.org/10.1029/2019MS001723, 2019. a, b, c
Dall'Amico, M., Endrizzi, S., Gruber, S., and Rigon, R.: A robust and energyconserving model of freezing variablysaturated soil, The Cryosphere, 5, 469–484, https://doi.org/10.5194/tc54692011, 2011. a, b, c, d, e
D'Amato, C.: BrokerGEO, GitHub [code], available at: https://github.com/geoframecomponents/BrokerGEO, last access: 20 April 2021. a
D'Amato, C., Tubini, N., Benettin, P., Rinaldo, A., Noto, V., and Rigon, R.: Investigation of a critical zone column with the GEOframe Plant Atmosphere Continuum Estimator (GEOSPACE), Hydrol. Earth Syst. Sci., in preparation, 2022. a
David, O.: CSV Data Files, codeBeamer [data set], available at: https://alm.engr.colostate.edu/cb/wiki/16970, last access: 1 April 2021. a
David, O., Ascough II, J. C., Lloyd, W., Green, T. R., Rojas, K., Leavesley, G. H., and Ahuja, L. R.: A software engineering perspective on environmental modeling framework design: The Object Modeling System, Environ. Modell. Softw., 39, 201–213, https://doi.org/10.1016/j.envsoft.2012.03.006, 2013. a, b, c, d, e, f, g, h, i, j, k
Di Nucci, C.: Theoretical derivation of the conservation equations for single phase flow in porous media: a continuum approach, Meccanica, 49, 2829–2838, https://doi.org/10.1007/s110120140022y, 2014. a
Dong, Y., McCartney, J. S., and Lu, N.: Critical review of thermal conductivity models for unsaturated soils, Geotech. Geol. Eng., 33, 207–221, https://doi.org/10.1007/s1070601598432, 2015. a
Dunne, T. and Black, R. D.: An experimental investigation of runoff production in permeable soils, Water Resour. Res., 6, 478–490, https://doi.org/10.1029/WR006i002p00478, 1970. a
Eckel, B.: Thinking in JAVA, Prentice Hall Professional, 2003. a
Eisenberg, D., Kauzmann, W., and Kauzmann, W.: The structure and properties of water, Oxford University Press, ISBN: 9780198570264, 2005. a, b
Endrizzi, S., Gruber, S., Dall'Amico, M., and Rigon, R.: GEOtop 2.0: simulating the combined energy and water balance at and below the land surface accounting for soil freezing, snow cover and terrain effects, Geosci. Model Dev., 7, 2831–2857, https://doi.org/10.5194/gmd728312014, 2014. a, b
Engeler, I., Franssen, H. H., Müller, R., and Stauffer, F.: The importance of coupled modelling of variably saturated groundwater flowheat transport for assessing river–aquifer interactions, J. Hydrol., 397, 295–305, https://doi.org/10.1016/j.jhydrol.2010.12.007, 2011. a, b
Farthing, M. W. and Ogden, F. L.: Numerical solution of Richards' equation: A review of advances and challenges, Soil Sci. Soc. Am. J., 81, 1257–1269, https://doi.org/10.2136/sssaj2017.02.0058, 2017. a, b, c, d, e, f, g, h
Formetta, G., Rigon, R., Chávez, J. L., and David, O.: Modeling shortwave solar radiation using the JGrassNewAge system, Geosci. Model Dev., 6, 915–928, https://doi.org/10.5194/gmd69152013, 2013. a, b, c, d
Formetta, G., Antonello, A., Franceschi, S., David, O., and Rigon, R.: Hydrological modelling with components: A GISbased opensource framework, Environ. Modell. Softw., 55, 190–200, https://doi.org/10.1016/j.envsoft.2014.01.019, 2014a. a, b, c, d
Formetta, G., Kampf, S. K., David, O., and Rigon, R.: Snow water equivalent modeling components in NewAgeJGrass, Geosci. Model Dev., 7, 725–736, https://doi.org/10.5194/gmd77252014, 2014b. a
Formetta, G., Bancheri, M., David, O., and Rigon, R.: Performance of sitespecific parameterizations of longwave radiation, Hydrol. Earth Syst. Sci., 20, 4641–4654, https://doi.org/10.5194/hess2046412016, 2016. a, b
Frampton, A., Painter, S. L., and Destouni, G.: Permafrost degradation and subsurfaceflow changes caused by surface warming trends, Hydrogeol. J., 21, 271–280, https://doi.org/10.1007/s100400120938z, 2013. a
Freeman, E., Robson, E., Bates, B., and Sierra, K.: Head first design patterns, O'Reilly Media, Inc., ISBN 9780596007126, 2004. a, b, c, d, e, f
Freeze, R. A. and Harlan, R. L.: Blueprint for a physicallybased, digitallysimulated hydrologic response model, J. Hydrol., 9, 237–258, 1969. a
Furman, A.: Modeling coupled surface–subsurface flow processes: A review, Vadose Zone J., 7, 741–756, https://doi.org/10.2136/vzj2007.0065, 2008. a, b, c, d, e
Gamma, E., Helm, R., Johnson, R., Vlissides, J., and Patterns, D.: Elements of Reusable ObjectOriented Software, Design Patterns, AddisonWesley Publishing Company, Massachusetts, 1995. a, b, c, d, e
GEOframe: FairUse&Publication Policy_v1.pdf, OSF, available at: https://osf.io/wgdyq/, last access: 21 December 2021a. a
GEOframe: WHETGEO1D, available at: https://geoframe.blogspot.com/2021/05/whetgeo1d.html, last access: 18 May 2021b. a
GEOframe Summer School: http://geoframe.blogspot.com/2021/06/geoframesummerschool2021gss2021.html, last access: 21 December 2021. a
Germann, P. and Beven, K.: Water flow in soil macropores I. An experimental approach, J. Soil Sci., 32, 1–13, https://doi.org/10.1111/j.13652389.1981.tb01681.x, 1981. a
Greenspan, D. and Casulli, V.: Numerical Analysis for Applied Mathematics, Science, and Engineering, Addison Wesley, https://doi.org/10.1201/9780429493393, 1988. a
Gugole, F., Dumbser, M., and Stelling, G.: An efficient threedimensional semiimplicit finite volume scheme for the solution of coupled freesurface and variably saturated subsurface flow, in: EGU General Assembly Conference Abstracts, 600, 2018. a
Hall, C. A., Saia, S. M., Popp, A. L., Dogulu, N., Schymanski, S. J., Drost, N., van Emmerik, T., and Hut, R.: A Hydrologist's Guide to Open Science, Hydrol. Earth Syst. Sci. Discuss. [preprint], https://doi.org/10.5194/hess2021392, in review, 2021. a
Hay, L. E., Leavesley, G. H., Clark, M. P., Markstrom, S. L., Viger, R. J., and Umemoto, M.: Step wise, multiple objective calibration of a hydrologic model for a snowmelt dominated basin 1, J. Am. Water Resour. As., 42, 877–890, https://doi.org/10.1111/j.17521688.2006.tb04501.x, 2006. a, b
Horton, R. E.: The role of infiltration in the hydrologic cycle, EOS T. Am. Geophys. Un., 14, 446–460, https://doi.org/10.1029/TR014i001p00446, 1933. a
Hrachowitz, M., Benettin, P., Van Breukelen, B. M., Fovet, O., Howden, N. J., Ruiz, L., Van Der Velde, Y., and Wade, A. J.: Transit times–The link between hydrology and water quality at the catchment scale, WIREs Water, 3, 629–657, https://doi.org/10.1002/wat2.1155, 2016. a
IDEAS: Methodologies and “How To”, available at: https://ideasproductivity.org/ideasclassic/howto/, last access: 20 October 2021. a
Jones, J. E. and Woodward, C. S.: Newton–Krylovmultigrid solvers for largescale, highly heterogeneous, variably saturated flow problems, Adv. Water Resour., 24, 763–774, https://doi.org/10.1016/S03091708(00)000750, 2001. a
Kelley, C. T.: Solving nonlinear equations with Newton's method, SIAM, https://doi.org/10.1137/1.9780898718898, 2003. a
Kennedy, J. and Eberhart, R.: Particle swarm optimization, in: Proceedings of ICNN'95international conference on neural networks, Perth, WA, Australia, 27 November–1 December 1995, IEEE, vol. 4, 1942–1948, 1995. a
Kollet, S. J. and Maxwell, R. M.: Integrated surface–groundwater flow modeling: A freesurface overland flow boundary condition in a parallel groundwater flow model, Adv. Water Resour., 29, 945–958, https://doi.org/10.1016/j.advwatres.2005.08.006, 2006. a
Kosugi, K.: General model for unsaturated hydraulic conductivity for soils with lognormal poresize distribution, Soil Sci. Soc. Am. J., 63, 270–277, https://doi.org/10.2136/sssaj1999.03615995006300020003x, 1999. a
Kurylyk, B. L. and Watanabe, K.: The mathematical representation of freezing and thawing processes in variablysaturated, nondeformable soils, Adv. Water Resour., 60, 160–177, https://doi.org/10.1016/j.advwatres.2013.07.016, 2013. a
Liu, S., Lu, L., Mao, D., and Jia, L.: Evaluating parameterizations of aerodynamic resistance to heat transfer using field measurements, Hydrol. Earth Syst. Sci., 11, 769–783, https://doi.org/10.5194/hess117692007, 2007. a
Lloyd, W., David, O., Ascough II, J. C., Rojas, K. W., Carlson, J. R., Leavesley, G. H., Krause, P., Green, T. R., and Ahuja, L. R.: Environmental modeling framework invasiveness: Analysis and implications, Environ. Modell. Softw., 26, 1240–1250, https://doi.org/10.1016/j.envsoft.2011.03.011, 2011. a
Lu, N.: Generalized soil water retention equation for adsorption and capillarity, J. Geotech. Geoenviron., 142, 04016051, https://doi.org/10.1061/(ASCE)GT.19435606.0001524, 2016. a
Martin, R. C.: Clean code: a handbook of agile software craftsmanship, Pearson Education, ISBN13: 9780132350884, 2009. a, b
Mcdonough, J. M.: Introductory Lectures on Turbulence: Physics, Mathematics and Modeling, Mechanical Engineering Textbook Gallery, 2, available at: https://uknowledge.uky.edu/me_textbooks/2 (last access: 21 December 2021), 2007. a
Menard, C. B., Essery, R., Krinner, G., Arduini, G., Bartlett, P., Boone, A., BrutelVuilmet, C., Burke, E., Cuntz, M., Dai, Y., Decharme, B., Dutra, E., Fang, X., Fierz, C., Gusev, Y., Hagemann, S., Haverd, V., Kim, H., Lafaysse, M., Marke, T., Nasonova, O., Nitta, T., Niwano, M., Pomeroy, J., Schädler, G., Semenov, V. A., Smirnova, T., Strasser, U., Swenson, S., Turkov, D., Wever, N., and Yuan, H.: Scientific and human errors in a snow model intercomparison, B. Am. Meteorol. Soc., 102, E61–E79, https://doi.org/10.1175/BAMSD190329.1, 2020. a
Mualem, Y.: A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res., 12, 513–522, https://doi.org/10.1029/WR012i003p00513, 1976. a
Muskat, M. and Meres, M. W.: The flow of heterogeneous fluids through porous media, Physics, 7, 346–363, https://doi.org/10.1063/1.1745403, 1936. a
National Research Council: Basic Research Opportunities in Earth Science, The National Academies Press, Washington, DC, https://doi.org/10.17226/9981, 2001. a
Newman, S.: Building Microservices: Designing FineGrained Systems, ISBN10: 1491950358, 2015. a
Nicolsky, D. and Romanovsky, V. E.: Modeling longterm permafrost degradation, J. Geophys. Res.Earth, 123, 1756–1771, https://doi.org/10.1029/2018JF004655, 2018. a
Nobel, P. S.: Physicochemical & environmental plant physiology, Academic Press, ISBN 0125200250, 1999. a
Ochsner, T. E., Horton, R., and Ren, T.: A new perspective on soil thermal properties, Soil Sci. Soc. Am. J., 65, 1641–1647, https://doi.org/10.2136/sssaj2001.1641, 2001. a
Paniconi, C. and Putti, M.: A comparison of Picard and Newton iteration in the numerical solution of multidimensional variably saturated flow problems, Water Resour. Res., 30, 3357–3374, https://doi.org/10.1029/94WR02046, 1994. a, b
Paniconi, C. and Putti, M.: Physically based modeling in catchment hydrology at 50: Survey and outlook, Water Resour. Res., 51, 7090–7129, https://doi.org/10.1002/2015WR017780, 2015. a
Paniconi, C. and Wood, E. F.: A detailed model for simulation of catchment scale subsurface hydrologic processes, Water Resour. Res., 29, 1601–1620, https://doi.org/10.1029/92WR02333, 1993. a
PC Progress Discussion Forums: https://www.pcprogress.com/forum/viewtopic.php?f=3&t=3632, last access: 20 April 2021. a
Peckham, S. D., Hutton, E. W., and Norris, B.: A componentbased approach to integrated modeling in the geosciences: The design of CSDMS, Comput. Geosci., 53, 3–12, https://doi.org/10.1016/j.cageo.2012.04.002, 2013. a
Penna, D., Hopp, L., Scandellari, F., Allen, S. T., Benettin, P., Beyer, M., Geris, J., Klaus, J., Marshall, J. D., Schwendenmann, L., Volkmann, T. H. M., von Freyberg, J., Amin, A., Ceperley, N., Engel, M., Frentress, J., Giambastiani, Y., McDonnell, J. J., Zuecco, G., Llorens, P., Siegwolf, R. T. W., Dawson, T. E., and Kirchner, J. W.: Ideas and perspectives: Tracing terrestrial ecosystem water fluxes using hydrogen and oxygen stable isotopes – challenges and opportunities from an interdisciplinary perspective, Biogeosciences, 15, 6399–6415, https://doi.org/10.5194/bg1563992018, 2018. a
Priestley, C. H. B. and Taylor, R. J.: On the assessment of surface heat flux and evaporation using largescale parameters, Mon. Weather Rev., 100, 81–92, https://doi.org/10.1175/15200493(1972)100<0081:OTAOSH>2.3.CO;2, 1972. a
Quarteroni, A., Sacco, R., and Saleri, F.: Numerical mathematics, vol. 37, Springer Science and Business Media, ISBN10: 3642071015, 2010. a
Raupach, M. and Thom, A. S.: Turbulence in and above plant canopies, Annu. Rev. Fluid Mech., 13, 97–129, https://doi.org/10.1146/annurev.fl.13.010181.000525, 1981. a
Regenass, D., Schlemmer, L., Jahr, E., and Schär, C.: It rains and then? Numerical challenges with the 1D Richards equation in kilometerresolution land surface modelling, Hydrol. Earth Syst. Sci. Discuss. [preprint], https://doi.org/10.5194/hess2021426, in review, 2021. a
Richards, L. A.: Capillary conduction of liquids through porous mediums, Physics, 1, 318–333, https://doi.org/10.1063/1.1745010, 1931. a
Richardson, L. F.: Weather prediction by numerical process, Cambridge University Press, https://doi.org/10.1002/qj.49704820311, 1922. a
Rigon, R., Bertoldi, G., and Over, T. M.: GEOtop: A distributed hydrological model with coupled water and energy budgets, J. Hydrometeorol., 7, 371–388, https://doi.org/10.1175/JHM497.1, 2006a. a, b, c
Rigon, R., Ghesla, E., Tiso, C., and Cozzini, A.: The HORTON machine: a system for DEM analysis The reference manual, Tech. rep., University of Trento, ISBN 9788884431479, available at: http://hdl.handle.net/11572/48576 (last access: 21 December 2021), 2006b. a
Rigon, R., Bancheri, M., Formetta, G., and de Lavenne, A.: The geomorphological unit hydrograph from a historicalcritical perspective, Earth Surf. Proc. Land., 41, 27–37, https://doi.org/10.1002/esp.3855, 2016a. a
Rigon, R., Bancheri, M., and Green, T. R.: Ageranked hydrological budgets and a travel time description of catchment hydrology, Hydrol. Earth Syst. Sci., 20, 4929–4947, https://doi.org/10.5194/hess2049292016, 2016b. a
Rizzoli, A., Svensson, M., Rowe, E., Donatelli, M., Muetzelfeldt, R., van der Wal, T., van Evert, F., and Villa, F.: Modelling framework (SeamFrame) requirements, Tech. rep., SEAMLESS, 2006. a
Roe, P. L.: Approximate Riemann solvers, parameter vectors, and difference schemes, J. Comput. Phys., 43, 357–372, https://doi.org/10.1016/00219991(81)901285, 1981. a
Romano, N., Brunone, B., and Santini, A.: Numerical analysis of onedimensional unsaturated flow in layered soils, Adv. Water Resour., 21, 315–324, https://doi.org/10.1016/S03091708(96)000590, 1998. a, b
Ronan, A. D., Prudic, D. E., Thodal, C. E., and Constantz, J.: Field study and simulation of diurnal temperature effects on infiltration and variably saturated flow beneath an ephemeral stream, Water Resour. Res., 34, 2137–2153, https://doi.org/10.1029/98WR01572, 1998. a, b
Saito, H., Šimůnek, J., and Mohanty, B. P.: Numerical analysis of coupled water, vapor, and heat transport in the vadose zone, Vadose Zone J., 5, 784–800, https://doi.org/10.2136/vzj2006.0007, 2006. a
Serafin, F.: Enabling modeling framework with surrogate modeling capabilities and complex networks, PhD thesis, University of Trento, 2019. a, b, c, d, e, f
Shewchuk, J. R.: An Introduction to the Conjugate Gradient Method Without the Agonizing Pain, available at: https://www.bibsonomy.org/bibtex/2175fea4e807e879ebcb76b79a9992437/timo (last access: 21 December 2021), 1994. a
Šimůnek, J., Van Genuchten, M. T., and Sejna, M.: The HYDRUS1D software package for simulating the onedimensional movement of water, heat, and multiple solutes in variablysaturated media, University of CaliforniaRiverside Research Reports, 3, 1–240, 2005. a, b
Šimůnek, J., Van Genuchten, M. Th., and Šejna, M.: The HYDRUS1D software package for simulating the onedimensional movement of water, heat, and multiple solutes in variablysaturated media. Version 3.0, HYDRUS Softw., Ser. 1, Dep. of Environ. Sci., Univ. of California, Riverside, CA, 2012. a
Sophocleous, M.: Analysis of water and heat flow in unsaturatedsaturated porous media, Water Resour. Res., 15, 1195–1206, https://doi.org/10.1029/WR015I005P01195, 1979. a
Srivastava, R. and Yeh, T.C. J.: Analytical solutions for onedimensional, transient infiltration toward the water table in homogeneous and layered soils, Water Resour. Res., 27, 753–762, https://doi.org/10.1029/90WR02772, 1991. a, b, c, d
Tomin, P. and Lunati, I.: Investigating Darcyscale assumptions by means of a multiphysics algorithm, Adv. Water Resour., 95, 80–91, https://doi.org/10.1016/j.advwatres.2015.12.013, 2016. a
Tubini, N. and Rigon, R.: GEOframePy, available at: https://pypi.org/project/geoframepy/, last access: 22 October 2021a. a
Tubini, N. and Rigon, R.: OMS project for WHETGEO1D, v1.0beta, Zenodo [data set], https://doi.org/10.5281/zenodo.4749319, 2021b (data available at: https://github.com/GEOframeOMSProjects/OMS_Project_WHETGEO1D, last access: 21 December 2021). a, b
Tubini, N. and Rigon, R.: WHETGEO1D v1.0beta, Zenodo [code], https://doi.org/10.5281/zenodo.5112727, 2021c (data available at: https://github.com/geoframecomponents/WHETGEO1D, last access: 21 December 2021). a
Tubini, N., Formetta, G., Rigon, R., Bancheri, M., Serafin, F., Bottazzi, M., Tasin, S., Dalla Torre, D., and D'Amato, C.: Day 12 WHETGEO1D, OSF [code], https://doi.org/10.17605/OSF.IO/Z53C6, 2021a. a
Tubini, N., Gruber, S., and Rigon, R.: A method for solving heat transfer with phase change in ice or soil that allows for large time steps while guaranteeing energy conservation, The Cryosphere, 15, 2541–2568, https://doi.org/10.5194/tc1525412021, 2021b. a, b, c, d, e, f, g
Unidata: Chunking Data with NetCDF4, available at: https://www.unidata.ucar.edu/software/netcdf/workshops/2008/nc4chunking/index.html, last access: 1 April 2021a. a
Unidata: Formats and Performance, available at: https://www.unidata.ucar.edu/software/netcdf/workshops/2008/performance/index.html, last access: 1 April 2021b. a
Unidata: NetCDF version 4.3.22, UCAR/Unidata Program Center, Boulder, CO, https://doi.org/10.5065/D6H70CW6, 2021. a
Van Genuchten, M. T.: A closedform equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44, 892–898, https://doi.org/10.2136/sssaj1980.03615995004400050002x, 1980. a
Vanderborght, J., Kasteel, R., Herbst, M., Javaux, M., Thiéry, D., Vanclooster, M., Mouvet, C., and Vereecken, H.: A set of analytical benchmarks to test numerical models of flow and transport in soils, Vadose Zone J., 4, 206–221, https://doi.org/10.2113/4.1.206, 2005. a, b, c
Walvoord, M. A. and Kurylyk, B. L.: Hydrologic impacts of thawing permafrost–A review, Vadose Zone J., 15, 1–20, https://doi.org/10.2136/vzj2016.01.0010, 2016. a
Walvoord, M. A., Voss, C. I., and Wellman, T. P.: Influence of permafrost distribution on groundwater flow in the context of climatedriven permafrost thaw: Example from Yukon Flats Basin, Alaska, United States, Water Resour. Res., 48, W07524, https://doi.org/10.1029/2011WR011595, 2012. a
Whitaker, S.: Flow in porous media I: A theoretical derivation of Darcy's law, Transport Porous Med., 1, 3–25, https://doi.org/10.1007/BF01036523, 1986. a
Wierenga, P., Hagan, R. M., and Nielsen, D.: Soil temperature profiles during infiltration and redistribution of cool and warm irrigation water, Water Resour. Res., 6, 230–238, https://doi.org/10.1029/WR006i001p00230, 1970. a
Yang, Y., Wu, J., Zhao, S., Han, Q., Pan, X., He, F., and Chen, C.: Assessment of the responses of soil pore properties to combined soil structure amendments using Xray computed tomography, Sci. Rep.UK, 8, 695–705, https://doi.org/10.1038/s41598017189971, 2018. a
Zha, Y., Yang, J., Zeng, J., Tso, C.H. M., Zeng, W., and Shi, L.: Review of numerical solution of Richardson–Richards equation for variably saturated flow in soils, WIREs Water, 6, e1364, https://doi.org/10.1002/wat2.1364, 2019. a, b, c, d, e, f
Zhang, S., Meurey, C., and Calvet, J.C.: Identification of soilcooling rains in southern France from soil temperature and soil moisture observations, Atmos. Chem. Phys., 19, 5005–5020, https://doi.org/10.5194/acp1950052019, 2019. a
 Abstract
 Introduction
 Mathematical issues and numerics
 Design and deployment of WHETGEO1D
 Information for users and developers
 Conclusions
 Appendix A: GEOframe
 Appendix B: OMS3
 Appendix C: R^{2} test cases
 Appendix D: Energy budget
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Mathematical issues and numerics
 Design and deployment of WHETGEO1D
 Information for users and developers
 Conclusions
 Appendix A: GEOframe
 Appendix B: OMS3
 Appendix C: R^{2} test cases
 Appendix D: Energy budget
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement