HydroGeoSphere: A Fully Integrated, Physically Based Hydrological Model

Introduction The importance of a quantitative understanding of the hydrological cycle increases with the ever-growing demand for water for anthropogenic needs. Numerical models are inevitable tools in this undertaking. A wide range of numerical models of different complexity have been developed for this purpose, ranging from simple, lumped parameter models to more complex, physically based models. The foundation of physically based models is the blueprint paper by Freeze and Harlan (1969), and numerous physically based models have been developed following this blueprint. HydroGeoSphere (HGS), the code discussed in this review, is one of them. The origin of HGS is the code FRAC3DVS, developed by R. Therrien at the University of Waterloo as part of his doctoral work under the supervision of E.A. Sudicky (Therrien 1992). FRAC3DVS was designed to simulate variably saturated groundwater flow and advective-dispersive solute transport in porous or discretely fractured porous media. In 2002, a two-dimensional (2D) surface water flow and transport component were implemented in FRAC3DVS and the code was renamed HydroGeoSphere. Until recently, the code was free for academic research, while commercial users paid a license fee between 3000 and 6000 US dollars depending on the number of CPU cores the code will use in a parallel computational platform. The code can be downloaded by contacting the developers through the website: http://hydrogeosphere.org/. HGS has been designed to solve simple problems (e.g., regular geometry, steady-state saturated flow etc.) as well as very complex problems (e.g., integrated flow,


Introduction
The importance of a quantitative understanding of the hydrological cycle increases with the ever-growing demand for water for anthropogenic needs. Numerical models are inevitable tools in this undertaking. A wide range of numerical models of different complexity have been developed for this purpose, ranging from simple, lumped parameter models to more complex, physically based models. The foundation of physically based models is the blueprint paper by Freeze and Harlan (1969), and numerous physically based models have been developed following this blueprint. HydroGeoSphere (HGS), the code discussed in this review, is one of them. The origin of HGS is the code FRAC3DVS, developed by R. Therrien at the University of Waterloo as part of his doctoral work under the supervision of E.A. Sudicky (Therrien 1992). FRAC3DVS was designed to simulate variably saturated groundwater flow and advective-dispersive solute transport in porous or discretely fractured porous media. In 2002, a two-dimensional (2D) surface water flow and transport component were implemented in FRAC3DVS and the code was renamed HydroGeoSphere. Until recently, the code was free for academic research, while commercial users paid a license fee between 3000 and 6000 US dollars depending on the number of CPU cores the code will use in a parallel computational platform. The code can be downloaded by contacting the developers through the website: http://hydrogeosphere.org/.
HGS has been designed to solve simple problems (e.g., regular geometry, steady-state saturated flow etc.) as well as very complex problems (e.g., integrated flow, solute, and heat transport in regional-scale catchments). It has been applied to a wide range of hydrological problems including large-scale watershed modeling (Unger et al. 2008), the disconnection between surface water and groundwater (Brunner et al. 2009), fractured rock hydrology (Park et al. 2004), the simulation of the fate of contaminants in groundwater ), atmospheric/surface/subsurface heat exchange processes along river reaches (Brookfield et al. 2009) the behavior of solutes in surface water (McCallum et al. 2010), the dynamics of river bank storage processes (Doble et al. 2011), low-temperature geothermal energy (Raymond and Therrien 2008), the simulation of brines and seawater intrusion (Graf and Therrien 2008), and dual permeability systems (Schwartz et al. 2010). The code has also been applied to glaciation processes (Lemieux et al. 2008) and oil sands reclamation (Carrera-Hernandez et al. 2011). In all of these fields, HGS has been used as a tool to model field sites or catchments, but also as a tool to understand the basic physics of the underlying hydrological phenomena. It is our own observation that the users of HGS are mainly found in academia, but the number of users in industry is certainly growing.

Modeling Capabilities
The most important feature of HGS is its ability to simulate water flow in a fully integrated mode, thus allowing precipitation to partition into all key components of the hydrologic cycle. In a physically based, natural way, rainfall input will partition into overland and streamflow, evaporation, transpiration, groundwater recharge, or subsurface discharge into surface water bodies such as rivers or lakes. However, the program's capabilities are not limited to the "classical" processes of the hydrological cycle. Variable-density flow and transport, straight, or branching chain first-order decay reactions, reactive chemical species transport, heat transport, unsaturated flow and flow through fractures can also be simulated by HGS in a physically based way. In the following, these features are discussed in greater detail.
The numerical formulation of HGS is based on the assumption that a subsurface flow equation for a porous medium (either saturated or unsaturated) is always solved. A three-dimensional (3D) modified formulation of the Richards' equation is implemented. To relate pressure head to saturation and relative hydraulic conductivity, either lookup tables, the van Genuchten (1980) or Brooks Corey (1964) relations can be used. The simulation of variably saturated flow in a dual continuum is also possible, based on the formulation of Gerke and van Genuchten (1993). The second continuum could represent fractures or macropores. A range of subsurface boundary conditions are available, including boundaries of prescribed hydraulic head (Dirichlet), prescribed flux (Neuman), drains, and seepage faces. Seepage faces are implemented following Forsyth (1988). HGS simulates the abstraction of groundwater through wells using the formulation of Therrien and Sudicky (2000). Flow in wells (penetrating a variably saturated aquifer) is described with a 1d free surface flow equation formulated along the axis of the well. Flow in tile drains is based on the continuity equation of open channel flow. The fluid flux for the tile drain can be approximated with the friction formulas or the diffusive wave approximation to the shallow water flow equations.
Surface water flow is simulated using a 2D, depthaveraged flow equation (the diffusion-wave approximation of the Saint Venant equations). To account for smallscale topographic variations the concepts of rill storage and storage exclusion are implemented. Rill storage is the amount of storage that must be filled before lateral surface flow can occur. Storage exclusion allows accounting for a reduced surface domain porosity (e.g., through the presence of buildings). Channel flow can also be simulated and the model formulation is analogous to one of tile drains. Two approaches for coupling the surface to the subsurface are available. The first one is the common node approach and it is based on the assumption of continuity of hydraulic head between the two domains (Therrien and Sudicky 1996). The second approach is the so-called dual node approach and is based on a first-order exchange coefficient. This is implemented in HGS using the concept of a coupling length described by Park et al. (2009). For an ongoing discussion of this approach refer to Ebel et al. (2009).
A range of boundary conditions can be assigned to the surface domain. Areal rainfall is implemented as an input water flux multiplied by the contributing surface area. Interception, evaporation, and transpiration are modeled following the conceptualization of Kristensen and Jensen (1975). In this approach, transpiration and evaporation are calculated as a function of moisture content and a range of parameters describing the type of vegetation (e.g., leaf area index). Boundaries constraining surface water flow include a critical depth boundary and a zero-depth gradient boundary. The zero-depth gradient condition forces the slope of the water level to equal the bed slope (defined by the user). The discharge at this boundary is calculated using Manning's equation.
HGS incorporates the capability of simulating nonreactive as well as reactive chemical species transport in the surface and the subsurface domain as well as in discrete fractures, a second porous continuum, wells, and tile drains. Chain decay reactions and isotopic fractionation can also be simulated.
HGS is fully capable of simulating thermal energy transport. The implementation of heat transport in the fully saturated domain is described by Graf (2005) and Graf and Therrien (2007). The user can choose if changes in heat and solute concentrations affect fluid properties such as density or viscosity. Density-driven flow is now coupled with unsaturated flow, following the conceptualization of Graf and Boufadel (2011). Brookfield et al. (2009) further extended the heat-transport capabilities in HGS to include energy transport in the subsurface under variably saturated conditions as well as in the surface water domain. The ambient air temperature and incoming short-and long wave solar radiation are used to drive the thermal system. The approach is based on the formulation of Verseghy (1991). If energy exchange with the atmosphere is simulated, the user can either predefine potential evaporation or simulate it based on the energy exchange fluxes with the atmosphere. The full details on the implementation of energy and solute transport are discussed by Graf and Therrien (2007) and Brookfield et al. (2009).
HGS uses numerical methods that are designed for computational efficiency and robustness. For both the surface and the subsurface domains, a control volume finite element method (Forsyth 1991) is employed to discretize the surface and subsurface flow equations. Transport equations are solved with the standard Galerkin element method or the control volume finite element method. Alternatively, HGS can mimic a 7-point finite difference scheme for either flow or transport by manipulating the influence coefficients that arise from a mesh composed of hexahedral elements. Different methods to solve the discretized flow and transport equations are available. The system of algebraic equations (matrix equation) arising from discretization is solved by a preconditioned iterative solver. Preconditioning is done by performing an incomplete LU factorization of the coefficient matrix. The ORTHOMIN accelerator by Behie and Forsyth (1984) is used during the iteration. Acceleration techniques available include conjugate gradient (CG), generalized minimal residual (GMRES), and biconjugate gradient stabilized (BiCGSTAB) methods.
Nonlinear equations are linearized using the Newton-Raphson technique as described by Huyakorn and Pinder (1983). A noteworthy feature of HGS is the primary variable substitution in solving the Richards' equation. The Richards' equation has improved convergence properties if the primary variable is saturation rather than pressure head in dry regions of the domain (Forsyth et al. 1995). The approach used in HGS is to switch the primary variable in wet or fully saturated regimes to pressure head, while retaining saturation as the primary variable in drier portions of the unsaturated zone. HGS allows for adaptive time stepping following the approach by Forsyth and Sammon (1986). HGS is written in FORTRAN. It runs on all versions of Windows and Linux systems and is available in both 32 and 64 bit. A program installer can be downloaded from the University of Waterloo, and a password to access the download area is required. The installer contains four executables, as well as documentation, test, and verification cases. Once installed, the folder containing HGS and all supporting files is only about 20 MB in size. Along with the executables, a manual is provided with comprehensive information on the conceptualization of the model, the input instructions, as well as the detailed description of 30+ illustration examples. All input files required to run the examples are provided as well. The illustration examples cover a range of typical benchmark problems.
The code is being continuously developed. According to a presentation by the developers at the HGS-user conference 2011 in Hannover, a range of new features are currently being implemented. Notable developments include a full coupling to atmospheric models, OpenMP parallelization, the development of a multithread parallel version for Windows-based PCs, subgridding and subtiming, the simulation of snowmelt/thaw capability, 1-d channel flow with flow-control features such as dams and weirs and additional options to directly include ARC-GISbased data in the simulation. The extraction of baseflow following the method of Partington et al. (2011) is currently being implemented. HGS continues to advance the range of physical processes which can be simulated, and is taking full advantage of the very latest numerical and computational advances to enhance computational efficiency and numerical robustness.

Pre-and Postprocessing
HGS does currently not have a graphical user interface (GUI). All model parameters, grid structures, material properties, or numerical parameters are written in text files. A preprocessor (called GROK) prepares the input files for HGS. The main input file is an instructiondriven text file that only requires a text editor. It is very short for simple problems but can also be designed for complex problems. Once the preprocessing is completed, HGS is run. Once the simulation has terminated, a postprocessing program called HSPLOT converts the output data to a format that can be read by third-party visualization packages such as TECPLOT (Tecplot Inc. 2011) or GMS (AquaVeo 2011).
Key tasks carried out by GROK are memory allocation, grid generation, and the assignment of properties to all parameters. GROK does not allocate memory automatically, but the user can easily increase the upper limits assigned to memory by manipulating a text file. A range of options to construct and manipulate the grid are available. GROK can construct grids composed of hexahedral blocks, triangular prisms, and tetrahedra. 2D triangular grids can be imported from GRIDBUILDER or GMS. GRIDBUILDER is an efficient tool to generate complex meshes and has been developed in parallel to HGS at the University of Waterloo by Rob McLaren (GRIDBUILDER is also used by the model FEFLOW). It offers a variety of options to read georeferenced raster data and vector data. Such data can, for example, be used to define the geometric boundaries of the model domain.
Once the grid has been generated by GRIDBUILDER, HGS offers the option to re-project the coordinate system, for example, from a spherical to a planar 2d system. To manipulate properties of the grid or to assign boundary conditions, subsets of the grid have to be selected through the input files. A range of options to select nodes, elements, or faces are available to the user. Nodes or elements can, for example, be selected according to geometric criteria (e.g., nodes/elements within a given block or a defined sphere are selected). Alternatively, nodes, elements, or faces can be selected in GRIDBUILDER. This is especially convenient because GRIDBUILDER can handle georeferenced ARC-GIS data. The chosen elements can be exported by GRIDBUILDER, and specific properties can then be assigned to the chosen grid components. In case the user does not define a required parameter value explicitly, GROK assigns default properties. Once all input files are prepared, GROK is run to prepare the input files for HGS. This process is normally fast (e.g., a few seconds) for grids smaller than 100,000 nodes. Memory is allocated dynamically in HGS based on the problem specifications written in the GROK output files, which minimizes the amount of memory required by HGS. After preparing the input files for HGS, GROK writes a file that summarizes all the relevant information. In case of an error in the input file, the summary file stops at the instruction that caused the problem. A sample input file is presented in Figure 1.
After completing the preprocessing step, HGS can be executed. A very useful option in HGS is that the user can change a range of numerical parameters (e.g., the parameters controlling the time stepping, the data output frequency, number of allowed iterations etc.) during the simulation. During the simulation HGS writes a series of ASCII output files, such as detailed numerical information for all timesteps, the global water balance or simulation data such as hydraulic heads and solute concentrations. The user can define observation boreholes where HGS writes detailed information on all simulation data. Data such as hydraulic heads, flow velocities, or solute concentrations are stored in binary files. In order to visualize and further interpret these results, a program called HSPLOT must be executed. HSPLOT converts the binaries to ASCII files that can be imported by third-party software packages such as TECPLOT or GMS. Besides allowing for the visualization of data, such packages allow for various postprocessing steps such as the integration of fluid fluxes over time.
Having described both the model capabilities and pre-and postprocessing steps, we are now in a position to reflect on the learning curve for HGS users. In our own experience, we believe that the learning curve for fairly simple problems is not too steep. However, the learning curve is much steeper in the case of complex regional scale models and situations which extend beyond only groundwater flow modelling e.g., fully coupled surface and groundwater flow, or thermohaline transport processes. This general observation is expected to be true of any modelling environment whether it contains or does not contain a GUI and therefore, in our opinion, does not pertain specifically to HGS. In addition, new users have also to consider the time and knowledge required to generate grids through GRIDBUILDER or GMS. Additionally, for advanced postprocessing and visualization, the user has to learn how to use packages such as TECPLOT. While the TECPLOT package offers limitless options for visualization and interpretations of the results, a considerable amount of time is required to become familiar with the program.

What We Like
The capabilities of HGS are impressive. Through HGS, the user has a powerful and versatile tool that can be applied to an extremely wide field of hydrological problems. The numerical methods employed are fast, robust, and stable. The numerous, peer-reviewed papers (40+) demonstrate that the model stands up to its claims. Its integrated nature allows water to flow naturally, therefore the definition of "unnatural" boundary conditions is not required. For example, a river will form naturally within the model and interact with the groundwater in a physically based way. There is no need to predefine the river boundaries or its hydraulic head. Nevertheless, the user can impose such boundary conditions if he/she chooses to do so. Because of its original development, the model is very well suited for subsurface simulations in fractured rock (porous or nonporous) and it incorporates the most commonly used conceptual models for fractured rock (i.e., equivalent porous medium, dual porosity, discrete fracture network). This is one of the many features of HGS that is not included in other fully integrated models. That such a powerful tool can be provided in a single platform is remarkable.
The way the input files are written might require some time for users to get used to. However, they are intuitive and very flexible, because no predefined structure of files is required. For example, once the grid is generated, the sections as shown Figure 1 can be arranged in any order. The default simulation setup is for steady-state saturated groundwater flow and the user activates modeling options by adding simple instructions in the GROK file such as "transient flow" (Figure 1) or "solute transport." As a result, the user does not need to be concerned with input data for processes that are not simulated, as shown by the pumping test example in Figure 1 where there is no input data for options such as solute transport or surface water flow. The preprocessor writes all relevant information to clearly structured ASCII files, and the user can rapidly check all aspects of the problem. Another positive feature is the highly flexible way to assign properties and boundary conditions. Because of the intuitive nature of the input instructions, the lack of a GUI is, in our opinion, not a negative point at all because it allows the user to take full control of the input. Therefore, a very high level of transparency is achieved throughout the modeling process. It also forces new users to understand how the code is structured and works. The text file structure also allows implementing of model parameter changes rapidly, and also makes the coupling to parameter estimation codes such as PEST (Doherty 2010) straightforward.
There is a growing community of HGS users worldwide. Since 2010, annual HGS meetings have been organized by the user community. The meetings offer a good opportunity to exchange experiences with the code and discuss user requirements with the developers who are very accessible. The recently opened user forum on the HGS homepage (http://hydrogeosphere.org/) will hopefully benefit from the rapidly increasing number of users. The supporting staff at Waterloo are incredibly efficient and helpful when it comes to answering e-mails and giving support on complex questions. They can be directly contacted through the HGS homepage. They also regularly release updated versions of the code that can be downloaded by users from an ftp site.
Another positive point is that the software can be installed without any problems whatsoever. Also, the stability of Windows is not affected through use of HGS. We have not observed a single-forced termination of the executables, even during very long execution times (e.g., several days). HGS requires significant computational power, which can slow down single core machines considerably. However, on dual or quad core machines the user can work in parallel to HGS.

What We Do Not Like
One of the weaker points of HGS is that the manual is not updated synchronously with the code developments. A number of the code's capabilities are not documented (e.g., permafrost), and parts of the code that are no longer supported are still documented in the manual (e.g., options on how evaporation is calculated). Very little information on the numerical capabilities is provided, and it is therefore unlikely that users can take full advantage of the numerous numerical techniques included in the code. (A good overview to some of the numerical concepts is given by Rausch et al. 2005.) Several references that are listed in the body of the manual are not found in the reference list. The documentation of GRIDBUILDER suffers from the same shortcomings. For example, GRIDBUILDER provides extremely useful functions to manipulate the digital elevation model along a river, but these functions are not documented.
Another challenging point is that multiple versions of HGS exist and it is hard for the user to stay on top of the changes, new bugs, bug fixes, and updates. Also, some functions no longer work in updated versions. For example, the axisymmetric grid generation does not function as described in the manual, and the option of primary variable substitution does not work in certain configurations. A comprehensive and continuously updated list with known bugs and limitations would really help users avoid traps, save time, and avoid potential frustration.
Although we do not consider the "outsourcing" of visualization and postprocessing as a limitation, the user has to keep in mind that these third-party products are quite expensive, and that some time to get familiar with TECPLOT or GMS is required. Beginners would also benefit from a short tutorial on how to import and manipulate HGS-output data in such programs.
Also, we miss the option to write zoned water budgets. While an option is provided to extract water fluxes through slices, this option does not work for the extraction of surface-subsurface exchange fluxes or for subareas of transpiration/evaporation rates. While these fluxes can be extracted via TECPLOT, a direct option would be required to use the output data in parameter estimation processes. Another small but in some cases immensely helpful option would be to allow the postprocessor HSPLOT to write HGS-output files directly to TECPLOT or GMS format, without taking the detour of converting the binaries to often extremely large ASCII files.

Discussion and Conclusions
The blueprint of Freeze and Harlan (1969) was based on the vision to integrate the critical components of the hydrological cycle in one conceptual framework. HGS is one of the most powerful codes realizing this vision, with modeling capabilities that go far beyond the original blueprint: it can be used in virtually any field of hydrological sciences. Planned developments on integrating additional components of the physical environment such as the atmosphere will ensure that the code continues to overcome the artificial and limiting conceptual separations between hydrological sub-disciplines.
In addition to the flexibility and features of the code outlined earlier, it is important to point out that the code developers are extremely accessible and helpful and very open-minded to discussions on weaknesses. It is our observation that they are actively improving HGS all the time on virtually a good will basis. We encourage new users to give it a try, and to invest the time required to master the handling of the code. They will be rewarded with a powerful tool that allows them to model and to explore the underlying physics of the hydrological cycle. The balancing of stronger and weaker points demonstrates that the strengths by far outweigh the weaknesses. Although it is positive that the code has been virtually free in the past for academic use, it also means that the resources to address some of the above mentioned limitations (most importantly developing a perpetually updated manual) are limited. We believe that increasing the fees and investing the additional resources to overcome these limitations would be highly beneficial for the HGS community and would potentially make the code more attractive for industry users.
To finalize this review, we think it is useful to briefly comment on the ongoing debate in the hydrological community on complexity and simplicity in hydrogeological modeling. Obviously, HGS is a complex code, and therefore a prime target to criticism (e.g., Beven 2002) aimed at the complexity of conceptual models. A justified and often raised criticism of complex models is that the amount of parameters required often exceeds the resources and capabilities to obtain them, and the long run times of physically based models make parameter estimation difficult or even impossible. On the other hand, simplified models may not offer sufficient insights into the underlying physics and accordingly, the errors associated with the simplifications are often unknown. In fact, to quantify the systematic errors associated with such simplifications, more complex physically based models are an extremely useful tool. We do not believe that anyone can win the debate on simplicity and complexity, as countless solid points can be made on both sides. We recall Einstein's famous quote "Everything should be made as simple as possible but never simpler." Knowing how simple is "too simple" is a key challenge for hydrology in the 21st century. This assessment will require rigorous and quantitative evaluation to determine appropriate model formulation (and appropriate complexity or simplicity) relative to model purpose. The model, whether complex or simple, must ultimately, reliably, and robustly answer a particular question that is the intended purpose of the modeling exercise. In this sense, having access to and being able to master a complex model puts the user in the comfortable position to adjust the level of complexity as required to address the problem. We prefer to have a powerful model with more functionality at our disposal, even though we might not use all the features all the time. This is not different to using a modern scientific calculator-it has many functions we may not always use, but they are there if we need them. The door to simplification is always open. The converse is not true of the ancient abacus.