Articles | Volume 13, issue 1
Model description paper
29 Jan 2020
Model description paper |  | 29 Jan 2020

The Canadian Hydrological Model (CHM) v1.0: a multi-scale, multi-extent, variable-complexity hydrological model – design and overview

Christopher B. Marsh, John W. Pomeroy, and Howard S. Wheater

Despite debate in the rainfall–runoff hydrology literature about the merits of physics-based and spatially distributed models, substantial work in cold-region hydrology has shown improved predictive capacity by including physics-based process representations, relatively high-resolution semi-distributed and fully distributed discretizations, and the use of physically identifiable parameters that require limited calibration. While there is increasing motivation for modelling at hyper-resolution (< 1 km) and snowdrift-resolving scales ( 1 to 100 m), the capabilities of existing cold-region hydrological models are computationally limited at these scales.

Here, a new distributed model, the Canadian Hydrological Model (CHM), is presented. Although designed to be applied generally, it has a focus for application where cold-region processes play a role in hydrology. Key features include the ability to do the following: capture spatial heterogeneity in the surface discretization in an efficient manner via variable-resolution unstructured meshes; include multiple process representations; change, remove, and decouple hydrological process algorithms; work at both a point and spatially distributed scale; scale to multiple spatial extents and scales; and utilize a variety of forcing fields (boundary and initial conditions). This paper focuses on the overall model philosophy and design, and it provides a number of cold-region-specific features and examples.

1 Introduction

Hydrological models are important tools for understanding past hydrological events, evaluating anthropogenic impacts on natural systems, and informing water resource and management decisions under contemporary and future climates (DeBeer et al., 2015; Freeze and Harlan, 1969; Milly et al., 2008; Mote et al., 2005; Nazemi et al., 2013; Wheater, 2015). Due to the significant role mountains play in the global water supply as “water towers” (Viviroli et al., 2007), the fragility of arctic and mountain ecosystems (Bring et al., 2016), and these regions' sensitivity to anthropogenic climate change (Duarte et al., 2012; Mote et al., 2005; Musselman et al., 2017; Rasouli et al., 2015), there is substantial motivation to provide timely and accurate simulations that can be used to address current and future management challenges in these cold regions. Although the need for multi-scale (Samaniego et al., 2017), hyper-resolution (sub-1 km) (Wood et al., 2011), and snowdrift-resolving scales (1 to 100 m) (Pomeroy and Bernhardt, 2017) is becoming clear, contemporary cold-region models suffer from shortcomings when run over large extents and high spatial resolutions, and they may be limited to the spatial scale at which they operate.

Numerous studies suggest that model performance is greatly improved in cold regions when including explicit spatial heterogeneity, identifiable parameter spaces, and a full range of cold-region hydrological processes, e.g., Pomeroy et al. (1998a, b), Bartelt and Lehning (2002), Bowling et al. (2004), Etchevers et al. (2004), Raderschall et al. (2008), Dornes et al. (2008b), Essery et al. (2013, 2009), Pomeroy et al. (2013), Fang et al. (2013), Fiddes and Gruber (2014), Kumar et al. (2013), Endrizzi et al. (2014), Mosier et al. (2016), and Painter et al. (2016). A better understanding of the physical system instead of solely focusing on parameter optimization (Bahremand, 2016) and ensuring that models are not needlessly constrained by the rigidity of the model structure, choice of parametrization, and representation of spatial variability and hydrological connectivity (Mendoza et al., 2015) are expected to further predictive capacity. Physics-based models may also limit the reliance upon calibrated effective values and decrease uncertainty due to requiring the use of physically identifiable parameters (Fatichi et al., 2016; Pomeroy et al., 2013). The use of lightly calibrated or uncalibrated models is increasingly important for simulating future conditions, as climate non-stationarity increases the uncertainty of calibrated models (Brigode et al., 2013; Vaze et al., 2010). Distributed, physics-based models are thus often the most appropriate type of hydrological model for simulating distributed state variables (Dornes et al., 2008a; Fatichi et al., 2016), simulating catchments with extreme heterogeneity (Kumar et al., 2013), or simulating process interactions (Dornes et al., 2008a; Horne and Kavvas, 1997; Maxwell and Kollet, 2008). These improvements motivate the continued development of spatially discrete, physics-based models.

Although there is uncertainty in the optimum levels of complexity required in cold-region hydrological models (Avanzi et al., 2016; Clark et al., 2017), such models have unique requirements and considerations; a brief summary follows. The largest discharge event of the year often results from the melt of the seasonal snowpack (Davies et al., 1987; Gray and Male, 1981), and therefore substantial effort has been invested in snow model development, e.g., Jordan (1991), Marks et al. (1998), Bartelt and Lehning (2002), Vionnet et al. (2012), Leroux and Pomeroy (2017), as well as flexible snow-cover modelling systems, e.g., the Factorial Snow Model (FSM) (Essery, 2015) and ES-CROC (Ensemble System Crocus) (Lafaysse et al., 2017). Streamflow discharge is impacted by snowmelt spatial heterogeneity due to variability in surface energetics (Carey and Woo, 1998; Dozier and Frew, 1990; Harder et al., 2019; Marks et al., 1992; Mott et al., 2013; Munro and Young, 1982; Olyphant, 1986; Pomeroy et al., 2003; Schlögl et al., 2018), precipitation spatial variability (Harder and Pomeroy, 2013; Lehning et al., 2008; Marks et al., 2013), vegetation canopy interception (Hedstrom and Pomeroy, 1998; Kuchment and Gelfan, 2004), and snow redistribution via wind processes (Essery et al., 1999; MacDonald et al., 2009; Mott et al., 2010; Pomeroy et al., 1993; Pomeroy and Li, 2000; Winstral et al., 2002). Snowmelt runoff is further complicated due to frozen soils that limit infiltration rates (McCauley et al., 2002; Zhao and Gray, 1999) such that standard infiltration representations are insufficient (Lundberg et al., 2016). Active layer depth above permafrost dramatically impacts surface characteristics (e.g., topography, vegetation, soils), streamflow seasonality, and water partitioning (Walvoord and Kurylyk, 2016). In cold regions, the numerous lakes and wetlands impact the local climate during ice-free periods (Latifovic and Pouliot, 2007; Rouse et al., 2005; Shook et al., 2015). In summary, cold-region hydrological models have unique challenges and must include a variety of process representations not considered in most temperate hydrological models.

There are, however, significant limitations in hydrological modelling ability. For instance, there are deficiencies due to substantial heterogeneity and difficulty in observing surface and subsurface parameters and processes (Freeze, 1974), no single scale at which the homogeneity of control volumes is achieved (Beven, 1989; Blöschl and Sivapalan, 1995; Klemeš, 1983; Shook and Gray, 1996), and mismatches between underlying theory and applied scales (Or et al., 2015). These limitations manifest as (1) uncertainties in model parameters, initial conditions, boundary conditions, forcing data; (2) incomplete process representations, selections, and linkages (Beven, 1993; Beven and Westerberg, 2011; Clark et al., 2008; Fatichi et al., 2016; Raleigh et al., 2015; Slater et al., 2013; Wagener and Montanari, 2011); and (3) issues of complexity including the degree of physics-based equations, the number of parameters, forcing data requirements, and spatial discretization requirements (Beven, 1993; Clark et al., 2008; Hrachowitz and Clark, 2017). In addition, these models may have limited structural flexibility for incorporating multiple modelling philosophies (e.g., Dornes et al., 2008b; Clark et al., 2011) or have limitations in incorporating next-generation data products such as unmanned aerial vehicle (UAV) imagery (Bühler et al., 2016; Harder et al., 2016; Spence and Mengistu, 2016). Without care, physically based, mechanistic approaches can result in over-parameterized models (Perrin et al., 2001) that are highly uncertain and difficult to verify (Beven, 1993) due to a mismatch in model element and observed scales as well as limited high-resolution spatially distributed data (Beven, 1989). Physically based models should be used critically, with a proper appreciation of the strengths and limitations, and should be dependent on the purpose of the modelling (Beven, 1993, 2006; Das et al., 2008; Perrin et al., 2001).

In order to address the scientific and societal demands placed on hydrological models, there is a need for a new generation of hydrological models that allow for the following.

  1. Multi-scale, spatially distributed process representation. Although semi-distributed schemes such as the group response unit (GRU) or hydrological response unit (HRU) approach have had substantial success in cold regions, e.g., Pietroniro et al. (2007), Pomeroy et al. (2007), and Clark et al. (2015), complex spatial behaviours cannot be modelled unless the HRUs are constructed a priori to produce the behaviours. This limits simulating cascading processes and emergent behaviours, e.g., the accumulation of non-linear process interactions leading to basin-wide behaviours. Representing mass and energy heterogeneities as well as interactions at multiple spatial scales (Hrachowitz and Clark, 2017; Samaniego et al., 2017) moving towards regional predictions (Sivapalan, 2018) has been suggested as a path to improving predictive capacity. Fully distributed, raster-based models are inefficient with the need for many raster cells, greatly limiting the applicability for both high resolution and over large extents. The deficiencies in HRU, GRU, and raster-based models point towards a need for an improved terrain representation that allows for both high resolution as needed and applicability for modelling over large extents.

  2. Flexible model structure. Many models use a rigid model structure that does not allow for easily changing model algorithms and all parameters or easily testing different algorithms or hypotheses. An improved approach is to allow process modularity for easily modifying aspects of a model's structure and complexity. Such model flexibility has been present in many rainfall–runoff models, e.g., MMS (Leavesley et al., 2002), FUSE (Clark et al., 2008), and SUPERFLEX (Fenicia et al., 2011), but to the authors' knowledge, such modularity in cold-region models has been limited to the Cold Regions Hydrological Model (CRHM) (Pomeroy et al., 2007) and SUMMA (Clark et al., 2015). These are both modular, physics-based, semi-distributed hydrological response unit (HRU) models with capability for cold-region hydrology. A flexible model structure should allow for easily scaling between temporal scales (i.e., time stepping), spatial extents, spatial resolutions, and process representations as required. Assumptions on explicit coupling between processes leads to difficulty in testing different process representations and limits the inclusion of existing code. Despite the rich set of cold-region snow and hydrological models, e.g., Alpine3-D (Lehning et al., 2006), iSnobal (Marks et al., 1998), GeoTOP (Endrizzi et al., 2014), MESH (Pietroniro et al., 2007), CRHM (Pomeroy et al., 2007), SUMMA (Clark et al., 2015), SRGM (Gelfan et al., 2004; Kuchment and Gelfan, 2004), ESCROC (Lafaysse et al., 2017), and VIC (Cherkauer et al., 2003), there are no explicitly distributed, modular cold-region models.

  3. Ease of changing model parameters, as well as initial and boundary conditions. Model parameters, initial conditions, and boundary conditions are uncertain in hydrological systems and are a significant constraint on model complexity and validity. Hard-coded parameters can be a significant source of uncertainty as they are effectively treated as physical constants (Mendoza et al., 2015). Modern models must be developed so that changing initial conditions, parameters, and all aspects of the model configuration are trivial and easily done within the context of an uncertainty framework. Due to the long temporal durations for which climate change scenarios are done, flexibility in changing surface parameters with time, e.g., vegetation cover, also needs to be possible.

  4. Efficient use of computational resources. Unlike GRU- or HRU-based models, distributed models are generally discretized using a raster approach with a fixed spatial resolution. This can lead to either increased computational requirements or non-optimum use of computer resources due to the overrepresentation of the surface (e.g., homogenous locations), while choosing a coarser-sized mesh may result in failure to capture quickly varying and extremely important heterogeneity. Because of the general overrepresentation of topography via a fixed-resolution raster, these distributed models become difficult to parametrize and computationally expensive to run, limiting their applicability to large spatial extents. Using more efficient terrain representations as well as modern high-performance computing paradigms can reduce this wasted computational effort.

  5. Allow appropriate model complexity. Raster-based models with high-resolution grid cells, and wasted computational effort as noted above, often lead to arbitrary complexity reduction and process removal due to computational constraints. Reducing the model runtime is often a justification for simpler conceptual models, for simpler landscape representations, and for fewer computational elements. Hydrological model complexity should be warranted based upon the simulation results and needs but not for simplicity's sake.

Although there are substantial advantages to using physics-based, fully distributed models, data (forcing and validation) and computational limitations have slowed their development and adoption. However, recent technological progress has been progressively removing some of these limitations. For example, unmanned aerial vehicle (UAV) imagery is providing submetre digital surface and elevation maps (Bühler et al., 2016; Harder et al., 2016), vegetation classification (Spence and Mengistu, 2016), hydrological features (Spence and Mengistu, 2016), and initial conditions, e.g., snow cover (Bühler et al., 2016; Harder et al., 2016). Surface geophysical methods are improving the characterization of large-scale subsurface properties (Hubbard et al., 2013). Remote sensing products of soil properties are of increasingly higher quality (Mohanty, 2013), and high-resolution satellite imagery can be used to diagnose spatial patterns of snow cover (Wayand et al., 2018). Widespread access to high-performance computing (HPC) resources, e.g., Compute Canada (Canada), Extreme Science and Engineering Discovery Environment (XSEDE; United States), National Computational Infrastructure (NCI; Australia), and the Horizon 2020 initiative (European Union), can help offset the increased computational cost of the simulations and of the uncertainty analysis needed to constrain a priori estimated physically based parameters (Paniconi and Putti, 2015). Lastly, efficient uncertainty analysis frameworks such as VARS (Razavi and Gupta, 2016) can also decrease the total number of required simulations to estimate uncertainty, further reducing the computational burden. However, estimates of critical subsurface properties such as hydraulic conductivity cannot be represented a priori with sufficient confidence or at the correct scale, viz. effective model element parameters (Binley et al., 1989), to avoid calibration (Freeze, 1974).

In summary, models will always require a trade-off between computational complexity (e.g., algorithms, landscape representation, initial conditions, parameters, and terrain discretization) and model performance (e.g., modelled versus observed). Cold-region hydrological models have unique requirements that motivate the inclusion of explicit spatial heterogeneity via semi-distributed and fully distributed discretizations. To simulate the complex inter-process interactions that lead to important hydrological features, a variety of features must exist within a distributed, process-based modelling framework.

This paper outlines the philosophy and details of a new hydrological model, the Canadian Hydrological Model (CHM), and how the development of this modelling framework addresses the above outlined limitations of many existing hydrological models and contributes to cold-region modelling. This paper focuses on the overall model philosophy and design, and it provides a small number of cold-region-specific features and examples.

2 Design and overview

2.1 Overview

The Canadian Hydrological Model (CHM) is a spatially distributed, modular modelling framework. Although not restricted to cold regions, it is designed with both cold-region and temperate zone processes in mind and has various capabilities that facilitate the modelling of these domains. The design goal of CHM is to use existing high-quality open-source libraries and modern high-performance computing (HPC) paradigms. By providing a framework that allows for as loose or tight a coupling between processes as required, CHM allows for the integration of current state-of-the-art process representations and makes no assumptions about the complexity of these process representations. For example, it allows for the testing of the representations in a consistent manner, diagnosing model behaviour due to parameter changes, process representation changes, and basin discretization. Spatially, it allows for domains at point (10−6 km2), hillslope (1 to 10 km2), basin (100 km2), regional (8000 km2), and provincial or state (> 1 000 000 km2) scales. The following sections outline the framework features, including terrain representation, surface parameterization, process representation, meteorological inputs, parallelism, uncertainty analysis, visualization and analysis, and adaptation of raster algorithms. Although the CHM will eventually include the entirety of the hydrological cycle, at this time only snow accumulation and surface meteorology processes are implemented. Additional model components are being developed and will be available in future versions of CHM.

2.2 Terrain representation

The spatial variability of terrain is a key component to any model and is an important component of model complexity. Regardless of how sophisticated, physics-based, and spatially explicit a hydrological model may be, at some level the hydrological system is conceptualized and aggregated into a control volume (Vrugt et al., 2008). Structured meshes, also known as rasters and grids, are a landscape discretization whereby the landscape is discretized by uniformly sized cells. Raster-based hydrological models are common (Tucker, 2001) because their computer representation is trivial, and the widespread use of rasters, such as in remotely sensed data, makes using them a natural choice in hydrological models. However, rasters have a number of significant limitations, the most limiting being a fixed spatial resolution over the entire basin (Tucker, 2001). This results in potentially large computational inefficiencies due to the overrepresentation of topography. This arises as a result of requiring small raster cells (elements) to capture the spatial variability in areas of high topographic variability or (sub)surface variability (e.g., vegetation, soils), which results in the overrepresentation of areas that have limited spatial variability. Coarse-resolution rasters also have discontinuities in the elevation data, and adjacent cells may have large elevation differences.

Figure 1Example of a variable-resolution triangulation mesh as produced by Mesher for the Bow River Basin west of Calgary in the Canadian Rocky Mountains and foothills. The triangular edges are shown as grey lines overlain on the original digital elevation model (DEM).

Unstructured triangular meshes, sometimes referred to as triangulated irregular networks (TINs), represent the topography via a set of irregularly sized, non-overlapping connected triangles, with each triangle face being of a constant slope (Chang, 2008). Areas of large topographic variability can have a higher density of small triangles in order to capture the spatial variability, and areas of relatively homogeneous topography have fewer large triangles. This a more efficient terrain representation than rasters (Shewchuk, 1996) and may have up to a 90 % reduction in computation elements (Ivanov et al., 2004; Marsh et al., 2018). Despite these computational advantages, a practical downside is that due to the widespread availability of raster data, conversion to an unstructured mesh is required. This results in increased uncertainty due to aggregation of the landscape into control volumes. The CHM uses a novel multi-objective approach for unstructured triangular mesh generation, Mesher, detailed in Marsh et al. (2018). A brief summary follows: quality Delaunay meshes are generated, ensuring a smooth graduation between small and large triangles; triangles are bounded with minimum and maximum triangle areas to ensure that process representations match the physical scale; and triangles are generated to fulfil tolerances (e.g., root mean square error – RMSE) to the underlying topographic raster and other important landscape features such as vegetation and soils. This mesh generation attempts to limit the amount of error introduced by the approximating surface given by the unstructured mesh and provides mechanisms to ensure that spatial heterogeneity in the landscape is correctly preserved.

Using this mesh generation, simulation domains can be constructed at a variety of spatial extents and, importantly, spatial scales. An example of this variable-resolution triangulation mesh for a part of the Bow River Basin in the Canadian Rockies and foothills west of Calgary, Alberta, Canada, is shown in Fig. 1. The triangular edges are shown in grey lines. The variable resolution produces larger triangles in the valley bottoms, where topographic variability is limited, and small triangles in the mountains, where the heterogeneity is greater. This allows for diagnosing the impact of scale on model performance as well as matching the process representation to the correct model length scale. Further constraints could ensure that streams are accurately defined.

Figure 2Directed acyclic graph showing module dependencies. Lines point to the module that requires the listed dependency. In this example, a snow-cover model, Snobal, is being driven by meteorology in order to drive a frozen soil infiltration model (Gray_inf).


2.3 Triangle parameterization

Setting values of parameters for the triangles, such as assigning vegetation or soil type to the triangle, is done during the mesh generation phase. The parameter values are stored in a file separate from the underlying mesh and can thus be easily changed at runtime. This allows for easily investigating the impact of parameter values on outputs. The parameterization of the triangles is done by (a) determining the valid raster cells under each triangle, (b) calculating an error metric for these cells, and assigning this value to the triangle. Maximum and mean are the two most commonly used methods, but it can be any user-defined function. For classified data, the mode is used. This would allow, for example, for the selection of the most dominant land-cover class. In addition, a user-specified classifier function can be given to easily classify continuous input parameters; e.g., classifying vegetation heights into vegetation classes. Lastly, CHM provides mechanisms to write model output to a format that can be used as input; that is, CHM can use its output to set triangle values for future simulations.

2.4 Modular process representation structure

A hydrological model is a hypothesis based on assumptions of how a hydrological system works (Savenije, 2009). Modular model structures allow for rigorously testing process representations and have been used with success in cold-region hydrology, e.g., the Cold Regions Hydrological Model (CRHM) (Pomeroy et al., 2007) and Structure for Unifying Multiple Modeling Alternatives (SUMMA) (Clark et al., 2015). A feature of CHM is that it provides a modular process representation that is suitable for distributed modelling, while maintaining high computational performance and flexibility.

In CHM, process representations are conceptualized into modules. Selecting various combinations of these modules in the CHM framework defines the overall model. A principal design goal of the module system is that a module has an enforced set of preconditions and post-conditions. Preconditions represent the variables that must be computed prior to a given module running, and post-conditions encapsulate variables that must be computed by the currently running module so as to be available as input for other modules. At runtime, the user-selected modules are linked together into a directed acyclic graph based on these variable dependencies, and module execution order is determined via a topological sort of this graph. This sort ensures that modules are run in an order so as to fulfil the precondition (i.e., the variable dependencies). Linkages between modules showing these dependencies are shown in Fig. 2. The lines with arrows show how variable dependencies are resolved between modules. The lines going from a module are the post-conditions that satisfy the preconditions of the next-to-be-run module. In this example, a snow-cover model, Snobal, is being driven by meteorology, with the output of Snobal being used as input to a frozen soil infiltration model (Gray_inf).

The hydrological literature has a diverse set of process representations that are either one-dimensional with no lateral exchange between elements (point scale) or are explicitly coupled with surrounding elements (Todini, 1988). CHM makes no assumption about either, and modules may either operate on a single triangle or on the entire domain. If only point-scale modules are selected, then CHM may be optionally run at a point scale, effectively disabling the rest of the distributed framework. As there are substantial merits to mixing top-down and bottom-up process representations (Hrachowitz and Clark, 2017; Pomeroy et al., 2004), CHM makes no assumptions on the complexity or type of process representation in a module – modules may be a mix of complex physics-based representations and conceptual representations. This also applies to process coupling. For example, a module could be a single process (e.g., a snow model), a coupled set of processes (e.g., coupled heat and energy snow model + frozen soil routine), or an entire existing model.

Table 1Cold-region surface process representations currently available in CHM.

Download Print Version | Download XLSX

Due to the strict preconditions and post-conditions required for module dependency resolution and the abstraction used in CHM, existing libraries and code can be used in a model. There is no need to rewrite the code. Therefore, any code that may be called via a C interface (e.g., Fortran, R, Python, MATLAB) is suitable to be used (with a few considerations) as a CHM module.

Summarized in Table 1 and described in brief below is a list of the processes currently available in CHM. Two energy balance snowpack models are available, Snobal and SNOWPACK. Snobal (Marks et al., 1999) is a two-layer energy balance model with a fixed upper layer that is used for the estimation of outgoing longwave radiation and atmosphere–snow temperature gradients for the turbulent heat flux. SNOWPACK (Bartelt and Lehning, 2002) is a multilayer finite-element energy balance snow model originally developed for avalanche hazard forecasting. In addition to the snow-cover albedo estimates provided by Snobal and SNOWPACK, the Canadian Land Surface Scheme (CLASS) albedo routine is available (Verseghy, 1991). Frozen soil infiltration is calculated using the parametric form of Gray et al. (2001). Horizontal snow mass is redistributed using a 3-D advection–diffusion blowing snow model derived for unstructured meshes (Marsh et al., 2020). Blowing snow saltation (Pomeroy and Gray, 1990), turbulent suspension (Pomeroy and Male, 1992), sublimation (Pomeroy et al., 1993), threshold shear stress for saltation (Li and Pomeroy, 1997), shear stress partitioning by vegetation and snow, and probabilistic upscaling (Pomeroy and Li, 2000) parameterizations comprise the blowing snow model. The vertical redistribution of mass in steeply sloping terrain is calculated using Snowslide (Bernhardt and Schulz, 2010) with a threshold slope and mass exceedance to transport mass downslope (i.e., it is not a prognostic avalanche model). The forest canopy is conceptualized into open and forest areas and uses the snow interception algorithm of Hedstrom and Pomeroy (1998) coupled to the intercepted snow sublimation and unloading algorithms of Pomeroy et al. (1998b) as well as the drip and rapid unloading formulations in Ellis et al. (2010). Subcanopy shortwave and longwave irradiance and turbulent transfer algorithms from Ellis and Pomeroy (2007), Pomeroy et al. (2009), and Ellis et al. (2010) are also included in the canopy module.

2.5 Input meteorology

Input meteorology is prescribed as a point source (herein, “virtual station”) defined by latitude, longitude, and elevation. However, a virtual station may have an arbitrary location and elevation, and it need not be within the simulation domain or correspond to a real meteorological station. This allows a virtual station to be located at, for example, the centroid of a numerical weather prediction output grid cell. Because all input meteorology is given as a point source, various spatial interpolants are present in CHM to provide a distributed field across all triangles.

Spatial interpolates are present as inverse distance weighting (IDW) and thin plate spine with tension. In some cases, no interpolation is desired, and therefore a third option called “nearest” is available – this uses the nearest virtual station without any spatial interpolation. Over large domains, such as when using numerical weather prediction output, every virtual station in the simulation domain should not be used in the interpolation to every triangle. Therefore, interpolants may query a list of either (a) virtual stations within some distance of the triangle or (b) the closest n virtual stations. This ensures that only nearby virtual stations are used to form the interpolant. Vertical elevation correction is provided by a set of specialty modules. All virtual stations are corrected to a common reference level using these modules prior to spatial interpolation. A list of these algorithms is summarized in Table 2.

Input meteorology may be given as either text files or as NetCDF files (Rew and Davis, 1990). When NetCDF files are used, the time step data are lazy-loaded such that only the current time step is read. This decreases the up-front load time and decreases total memory usage.

Table 2List of available meteorology interpolants.

Download Print Version | Download XLSX

2.6 Input filters

Input filters provide a mechanism to modify input meteorology during runtime. This is similar to the filter feature in CRHM (Pomeroy et al., 2007) and MeteoIO (Bavay and Egger, 2014). Filters are assigned to each virtual station, and each virtual station may have an arbitrary number of filters. The purpose of filters is to allow, for example, values outside a certain range to be filtered or to perform a correction such as taking an observed wind speed at 2 m and changing it to 10 m for use later in a process module. Filters operate per time step and can therefore consider the previous model time step for use in the correction; e.g., including snow depth to perform vertical wind speed height correction.

2.7 Point mode

Due to the difficulty in validating spatial models due to limited spatial observations, evaluation is generally performed using point observations. CHM may be run in point mode, which allows for simulating a single triangle without lateral interactions, using a specialized input module to pass a hydrometeorological station's observation data directly to the underlying process models. This is intended to simulate a point collocated with an input observation meteorology dataset and allows for traditional point simulations.

2.8 High-performance computing

In CHM, parallelism is currently implemented via the shared memory OpenMP library. Coding a process representation into a module will generally result in either a point-scale module (e.g., point-scale snow-cover model) or it will be a spatially coupled model (coupled advection–diffusion equation). The first type, owing to the fact it does not require knowledge of its neighbours to compute a value, corresponds to an embarrassingly parallel problem – that is, a problem that does not require any communication between threads. Herein, these are referred to as data parallel. Spatially coupled models require the solution at their neighbour triangles in order to compute a solution. These neighbours, in turn, require solutions at their neighbours, and so on. Therefore, this is a much more challenging type of problem to introduce parallelism to. Herein, these are referred to as domain parallel. Data parallel modules automatically have the parallelism implemented and require no special consideration from the developer. Domain parallel modules, however, require the module developer to implement parallelism as appropriate for the module.

Mixing these two types of parallelism complicates the implementation of parallel code. To provide as much seamless parallelism as possible, each module declares the type of algorithm it is: data parallel or domain parallel. After the topological sort is performed to determine module execution order, the modules are scheduled together into groups that share a parallelism type. For example, consider the following sorted list of modules, with their parallelism type in brackets.

  • mod_A (parallel::data)

  • mod_B (parallel::data)

  • mod_C (parallel::data)

  • mod_D (parallel::domain)

  • mod_E (parallel::data)

These would then be scheduled together into three groups.

  • Group 1

    • mod_A (parallel::data)

    • mod_B (parallel::data)

    • mod_C (parallel::data)

  • Group 2

    • mod_D (parallel::domain)

  • Group 3

    • mod_E (parallel::data)

Figure 3Output from CHM is in the ParaView format, allowing for time series analysis and full 3-D visualization in ParaView. Shown is shadowing over Marmot Creek Research Basin, Alberta, Canada.

The modules in group 1 are run in parallel together. Because they are data parallel, only one iteration over the mesh is required. Then, groups 2 and 3 are run. This scheduling mechanism reduces the overhead of a modular approach by limiting total iterations over the mesh and minimizing thread creation. Further, as most hydrological process representations are point scale, it allows for abstracting parallelism, resulting in “free” parallelism for the developer.

2.9 Uncertainty analysis

CHM provides a mechanism to easily allow modules to obtain parameter values from configuration files (JSON format), overriding the default hard-coded value. Changes to the model structure (i.e., choosing modules), initial conditions, and parameter files (e.g., land cover) are also done via this mechanism. Users may, via the command line, change any configuration value – thus simplifying uncertainty testing. This mechanism reduces situations in which changes require recompilation.

The Python code snippet shown in Listing 1 demonstrates changing values on the command line (via Python). This code is setting the name of three output files and adding a new module to be run.

Listing 1Example for setting output file names and adding a new module.

Figure 4Marmot Creek Research Basin, Kananaskis Valley, Alberta, in the Canadian Rocky Mountains. The basin outline is given as solid black, 100 m contour lines are shown in brown, stream channels are shown in blue, and man-made clearings are shown as hatched areas. The meteorological stations used for this study are shown as crosses. The southernmost set of clearings are ski runs in the Nakiska Ski Resort.

2.10 Visualization and analysis

The output format used is the ParaView (Ahrens et al., 2005) unstructured mesh format. This allows for the visualization of the simulation results in full 3-D, with time series analysis in ParaView, as shown in Fig. 3. The addition of a ParaView plug-in for CHM allows for displaying the date and time of the output. The animation view allows for exploring the spatio-temporal results. It also allows for immediate diagnosis of modelling errors, especially if the spatial pattern of an output variable is clearly incorrect: for example, if a coding error resulted in a patchwork of air temperatures instead of an expectedly smooth gradient with elevation, snowdrifts being formed in locations that were known to be incorrect such as the top of a ridge instead of in the lee, or Northern Hemisphere north-facing slopes receiving the most shortwave irradiance. There are many post-processing filters and tools available in ParaView, such as plotting an individual triangle's values over time. Because ParaView uses the Visualization Toolkit (VTK) library (Schroeder et al., 2006), the ParaView files can easily be loaded and post-processed using the Python VTK library in conjunction with traditional Python libraries such as NumPy (Oliphant, 2006) and SciPy (Jones et al., 2018).

In addition to the ParaView output, CHM provides a set of post-processing scripts that allows for converting the ParaView file to a rasterized GeoTiff or NetCDF file. This allows for using the output in post-processing algorithms that require arrays, or in GIS.

2.10.1 Adaptation of raster-based algorithms

The adaptation of raster-based algorithms is an important aspect of CHM as many existing algorithms are raster-based. Frequently, raster-based algorithms employ logic that performs queries such “look X length units in direction Y”. This is easily done on a structured mesh; however, on an unstructured grid, this process is non-obvious. Iterating over each triangle's neighbours results in a random walk across the domain, and brute-force iteration search methods are needlessly slow. CHM uses the k-d spatial search tree available within the dD Spatial Searching (Tangelder and Fabri, 2018) package in the Computational Geometry Algorithms Library (CGAL) to optimize spatial queries. Briefly, a k-d tree is a generalization of a binary search tree in high dimensions that decomposes the search domain into a set of small subdomains (Bentley, 1975). This tree structure can then be recursively searched, resulting in efficient spatial lookups. The k-d tree implementation is how nearby stations are determined. This technique for spatial searching can also be used to calculate terrain parameters, such as the terrain curvature.

3 Model application

3.1 Overview

The following section describes the methodology for evaluating various features of CHM and provides examples of usage. Although the CHM will eventually include the entirety of the hydrological cycle, snow accumulation and surface meteorology processes are currently implemented. Marmot Creek Research Basin (MCRB) in the Canadian Rockies in Alberta, Canada, is used as a location to test the two snow modules and various models to provide the driving meteorological forcing for these models that are currently implemented in CHM. The meteorological interpolants are tested in a leave-one-out validation across the MCRB. In addition, an adaptation of a raster-based terrain shadowing for shortwave irradiance calculation is presented, demonstrating the conversion of an algorithm from a raster to unstructured mesh. Finally, the parallel computation aspect of CHM is tested by performing a scaling analysis using different numbers of CPUs.

3.2 Study site

3.2.1 Marmot Creek

Marmot Creek Research Basin (MCRB) (Golding, 1970) is located in the Kananaskis River Valley of the Canadian Rockies, as shown in Fig. 4. It is a 9.4 km2 basin covered predominately by needle-leaf forest (Fang et al., 2013; Pomeroy et al., 2012). The climate is dominated by continental air masses with long and cold winters; however, these are interrupted by frequent chinooks (Foehns) in midwinter (DeBeer and Pomeroy, 2009). It spans an elevation range from 1700 to 2886 m (Rothwell et al., 2016), and snow covers the upper elevations of the basin from October to June. The average seasonal precipitation is approximately 600 mm at low elevations, increasing to over 1140 mm at the treeline (Rothwell et al., 2016).

3.2.2 Meteorological observations

Meteorological observations for air temperature, relative humidity, wind speed, precipitation, soil temperature, and incoming shortwave radiation for the Upper Clearing (1860 m), Vista View (1956 m), and Fisera Ridge (2325 m) sites, shown as crosses in Fig. 4, were used. Gap-filled, quality-corrected 15 min data for the water years 2007 to 2016 (inclusive) were used. Please see Fang et al. (2019) for further details. Precipitation was measured with Alter-shielded Geonor weighing precipitation gauges and corrected for wind-induced under-catch (Smith, 2009). Precipitation phase was determined via the psychrometric energy balance method of Harder and Pomeroy (2013). Longwave irradiance was calculated following Sicart et al. (2006). This was developed for mountainous terrain and was shown to have an error of less than 10 % over the snowmelt season. This method has been used with success at the MCRB.

Periodic snow surveys of depth and snow water equivalent (SWE) on long transects at Upper Clearing were conducted by various members of the Centre for Hydrology and used to quantify snowpack density. For each transect, there were at least 25 snow depth measurements and at least 6 gravimetric snow density measurements using an ESC-30 snow tube (Fang et al., 2019).

3.3 Models

3.3.1 Snow models

Point-scale evaluation of the two snow models in CHM, Snobal and SNOWPACK, was done at the Upper Clearing site.

Snobal (Marks et al., 1999) is a physics-based, two-layer snowpack model designed specifically for deep mountain snowpacks and approximates the snowpack with two layers; the surface fixed-thickness active layer (taken here as 0.1 m) is used to estimate surface temperature for outgoing longwave radiation and the atmosphere–snow exchange of sensible and latent heat via turbulent transfer. Snobal features coupled energy and mass balance, internal energy tracking, and liquid water storage calculations. Turbulent fluxes are explicitly calculated via Marks et al. (1992), a bulk transfer approach that includes a Monin–Obukhov stability correction. The ground heat flux is calculated from conduction with a single soil layer of known temperature.

SNOWPACK (Bartelt and Lehning, 2002; Lehning et al., 2002) is a multilayer finite-element model of mountain snowpacks, with application for avalanche hazard forecasting. It describes the microphysical properties of a snowpack and includes the dynamic addition and/or removal of snow layers using a system of partial differential equations (PDEs). These are discretized vertically into an arbitrary number of snow layers in a Lagrangian coordinate system. It has coupled energy and mass balance, internal energy, and liquid water storage calculations with a bulk transfer turbulent flux scheme with Monin–Obukhov stability correction (Michlmayr et al., 2008). The default Michlmayr et al. (2008) scheme was used herein.

Both SNOWPACK and Snobal were configured to use the albedo routine of Verseghy et al. (1993). The snow models were driven with observed precipitation, shortwave irradiance, wind speed, air temperature, relative humidity, and soil temperature at a 15 min time interval. Because of the sheltered nature of the Upper Clearing, no blowing snow was simulated (Musselman et al., 2015). Snow model parameters, such as roughness length, were set following Pomeroy et al. (2012).

3.3.2 Mesh generation

The unstructured mesh was created using the Mesher software. A 1 m × 1 m input elevation lidar DEM (Hopkinson et al., 2011) was used. The resulting mesh was generated to have a minimum triangle area equivalent to a 25 m × 25 m raster and represented the topography to within 25 m RMSE. This resulted in approximately  100 000 triangles.

Figure 5Dozier and Frew (1990) horizon shadowing algorithm. For observer A, a search along the azimuth that corresponds to the solar vector S is performed such that if the slope of H is greater than that of S, A is in shadow.

3.3.3 Raster algorithm adaptation (shadowing)

An example of the adaptation of a raster algorithm to the unstructured mesh is shown for a terrain shadowing algorithm, illustrated in Fig. 5, that calculates the shadows cast from surrounding terrain. The “look X length units in Y direction” query is required for finding obstructing terrain (e.g., a tall mountain) by searching along the azimuthal direction towards the sun. As a demonstration of the k-d tree usage in CHM, the shadowing algorithm of Dozier and Frew (1990) (herein DF90) was implemented for unstructured triangular meshes. In brief, the DF90 algorithm searches along an azimuthal direction within some horizontal distance and attempts to find terrain that is above the solar elevation. For an observer A, a search along the azimuth that corresponds to the solar vector S is performed. For each terrain element found, a new vector (H) is calculated. If the slope of H is greater than that of S, A is in shadow. Terrain is searched from the observer towards some maximum search radius in steps of size dx. Specifically, this adaptation of DF90 required using the k-d tree to find the triangle at a distance X m from the source triangle (A) along an azimuth that corresponded to the solar vector, S.

The DF90 shadowing algorithm was run for all of Marmot Creek using the mesh described in Sect. 3.2. A maximum search radius of 1000 m was used, discretized into 10 steps. The guidelines for choosing these search values follow two criteria: (1) the radius should be large enough to cover the distance across a representative valley length distance such that shadows from mountains across the valley are included; and (2) the step should be about half of a triangle length scale such that steps do not pass over triangles. The DF90 implementation was compared to observed shadowed area (see below), the Marsh et al. (2012) shadowing model, and the Solar Analyst (Fu and Rich, 1999) shadow model. Solar Analyst is an extension in the ArcGIS software by the Environmental Systems Research Institute (ESRI). The observed shadowed areas are from time series images from the field campaign detailed in Marsh et al. (2012) and were orthorectified using the software of Corripio (2004). Shadow location for 1 February 2011 at 17:00 was used in this comparison. The output from CHM was rasterized from the unstructured mesh at a 1 m × 1 m spatial resolution.

3.4 Leave-one-out comparison

To test the efficacy of the meteorological interpolants, a leave-one-out comparison was conducted for the Upper Clearing, Vista View, and Fisera Ridge stations. This entailed using two of the three meteorological stations as input for CHM in order to predict the third. For example, Upper Clearing and Vista View were used to predict meteorological conditions at Fisera Ridge, and Vista View and Fisera Ridge were used to predict Upper Clearing.

A total of 10 water years using 15 min data were simulated. The following meteorological interpolants were used: terrain shadowing (Dozier and Frew, 1990), cloud fraction (Walcek, 1994), air temperature (Cullen and Marshall, 2011), relative humidity (Kunkel, 1989), precipitation phase (Harder and Pomeroy, 2013), precipitation (Thornton et al., 1997), and solar radiation transmittance estimated from observed incoming shortwave values.

3.5 Parallel scaling

The heterogenous Westgrid cluster Graham was used to investigate the scaling performance of the CHM code with various numbers of CPUs. The base nodes were used. These have two Intel E5-2683 v4 Broadwell CPUs at 2.1 Ghz for a total of 32 cores and 128 GB of RAM. The modules run include the data parallel Snobal snowpack module, as well as a domain parallel advection–diffusion blowing snow module (Marsh et al., 2020).

Simulations were run for a mesh with  100 000 triangles. The model was run with 1, 2, 4, 6, 8, 16, and 32 cores, and for each core-count scenario, the fastest of five runs was taken. File output was disabled for these runs. The speed-up for the n core run (coren) was computed relative to the one-core run (core1):

(1) speed-up = core 1 core n .

Figure 6Comparison of Snobal (red) and SNOWPACK (blue) run as a point simulation within CHM for the Upper Clearing site at Marmot Creek Research Basin for 10 water years. Manual snow course observations are shown as black dots.


Table 3Root mean square error (RMSE, mm) and mean bias (MB, mm) for the SNOWPACK and Snobal models at the Upper Clearing site for each water year.

Download Print Version | Download XLSX

4 Results

4.1 Point-scale snow model

Shown in Fig. 6 is the simulated snow water equivalent (SWE, mm) for SNOWPACK (blue) and Snobal (red). The water year is denoted above each plot. Snow course observations are shown as black dots. The RMSE and mean bias (MB) values for both models for each water year are shown in Table 3 and averaged over all years in Table 4.

In 2007, SNOWPACK overestimates peak SWE more than Snobal, although ablation timing between the two is identical. In 2008, early season SWE is overestimated by SNOWPACK, although late season SWE is better estimated by SNOWPACK. Water year 2009 is poorly simulated in general, especially by Snobal. It is not clear what causes this poor performance. During the cold winters of 2010 and 2011, both models perform well. In 2012, Snobal underestimates peak SWE versus SNOWPACK. For the years 2013 to 2015 SNOWPACK better captures peak snow and the ablation period than Snobal. In 2016 Snobal better estimates SWE, as SNOWPACK overestimates during accumulation and for peak SWE. SNOWPACK tends to be more consistent in its prediction capacity, although it tends to overestimate, whereas Snobal tends to underestimate total SWE. Overall, SNOWPACK tends to perform better than Snobal, although there are individual years in which Snobal edges out SNOWPACK.

Figure 7Incoming shortwave radiation for the Marmot Creek Research Basin for 1 February 2011, 17:00 local time. The shadowing algorithm of Dozier and Frew (1990) (DF90) has been implemented on the unstructured mesh. Uniform dark blue indicates shadowed areas.

4.2 Adaptation of raster-based algorithm

Shortwave irradiance corrected for slope and aspect, with horizon (cast) shadows via an adaptation of the Dozier and Frew (1990) shadowing algorithm for unstructured meshes for the Marmot Creek Research Basin, is shown in Fig. 7. The simulation is for 1 February 2011 at 17:00 local time. High irradiance is shown in red, and shadows are shown in dark blue; these areas are receiving only diffuse radiation. The region north of Fisera Ridge is shown in detail in Fig. 8. This figure shows an orthorectified terrestrial photo of a shadow passing over Mt. Collembola from Fisera Ridge. The location of the shadowed region for 1 February 2011 at 17:00 local time is shown for the DF90 algorithm described herein (green), the observed shadow (red), the ArcGIS Solar Analyst shadow (black), and the Marsh 2012 algorithm (blue); the white region is the region not covered by the photograph. The DF90 implementation agrees quite well with observed shadow locations, and a sensitivity test (not shown) shows improved agreement with increasingly small triangles. The performance of the other two shadowing algorithms is detailed in Marsh et al. (2012). In brief, the high-resolution Solar Analyst performed the best but at the cost of multiple hours of runtime. The Marsh 2012 algorithm overpredicted the shadow on the right-hand side of the domain, whereas DF90 underpredicted the shadow delineation. In both cases, the TIN algorithms had a few incorrect shadow classifications on the upper slope as a result of the reduced spatial resolution, causing a small false positive. The triangular-shaped bumps along the shadow line are from the unstructured triangular mesh elements.

Table 4Root mean square error (RMSE, mm) and mean bias (MB, mm) errors averaged over all years.

Download Print Version | Download XLSX

Figure 8This shows an orthorectified terrestrial photo of a shadow passing over Mt. Collembola from Fisera Ridge – details are found in Marsh et al. (2012). The location of the shadowed region for 1 February 2011, 17:00 local time, is shown for the DF90 algorithm described herein (green), the observed shadow (red), the ArcGIS implementation for a 1 m × 1 m lidar raster (black), and for the Marsh et al. (2012) algorithm (blue).

Figure 9Leave-one-out analysis for Vista View (top row), Upper Clearing (middle), and Fisera Ridge (bottom). The values have been binned into 100 hex bins and coloured using the log of the normalized per-bin count. Grey values are bins that have a normalized count of less than 0.01.


4.3 Leave-one-out validation

The leave-one-out validation is shown in Fig. 9 for Vista View (top row), Upper Clearing (middle row), and Fisera Ridge (bottom row). The dashed line is the 1 : 1 line, and the solid black line is a linear regression line of best fit. The r2 value for this fit is shown in the bottom right corner. Due to significant over-plotting of the data points, the values have been binned into 100 hex bins and coloured using the log of the normalized per-bin count. Hex bins divide the xy plane into six-sided bins and counts values in these bins. The hexes avoid the visual artefacts that can occur with square bins. Grey values are bins that have a normalized count of less than 0.01. Because of the significant number of low and zero values in the shortwave and precipitation time series, this resulted in the per-bin colouring being difficult to read. Values of ISWR < 50 W m−2 and p<1 mm were removed for the colouring. Please note that these data were not removed for the linear fit, r2, MBE, or RMSE metrics.

Temperature was well predicted at all sites with r2 values of 0.99, 0.97, and 0.92 for Vista View, Upper Clearing, and Fisera Ridge, respectively. Both mid-elevation sites were better predicted than the high-elevation (Fisera Ridge) site. The majority of the data lie close to the 1 : 1 line. Upper Clearing had a warm bias (MB = 1.11 C), whereas Fisera Ridge had a cold bias (MB =0.37 C). Less spread was observed in the summer months (not shown), matching the results of Cullen and Marshall (2011). Relative humidity was the most poorly predicted variable. Vista View was the most accurately predicted (r2=0.9) with a slight (1.09 %) positive bias. Upper Clearing had more spread with a distinct negative bias (−6.2 %) and decreased r2 (0.76). Fisera Ridge was the most poorly predicted (r2=0.55, MB = 6.12 %). A separate analysis that grouped the data into winter and summer periods (not shown) showed improved results and less spread during the summer months, especially for Fisera Ridge; this summer period had r2=0.7, MB = 7.84 %, and RMSE = 15.32 %. Due to the proximity to vegetation, summer evapotranspiration may result in less temporal variability, dampening the responses. The interpolation methods assume a free atmosphere and thus do not capture these canopy interactions. During the winter months, the observed RH is predominately dominated by synoptic-scale forcing (Cullen and Marshall, 2011) and may be influenced by the sublimation of intercepted snow in the canopy (Pomeroy et al., 2012), which is not captured by this interpolation. The Fisera Ridge data have had substantial infilling for the RH variable (Fang et al., 2019), and the poor fit of CHM to these infilled data may be a result of the infilled data using a higher-elevation exposed ridge that may not be representative of Fisera Ridge. Shortwave irradiance is generally well captured, although Fisera Ridge has a larger negative bias (−15.79 W m−2) than the other two sites. Precipitation at Vista View and Upper Clearing was well predicted, and Fisera Ridge is again the least well predicted.

Figure 10Speed-up for a  100 000 triangle mesh using 1, 2, 4, 6, 8, 16, and 32 cores.


4.4 Parallel scaling

Shown in Fig. 10 are the scaling results for 1, 2, 4, 6, 8, 16, and 32 cores. Good scaling is observed with a 1.97× speed-up with 2 cores, 7.23× speed-up with 8 cores, a 12.3× speed-up with 16 cores, and a 20.5× speed-up with 32 cores. A sublinear scaling is expected due to the mixing of domain and data parallel modules. As most compute nodes are approximately 32 cores, this shows good per-node scaling and thus demonstrates motivation for moving towards a distributed memory model, such as a message passing interface (MPI).

5 Outlook

As described above, this paper outlines the goals of the CHM framework and the motivation for its development. The process representations described herein and currently available in CHM are primarily surface processes. A key component for future development is the addition of the full hydrological cycle to CHM. The use of an irregular geometry for the surface discretization somewhat complicates surface and subsurface computations of lateral mass and energy fluxes. Although an open research topic is how to best incorporate these fluxes in CHM, there are examples of how triangular mesh elements have been used in other models. This section outlines some of the techniques used in these models and how they might be incorporated in CHM.

There are benefits to the improved representation of the irregular geometry for various flow and routing algorithms. If line segments (i.e., triangle edges) are used to represent river and stream channels (along with a sub-grid in-channel flow parameterization), the nonuniform spacing and the unstructured nature of the triangle vertices can more readily represent meandering features than structured meshes (Tucker, 2001). Simple overland flow-routing methods on unstructured meshes can be easily derived to be analogous to the commonly used grid-based flow methods (Tucker, 2001), e.g., D-8 (O'Callaghan and Mark, 1984) and D-inf (Tarboton, 1997). When simulating the landscape evolution of river channels using unstructured and structured meshes, the directional constraints of a maximum of eight directions of a structured mesh have been identified as causing nonnatural channel evolution (Braun and Sambridge, 1997). More complex PDE formulations, such as the shallow water equations, can be discretized on the unstructured mesh. However, the numerical formulations require the use of more sophisticated numerical methods such as the finite-element or finite-volume approaches. The use of these numerical techniques on triangular meshes is common when estimating shallow water flows (Hagen et al., 2002) and has been done with success in hydrological models (Kumar et al., 2009; Qu and Duffy, 2007).

The vertical coupling of the surface flow with subsurface flows requires subsurface extension of the surface mesh. An approach is to use vertically extruded triangles, forming 3-D prisms of some given height that can be stacked vertically. These prisms can then be used to discretize variably saturated flow calculations, such as the Richards equation, and can optionally simulate lateral subsurface flows for a full 3-D model (Hopp et al., 2016; Kumar et al., 2009; Qu and Duffy, 2007). The vertical stacking of prisms is currently being used for aboveground discretization in an in-development blowing snow model for use in CHM (Marsh et al., 2020). Alternatively, vertical 1-D models of variably saturated flow can be used with various lateral assumptions (Hopp et al., 2016). Thus, there are various possibilities for inclusion in CHM, ranging from 1-D models with no lateral or surface coupling to fully coupled 3-D surface–subsurface flows.

One avenue of future research should examine multi-mesh approaches. Depending on the surface and subsurface geology, refinement of the subsurface mesh based on surface characteristics may be inappropriate. A possible approach would result in the generation of one or more subsurface meshes, refined as appropriate to represent the heterogeneity in subsurface properties. This mesh would then be coupled to the surface mesh through various approaches, such as a gradation in triangles linking the subsurface and surface meshes, a nested mesh approach, or via interpolation between the meshes. All approaches involve various technical challenges as well as accuracy trade-offs that have not been intensively explored. Full investigation of the merit of multi-mesh approaches must be done; however, it likely presents an elegant solution to incorporating heterogeneity where appropriate across multiple processes and scales.

6 Conclusion

Simulations of hydrological phenomena are increasingly important for the management and prediction of the hydrological cycle under anthropogenic climate change impacts, and cold regions are some of the most sensitive regions to these impacts. Spatially distributed models are generally thought to produce improved predictions in cold regions when spatially explicit prognostic variables are required; however, substantial challenges including initial conditions, boundary conditions, parameterizations, and computational costs all conspire to limit their applicability. Despite this, hyper-resolution models are increasingly being applied for water management and design decisions. There is a significant opportunity for next-generation models to address challenges in existing models such as seamless prediction at various spatial and temporal scales, the utilization of hyper-resolution data obtained by new remote sensing platforms, quantifying structural uncertainty in distributed models, and the utilization of modern high-performance computing infrastructure.

In this paper, a new modelling framework, the Canadian Hydrological Model (CHM), was presented as a first step towards these goals. Key features of CHM include the ability to do the following: capture spatial heterogeneity in an efficient manner; include multiple process representations; change, remove, and decouple hydrological process algorithms; work at both a point and spatially distributed scale; scale to multiple spatial extents and scales; and utilize a variety of forcing fields (boundary and initial conditions). The efficient representation of spatial heterogeneity is due to the use of unstructured variable-resolution triangular meshes. These can represent key landscape heterogeneities such as vegetation and topography with 50 % to 95 % fewer computational elements versus a fixed-resolution mesh.

To demonstrate and test cold-region operations, two snowpack models, Snobal and SNOWPACK, were compared at a point scale in a mountain clearing. Both models performed well and demonstrated skill in simulating SWE. Although the irregular geometry of a triangular mesh can complicate the application of raster-derived methods, there are various mechanisms in CHM to facilitate the adaptation. A new unstructured mesh implementation of the well-known Dozier and Frew (1990) shadowing algorithm was derived to demonstrate the adaptation of raster-based algorithms and the use of these mechanisms in CHM. This method performed well compared to existing high-resolution raster-based algorithms (Solar Analyst) and other unstructured mesh shadowing algorithms. A leave-one-out validation was done for the meteorological processes and these results showed a high degree of accuracy in the spatial interpolation of meteorological forcing in CHM. Air temperature was the most accurately predicted forcing variable (r2=0.9) and relative humidity was the most poor (r2=0.5 to 0.7). Lastly, a parallel computation scaling test demonstrated a good but sublinear scaling with the number of CPUs and demonstrated a need for increased parallelism efficiency via distributed memory models, such as MPI.

In summary, CHM is a first step towards a variable-resolution explicitly distributed model with a focus for application where cold-region processes play a role in hydrology. Although it remains a work in progress and only snow accumulation and surface meteorology processes are currently implemented, CHM will ultimately include the entirety of the hydrological cycle. The inclusion of irregular geometries is not a significantly problematic aspect for computations of lateral mass and energy exchanges, and other models have used these geometries without issue. However, the novel use of multi-mesh approaches to couple various meshes that have been refined for surface and subsurface optimization is likely a way forward for including increased explicit heterogeneity with a lower computational burden.

Code availability

The code for the Canadian Hydrological Model is open source under the GPLv3 license and is available at (last access: 2 January 2020) and archived with DOI (last access: 2 January 2020). The mesh generation software, Mesher, is open source under the GPLv3 license. It is available at (last access: 2 January 2020).

Author contributions

CM was responsible for the initial idea, coding, analysis, and paper preparation. JP was responsible for idea refinement, analysis refinement, field data experiment design, and paper revision. HW was responsible for idea refinement and paper revision.

Competing interests

The authors declare that they have no conflict of interest.


The authors would like to acknowledge Marmot Creek logistical support from the Nakiska Ski Resort; fieldwork support from Michael Solohub, Xing Fang, May Guan, Angus Duncan, Greg Galloway, and many others is gratefully noted.

Financial support

This research has been conducted with the financial support of the Canada Excellence and Canada Research Chairs programmes, the Natural Sciences and Engineering Research Council of Canada’s Discovery Grants, Alexander Graham Bell Scholarships, Changing Cold Regions Network, Alberta 90 Innovates, and the Global Water Futures programme supported by the Canada First Research Excellence Fund.

Review statement

This paper was edited by Bethanna Jackson and reviewed by Ruzica Dadic and one anonymous referee.


Ahrens, J., Geveci, B., and Law, C.: ParaView: An End-User Tool for Large Data Visualization, in Visualization handbook, Elsevier, 2005. 

Avanzi, F., Michele, C. D., Morin, S., Carmagnola, C. M., and Lejeune, Y.: Model complexity and data requirements in snow hydrology : seeking a balance in practical applications, Hydrol. Proc., 30, 2106–2118,, 2016. 

Bahremand, A.: HESS Opinions: Advocating process modeling and de-emphasizing parameter estimation, Hydrol. Earth Syst. Sci., 20, 1433–1445,, 2016. 

Bartelt, P. and Lehning, M.: A physical SNOWPACK model for the Swiss avalanche warning Part I: numerical model, Cold Reg. Sci. Technol., 35, 123–145,, 2002. 

Bavay, M. and Egger, T.: MeteoIO 2.4.2: a preprocessing library for meteorological data, Geosci. Model Dev., 7, 3135–3151,, 2014. 

Bentley, J. L.: Multidimensional binary search trees used for associative searching, Commun. ACM, 18, 509–517,, 1975. 

Bernhardt, M. and Schulz, K.: SnowSlide: A simple routine for calculating gravitational snow transport, Geophys. Res. Lett., 37, 1–6,, 2010. 

Beven, K.: Changing ideas in hydrology – The case of physically-based models, J. Hydrol., 105, 157–172,, 1989. 

Beven, K.: Prophecy, reality and uncertainty in distributed hydrological modelling, Adv. Water Resour., 16, 41–51,, 1993. 

Beven, K.: A manifesto for the equifinality thesis, Water Resour. Res., 320, 18–36,, 2006. 

Beven, K. and Westerberg, I.: On red herrings and real herrings: Disinformation and information in hydrological inference, Hydrol. Proc., 25, 1676–1680,, 2011. 

Binley, A., Elgy, J., and Beven, K.: A physically based model of heterogeneous hillslopes: 1. Runoff production, Water Resour. Res., 25, 1219–1226,, 1989. 

Blöschl, G. and Sivapalan, M.: Scale issues in hydrological modelling: A review, Hydrol. Proc., 9, 251–290,, 1995. 

Bowling, L. C., Pomeroy, J. W., and Lettenmaier, D. P.: Parameterization of Blowing-Snow Sublimation in a Macroscale Hydrology Model, J. Hydrometeorol., 5, 745–762,<0745:pobsia>;2, 2004. 

Braun, J. and Sambridge, M.: Modelling landscape evolution on geological time scales: a new method based on irregular spatial discretization, Basin Res., 9, 27–52,, 1997. 

Brigode, P., Oudin, L., and Perrin, C.: Hydrological model parameter instability: A source of additional uncertainty in estimating the hydrological impacts of climate change?, Water Resour. Res., 476, 410–425,, 2013. 

Bring, A., Fedorova, I., Dibike, Y., Hinzman, L., Mård, J., Mernild, S. H., Prowse, T., Semenova, O., Stuefer, S. L., and Woo, M.: Arctic terrestrial hydrology: A synthesis of processes, regional effects, and research challenges, J. Geophys. Res.-Biogeo., 121, 621–649,, 2016. 

Burridge, D. M. and Gadd, A. J.: The Meteorological Office Operational 10-level Numerical Weather Prediction Model, Scientific paper – Meteorological Office, (34), 39 pp., 1975. 

Bühler, Y., Adams, M. S., Bösch, R., and Stoffel, A.: Mapping snow depth in alpine terrain with unmanned aerial systems (UASs): potential and limitations, The Cryosphere, 10, 1075–1088,, 2016. 

Carey, S. K. and Woo, M.-k.: Snowmelt Hydrology of Two Subarctic Slopes, Southern Yukon, CanadaPaper Presented at the 11th Northern Res. Basins Symposium/Workshop (Prudhoe Bay to Fairbanks, Alaska, USA, 18–22 August 1997), Hydrol. Res., 29, 331–346,, 1998. 

Chang, K.-T.: Introduction to Geographic Information Systems, McGraw-Hill, New York, New York, 2008. 

Cherkauer, K. a, Bowling, L. C., and Lettenmaier, D. P.: Variable infiltration capacity cold land process model updates, Global Planet. Change, 38, 151–159,, 2003. 

Clark, M. P., Slater, A. G., Rupp, D. E., Woods, R. A., Vrugt, J. A., Gupta, H. V., Wagener, T., and Hay, L. E.: Framework for Understanding Structural Errors (FUSE): A modular framework to diagnose differences between hydrological models, Water Resour. Res., 44, 1–14,, 2008. 

Clark, M. P., Kavetski, D., and Fenicia, F.: Pursuing the method of multiple working hypotheses for hydrological modeling, Water Resour. Res., 47, W09301,, 2011. 

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 process-based hydrologic modeling: 1. Modeling concept, Water Resour. Res., 51, 2498–2514,, 2015. 

Clark, M. P., Bierkens, M. F. P., Samaniego, L., Woods, R. A., Uijlenhoet, R., Bennett, K. E., Pauwels, V. R. N., Cai, X., Wood, A. W., and Peters-Lidard, C. D.: The evolution of process-based hydrologic models: historical challenges and the collective quest for physical realism, Hydrol. Earth Syst. Sci., 21, 3427–3440,, 2017. 

Corripio, J. G.: Snow surface albedo estimation using terrestrial photography, Int. J. Remote Sens., 25, 5705–5729,, 2004. 

Cullen, R. M. and Marshall, S. J.: Mesoscale Temperature Patterns in the Rocky Mountains and Foothills Region of Southern Alberta, Atmos.-Ocean, 49, 189–205,, 2011. 

Das, T., Bárdossy, A., Zehe, E., and He, Y.: Comparison of conceptual model performance using different representations of spatial variability, Water Resour. Res., 356, 106–118,, 2008. 

Davies, T. D., Brimblecombe, P., Tranter, M., Tsiouris, S., Vincent, C. E., Abrahams, P., and Blackwood, I. L.: Seasonal Snowcovers: Physics, Chemistry, Hydrology, D. Reidel Publishing Company, 1987. 

DeBeer, C. M. and Pomeroy, J. W.: Modelling snow melt and snowcover depletion in a small alpine cirque, Canadian Rocky Mountains, Hydrol. Proc., 23, 2584–2599,, 2009. 

DeBeer, C. M., Wheater, H. S., Quinton, W. L., Carey, S. K., Stewart, R. E., MacKay, M. D., and Marsh, P.: The Changing Cold Regions Network: Observation, diagnosis and prediction of environmental change in the Saskatchewan and Mackenzie River Basins, Canada, Science China Earth Sciences, 58, 46–60,, 2015. 

Dodson, R. and Marks, D.: Daily air temperature interpolated at high spatial resolution over a large mountainous region, Climate Res., 8, 1–20,, 1997. 

Dornes, P., Pomeroy, J. W., Pietroniro, A., and Verseghy, D. L.: Effects of Spatial Aggregation of Initial Conditions and Forcing Data on Modeling Snowmelt Using a Land Surface Scheme, J. Hydrometeorol., 9, 789–803,, 2008a. 

Dornes, P. F., Pomeroy, J. W., Pietroniro, A., Carey, S. K., and Quinton, W. L.: Influence of landscape aggregation in modelling snow-cover ablation and snowmelt runoff in a sub-arctic mountainous environment, Hydrol. Sci. J., 53, 725–740,, 2008b. 

Dozier, J. and Frew, J.: Rapid calculation of terrain parameters for radiation modeling from digital elevation data, IEEE T. Geosci. Remote, 28, 963–969,, 1990. 

Duarte, C. M., Lenton, T. M., Wadhams, P., and Wassmann, P.: Abrupt climate change in the Arctic, Nat. Clim. Change, 2, 60–62,, 2012. 

Ellis, C. R., Pomeroy, J. W., Brown, T., and MacDonald, J.: Simulation of snow accumulation and melt in needleleaf forest environments, Hydrol. Earth Syst. Sci., 14, 925–940,, 2010. 

Ellis, C. R. and Pomeroy, J. W.: Estimating sub-canopy shortwave irradiance to melting snow on forested slopes, Hydrol. Proc., 21, 2581–2593,, 2007. 

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,, 2014. 

Essery, R.: A factorial snowpack model (FSM 1.0), Geosci. Model Dev., 8, 3867–3876,, 2015. 

Essery, R., Li, L. and Pomeroy, J.: A distributed model of blowing snow over complex terrain, Hydrol. Proc., 13, 2423–2438,<2423::AID-HYP853>3.0.CO;2-U, 1999. 

Essery, R., Rutter, N., Pomeroy, J., Baxter, R., Gustafsson, D., Barr, A., Bartlett, P., Elder, E., and Stahli, M.: SNOWMIP2: An evaluation of forest snow process simulations, B. Am. Meteorol. Soc., 90, 1120–1135,, 2009. 

Essery, R., Morin, S., Lejeune, Y., and Ménard, C. B.: A comparison of 1701 snow models using observations from an alpine site, Adv. Water Resour., 55, 131–148,, 2013. 

Etchevers, P., Martin, E., Brown, R., Fierz, C., Lejeune, Y., Bazile, E., Boone, A., Dai, Y.-J., Essery, R., Fernandez, A., Gusev, Y., Jordan, R., Koren, V., Kowalczyk, E., Nasonova, N. O., Pyles, R. D., Schlosser, A., Shmakin, A. B., Smirnova, T. G., Strasser, U., Verseghy, D., Yamazaki, T., and Yang, Z.-L.: Validation of the energy budget of an alpine snowpack simulated by several snow models (SnowMIP project), Ann. Glaciol., 38, 150–158,, 2004. 

Fang, X., Pomeroy, J. W., Ellis, C. R., MacDonald, M. K., DeBeer, C. M., and Brown, T.: Multi-variable evaluation of hydrological model predictions for a headwater basin in the Canadian Rocky Mountains, Hydrol. Earth Syst. Sci., 17, 1635–1659,, 2013. 

Fang, X., Pomeroy, J. W., DeBeer, C. M., Harder, P., and Siemens, E.: Hydrometeorological data from Marmot Creek Research Basin, Canadian Rockies, Earth Syst. Sci. Data, 11, 455–471,, 2019. 

Fatichi, S., Vivoni, E. R., Ogden, F. L., Ivanov, V. Y., Mirus, B., Gochis, D., Downer, C. W., Camporese, M., Davison, J. H., Ebel, B., Jones, N., Kim, J., Mascaro, G., Niswonger, R., Restrepo, P., Rigon, R., Shen, C., Sulis, M., and Tarboton, D.: An overview of current applications, challenges, and future trends in distributed process-based models in hydrology, J. Hydrol., 537, 45–60,, 2016. 

Fenicia, F., Kavetski, D., and Savenije, H. H. G.: Elements of a flexible approach for conceptual hydrological modeling: 1. Motivation and theoretical development, Water Resour. Res., 47, 1–13,, 2011. 

Fiddes, J. and Gruber, S.: TopoSCALE v.1.0: downscaling gridded climate data in complex terrain, Geosci. Model Dev., 7, 387–405,, 2014. 

Freeze, R. A.: Streamflow generation, Rev. Geophys., 12, 627–647,, 1974. 

Freeze, R. A. and Harlan, R. L.: Blueprint for a physically-based, digitally-simulated hydrologic response model, J. Hydrol., 9, 237–258,, 1969. 

Fu, P. and Rich, P. M.: Design and implementation of the Solar Analyst: an ArcView extension for modeling solar radiation at landscape scales, in Proceedings of the 19th Annual ESRI User Conference, San Diego, USA, 1–33, 1999. 

Gelfan, A. N., Pomeroy, J. W., and Kuchment, L. S.: Modeling Forest Cover Influences on Snow Accumulation, Sublimation, and Melt, J. Hydrometeorol., 5, 785–803,<0785:mfcios>;2, 2004. 

Golding, D. L.: Research results from Marmot Creek experimental watershed, Alberta, Canada, in IASH Unesco – Symposium on the results of research on representative and experimental basins, Wellington, New Zealand, 397–404, 1970. 

Gray, D. M. and Male, D. H.: Handbook of Snow: Principles, Processes, Management, and Use, Pergamon Press, 1981. 

Gray, D. M., Toth, B., Zhao, L., Pomeroy, J. W., and Granger, R. J.: Estimating areal snowmelt infiltration into frozen soils, Hydrol. Proc., 15, 3095–3111,, 2001. 

Hagen, S. C., Horstmann, O., and Bennett, R. J.: An Unstructured Mesh Generation Algorithm for Shallow Water Modeling, Int. J. Comput. Fluid. D, 16, 83–91,, 2002. 

Harder, P. and Pomeroy, J.: Estimating precipitation phase using a psychrometric energy balance method, Hydrol. Proc., 27, 1901–1914,, 2013. 

Harder, P., Schirmer, M., Pomeroy, J., and Helgason, W.: Accuracy of snow depth estimation in mountain and prairie environments by an unmanned aerial vehicle, The Cryosphere, 10, 2559–2571,, 2016. 

Harder, P., Pomeroy, J. W., and Helgason, W. D.: A simple model for local-scale sensible and latent heat advection contributions to snowmelt, Hydrol. Earth Syst. Sci., 23, 1–17,, 2019. 

Hedstrom, N. R. and Pomeroy, J. W.: Measurements and modelling of snow interception in the boreal forest, Hydrol. Proc., 12, 1611–1625,<1611::aid-hyp684>;2-4, 1998. 

Hopkinson, C., Crasto, N., Marsh, P., Forbes, D., and Lesack, L.: Investigating the spatial distribution of water levels in the Mackenzie Delta using airborne LiDAR, Hydrol. Proc., 25, 2995–3011,, 2011. 

Hopp, L., Fatichi, S., and Ivanov, V. Y.: Simulating water flow in variably saturated soils: a comparison of a 3-D model with approximation-based formulations, Hydrol. Res., 47, 274–290,, 2016. 

Horne, F. E. and Kavvas, M. L.: Physics of the spatially averaged snowmelt process, J. Hydrol., 191, 179–207,, 1997. 

Hrachowitz, M. and Clark, M. P.: HESS Opinions: The complementary merits of competing modelling philosophies in hydrology, Hydrol. Earth Syst. Sci., 21, 3953–3973,, 2017. 

Hubbard, S. S., Gangodagamage, C., Dafflon, B., Wainwright, H., Peterson, J., Gusmeroli, A., Ulrich, C., Wu, Y., Wilson, C., Rowland, J., Tweedie, C., and Wullschleger, S. D.: Quantifying and relating land-surface and subsurface variability in permafrost environments using LiDAR and surface geophysical datasets, Hydrogeol. J., 21, 149–169,, 2013. 

Iqbal, M.: Prediction of hourly diffuse solar radiation from measured hourly global radiation on a horizontal surface, Solar Energy, 24, 491–503,, 1980. 

Ivanov, V., Vivoni, E., Bras, R., and Entekhabi, D.: Preserving high-resolution surface and rainfall data in operational-scale basin hydrology: a fully-distributed physically-based approach, Water Resour. Res., 298, 80–111,, 2004. 

Jones, E., Oliphant, T. and Peterson, P.: SciPy: Open Source Scientific Tools for Python, availablet at: (last access: 8 November 2019), 2018. 

Jordan, R.: A one-dimensional temperature model for a snow cover, Technical documentation for SNTHERM.89, 1991. 

Klemeš, V.: Conceptualization and scale in hydrology, J. Hydrol., 65, 1–23,, 1983. 

Kuchment, L. and Gelfan, A.: Physicomathematical model of snow accumulation and melt in a forest, Russ. Meteorol. Hydrol., 57–65, 2004. 

Kumar, M., Duffy, C. J., and Salvage, K. M.: A Second-Order Accurate, Finite Volume–Based, Integrated Hydrologic Modeling (FIHM) Framework for Simulation of Surface and Subsurface Flow, Vadose Z. J., 8, 873–890,, 2009. 

Kumar, M., Marks, D., Dozier, J., Reba, M., and Winstral, A.: Evaluation of distributed hydrologic impacts of temperature-index and energy-based snow models, Adv. Water Resour., 56, 77–89,, 2013. 

Kunkel, K. E.: Simple procedures for extrapolation of humidity variables in the mountainous western United States, J. Climate, 2, 656–669, 1989. 

Lafaysse, M., Cluzet, B., Dumont, M., Lejeune, Y., Vionnet, V., and Morin, S.: A multiphysical ensemble system of numerical snow modelling, The Cryosphere, 11, 1173–1198,, 2017. 

Latifovic, R. and Pouliot, D.: Analysis of climate change impacts on lake ice phenology in Canada using the historical satellite data record, Remote Sens. Environ., 106, 492–507,, 2007. 

Leavesley, G. H., Markstrom, S. L., Restrepo, P. J., and Viger, R. J.: A modular approach to addressing model design, scale, and parameter estimation issues in distributed hydrological modelling, Hydrol. Proc., 16, 173–187,, 2002. 

Lehning, M., Bartelt, P., Brown, B., Fierz, C., and Satyawali, P.: A physical SNOWPACK model for the Swiss avalanche warning Part II. Snow microstructure, Cold Reg. Sci. Technol., 35, 147–167,, 2002. 

Lehning, M., Völksch, I., Gustafsson, D., Nguyen, T. A., Stähli, M., and Zappa, M.: ALPINE3-D: a detailed model of mountain surface processes and its application to snow hydrology, Hydrol. Proc., 20, 2111–2128,, 2006. 

Lehning, M., Löwe, H., Ryser, M., and Raderschall, N.: Inhomogeneous precipitation distribution and snow transport in steep terrain, Water Resour. Res., 44, W07404,, 2008. 

Leroux, N. R. and Pomeroy, J. W.: Modelling capillary hysteresis effects on preferential flow through melting and cold layered snowpacks, Adv. Water Resour., 107, 250–264,, 2017. 

Li, L. and Pomeroy, J. W.: Estimates of Threshold Wind Speeds for Snow Transport Using Meteorological Data, J. Appl. Meteorol., 36, 205–213,<0205:eotwsf>;2, 1997. 

Liston, G. E. and Elder, K.: A Meteorological Distribution System for High-Resolution Terrestrial Modeling (MicroMet), J. Hydrometeorol., 7, 217–234,, 2006. 

Lundberg, A., Ala-Aho, P., Eklo, O., Klöve, B., Kværner, J., and Stumpp, C.: Snow and frost: implications for spatiotemporal infiltration patterns – a review, Hydrol. Proc., 30, 1230–1250,, 2016. 

MacDonald, M. K., Pomeroy, J. W., and Pietroniro, A.: Parameterizing redistribution and sublimation of blowing snow for hydrological models: tests in a mountainous subarctic catchment, Hydrol. Proc., 23, 2570–2583,, 2009. 

Marks, D., Dozier, J., and Davis, R. E.: Climate and energy exchange at the snow surface in the Alpine Region of the Sierra Nevada: 1. Meteorological measurements and monitoring, Water Resour. Res., 28, 3029–3042,, 1992. 

Marks, D., Kimball, J., Tingey, D., and Link, T.: The sensitivity of snowmelt processes to climate conditions and forest cover during rain-on-snow: a case study of the 1996 Pacific Northwest flood, Hydrol. Proc., 12, 1569–1587, 1998. 

Marks, D., Domingo, J., Susong, D., Link, T., and Garen, D.: A spatially distributed energy balance snowmelt model for application in mountain basins, Hydrol. Proc., 13, 1935–1959,<1935::aid-hyp868>;2-c, 1999. 

Marks, D., Winstral, A., Reba, M., Pomeroy, J., and Kumar, M.: An evaluation of methods for determining during-storm precipitation phase and the rain/snow transition elevation at the surface in a mountain basin, Adv. Water Resour., 55, 98–110,, 2013. 

Marsh, C. B., Pomeroy, J. W., and Spiteri, R. J.: Implications of mountain shading on calculating energy for snowmelt using unstructured triangular meshes, Hydrol. Proc., 26, 1767–1778,, 2012. 

Marsh, C. B., Spiteri, R. J., Pomeroy, J. W., and Wheater, H. S.: Multi-objective unstructured triangular mesh generation for use in hydrological and land surface models, Comput. Geosci., 119, 49–67,, 2018. 

Marsh, C. B., Pomeroy, J. W., Spiteri, R. J., and Wheater, H. S.: A finite volume blowing snow model for use with variable resolution meshes, Water Resour. Res., 55, e24400,, 2020. 

Marty, C., Philipona, R., Fröhlich, C., and Ohmura, A.: Altitude dependence of surface radiation fluxes and cloud forcing in the alps: results from the alpine surface radiation budget network, Theor. Appl. Climatol., 72, 137–155,, 2002. 

Mason, P. and Sykes, R.: Flow over an isolated hill of moderate slope, Q. J. Roy. Meteor. Soc., 105, 383–395,, 1979. 

Maxwell, R. M. and Kollet, S. J.: Interdependence of groundwater dynamics and land-energy feedbacks under climate change, Nat. Geosci., 1, 665–669,, 2008. 

McCauley, C. A., White, D. M., Lilly, M. R., and Nyman, D. M.: A comparison of hydraulic conductivities, permeabilities and infiltration rates in frozen and unfrozen soils, Cold Reg. Sci. Technol., 34, 117–125,, 2002. 

Mendoza, P. A., Clark, M. P., Barlage, M., Rajagopalan, B., Samaniego, L., Abramowitz, G., and Gupta, H.: Are we unnecessarily constraining the agility of complex process-based models?, Water Resour. Res., 51, 716–728,, 2015. 

Michlmayr, G., Lehning, M., Koboltschnig, G., Holzmann, H., Zappa, M., Mott, R., and Schöner, W.: Application of the Alpine 3-D model for glacier mass balance and glacier runoff studies at Goldbergkees, Austria, Hydrol. Proc., 22, 3941–3949,, 2008. 

Milly, P. C. D., Betancourt, J., Falkenmark, M., Hirsch, R. M., Kundzewicz, Z., Lettenmaier, D. P., and Stouffer, R. J.: Stationarity Is Dead: Whither Water Management?, Science, 319, 573–574, 2008. 

Mohanty, B. P.: Soil Hydraulic Property Estimation Using Remote Sensing: A Review, Vadose Z. J., 12,, 2013. 

Mosier, T. M., Hill, D. F., and Sharp, K. V.: How much cryosphere model complexity is just right? Exploration using the conceptual cryosphere hydrology framework, The Cryosphere, 10, 2147–2171,, 2016. 

Mote, P. W., Hamlet, A. F., Clark, M. P., and Lettenmaier, D. P.: Declining Mountain Snowpack in Western North America, B. Am. Meteorol. Soc., 86, 39–49,, 2005. 

Mott, R., Schirmer, M., Bavay, M., Grünewald, T., and Lehning, M.: Understanding snow-transport processes shaping the mountain snow-cover, The Cryosphere, 4, 545–559,, 2010. 

Mott, R., Gromke, C., Grünewald, T., and Lehning, M.: Relative importance of advective heat transport and boundary layer decoupling in the melt dynamics of a patchy snow cover, Adv. Water Resour., 55, 88–97,, 2013. 

Munro, D. S. and Young, G. J.: An operational net shortwave radiation model for glacier basins, Water Resour. Res., 18, 220–230,, 1982. 

Musselman, K. N., Pomeroy, J. W., and Link, T. E.: Variability in shortwave irradiance caused by forest gaps: Measurements, modelling, and implications for snow energetics, Agr. Forest Meteorol., 207, 69–82,, 2015. 

Musselman, K. N., Clark, M. P., Liu, C., Ikeda, K., and Rasmussen, R.: Slower snowmelt in a warmer world, Nat. Clim. Change, 7, 214–219,, 2017. 

Nazemi, A., Wheater, H. S., Chun, K. P., and Elshorbagy, A.: A stochastic reconstruction framework for analysis of water resource system vulnerability to climate-induced changes in river flow regime, Water Resour. Res., 49, 291–305,, 2013. 

O'Callaghan, J. F. and Mark, D. M.: The extraction of drainage networks from digital elevation data, Comput. Vision Graph., 28, 323–344,, 1984. 

Oliphant, T. E.: A guide to NumPy, Trelgol Publishing, USA, 2006. 

Olyphant, G. A.: Longwave Radiation in Mountainous Areas and Its Influence on the Energy Balance of Alpine Snowfields, Water Resour. Res., 22, 62–66,, 1986. 

Or, D., Lehmann, P., and Assouline, S.: Natural length scales define the range of applicability of the Richards equation for capillary flows, Water Resour. Res., 51, 7130–7144,, 2015. 

Painter, S. L., Coon, E. T., Atchley, A. L., Berndt, M., Garimella, R., Moulton, J. D., Svyatskiy, D., and Wilson, C. J.: Integrated surface/subsurface permafrost thermal hydrology: Model formulation and proof-of-concept simulations, Water Resour. Res., 52, 6062–6077,, 2016. 

Paniconi, C. and Putti, M.: Physically based modeling in catchment hydrology at 50: Survey and outlook, Water Resour. Res., 51, 7090–7129,, 2015. 

Perrin, C., Michel, C., and Andréassian, V.: Does a large number of parameters enhance model performance? Comparative assessment of common catchment model structures on 429 catchments, J. Hydrol., 242, 275–301,, 2001. 

Pietroniro, A., Fortin, V., Kouwen, N., Neal, C., Turcotte, R., Davison, B., Verseghy, D., Soulis, E. D., Caldwell, R., Evora, N., and Pellerin, P.: Development of the MESH modelling system for hydrological ensemble forecasting of the Laurentian Great Lakes at the regional scale, Hydrol. Earth Syst. Sci., 11, 1279–1294,, 2007. 

Pomeroy, J. and Bernhardt, M.: Project Report for the 2017 GEWEX GHP Meeting, Kathmandu, Nepal, 2017. 

Pomeroy, J. W. and Gray, D. M.: Saltation of snow, Water Resour. Res., 26, 1583–1594,, 1990. 

Pomeroy, J. W. and Li, L.: Prairie and arctic areal snow cover mass balance using a blowing snow model, J. Geophys. Res.-Atmos., 105, 26619–26634,, 2000. 

Pomeroy, J. W. and Male, D. H.: Steady-state suspension of snow, J. Hydrol., 136, 275–301,, 1992. 

Pomeroy, J. W., Gray, D. M., and Landine, P. G.: The Prairie Blowing Snow Model: characteristics, validation, operation, J. Hydrol., 144, 165–192,, 1993. 

Pomeroy, J. W., Gray, D. M., Shook, K. R., Toth, B., Essery, R., Pietroniro, A., and Hedstrom, N.: An evaluation of snow accumulation and ablation processes for land surface modelling, Hydrol. Proc., 12, 2339–2367,<2339::AID-HYP800>3.0.CO;2-L, 1998a. 

Pomeroy, J. W., Parviainen, J., Hedstrom, N., and Gray, D. M.: Coupled modelling of forest snow interception and sublimation, Hydrol. Proc., 12, 2317–2337,<2317::aid-hyp799>;2-x, 1998b. 

Pomeroy, J. W., Toth, B., Granger, R. J., Hedstrom, N. R., and Essery, R. L. H.: Variation in surface energetics during snowmelt in a subarctic mountain catchment, J. Hydrometeorol., 4, 702–719,<0702:viseds>;2, 2003. 

Pomeroy, J. W., Granger, R. J., Hedstrom, N. R., Gray, D. M., Elliott, J., Pietroniro, A., and Janowicz, J. R.: The Process Hydrology Approach to Improving Prediction of Ungauged Basins in Canada, in Prediction in Ungauged Basins: Approaches for Canada's Cold Regions, Environment Canada, 67–100, 2004. 

Pomeroy, J. W., Gray, D. M., Brown, T., Hedstrom, N. R., Quinton, W. L., Granger, R. J., and Carey, S. K.: The cold regions hydrological model: a platform for basing process representation and model structure on physical evidence, Hydrol. Proc., 21, 2650–2667,, 2007. 

Pomeroy, J. W., Marks, D., Link, T., Ellis, C., Hardy, J., Rowlands, A., and Granger, R.: The impact of coniferous forest temperature on incoming longwave radiation to melting snow, Hydrol. Proc., 23, 2513–2525,, 2009. 

Pomeroy, J., Fang, X., and Ellis, C.: Sensitivity of snowmelt hydrology in Marmot Creek, Alberta, to forest cover disturbance, Hydrol. Proc., 26, 1891–1904,, 2012. 

Pomeroy, J. W., Fang, X., Shook, K., and Whitfield, P. H.: Predicting in ungauged basins using physical principles obtained using the deductive, inductive, and abdyctive reasoning approach, in: Putting prediction in ungauged basins into practice, Canadian Water Resources Association, 41–62, 2013. 

Qu, Y. and Duffy, C. J.: A semidiscrete finite volume formulation for multiprocess watershed simulation, Water Resour. Res., 43, 1–18,, 2007. 

Raderschall, N., Lehning, M., and Schär, C.: Fine-scale modeling of the boundary layer wind field over steep topography, Water Resour. Res., 44, 1–18,, 2008. 

Raleigh, M. S., Livneh, B., Lapo, K., and Lundquist, J. D.: How does availability of meteorological forcing data impact physically-based snowpack simulations?, J. Hydrometeorol., 99–120,, 2015. 

Rasouli, K., Pomeroy, J. W., and Marks, D. G.: Snowpack sensitivity to perturbed climate in a cool mid-latitude mountain catchment, Hydrol. Proc., 29, 3925–3940,, 2015. 

Razavi, S. and Gupta, H. V.: A new framework for comprehensive, robust, and efficient global sensitivity analysis: 1. Theory, Water Resour. Res., 52, 423–439,, 2016. 

Rew, R. and Davis, G.: NetCDF: An Interface for Scientific Data Access, IEEE Comput. Graph., 10, 76–82,, 1990. 

Rothwell, R., Hillman, G., and Pomeroy, J. W.: Marmot Creek Experimental Watershed Study, Forestry Chron., 92, 32–36,, 2016. 

Rouse, W. R., Oswald, C. J., Binyamin, J., Spence, C., Schertzer, W. M., Blanken, P. D., Bussières, N., and Duguay, C. R.: The Role of Northern Lakes in a Regional Energy Balance, J. Hydrometeorol., 6, 291–305,, 2005. 

Samaniego, L., Kumar, R., Thober, S., Rakovec, O., Zink, M., Wanders, N., Eisner, S., Müller Schmied, H., Sutanudjaja, E. H., Warrach-Sagi, K., and Attinger, S.: Toward seamless hydrologic predictions across spatial scales, Hydrol. Earth Syst. Sci., 21, 4323–4346,, 2017. 

Savenije, H. H. G.: HESS Opinions “The art of hydrology”, Hydrol. Earth Syst. Sci., 13, 157–161,, 2009. 

Schlögl, S., Lehning, M., Fierz, C., and Mott, R.: Representation of Horizontal Transport Processes in Snowmelt Modeling by Applying a Footprint Approach, Front. Earth Sci., 6, 120,, 2018. 

Schroeder, W., Martin, K. and Lorensen, B.: The Visualization Toolkit, 4th edn., Kitware, 2006. 

Shewchuk, J.: Triangle: engineering a 2D quality mesh generator and Delaunay triangulator, in Applied computational geometry towards geometric engineering, Springer Berlin/Heidelberg, 203–222, 1996. 

Shook, K. and Gray, D. M.: Small-scale Spatial Structure Of Shallow Snowcovers, Hydrol. Proc., 10, 1283–1292,<1283::aid-hyp460>;2-m, 1996. 

Shook, K., Pomeroy, J., and van der Kamp, G.: The transformation of frequency distributions of winter precipitation to spring streamflow probabilities in cold regions; case studies from the Canadian Prairies, J. Hydrol., 521, 395–409,, 2015. 

Sicart, J. E., Pomeroy, J. W., Essery, R. L. H., and Bewley, D.: Incoming longwave radiation to melting snow: observations, sensitivity and estimation in Northern environments, Hydrol. Proc., 20, 3697–3708,, 2006. 

Sivapalan, M.: From engineering hydrology to Earth system science: milestones in the transformation of hydrologic science, Hydrol. Earth Syst. Sci., 22, 1665–1693,, 2018. 

Slater, A. G., Barrett, A. P., Clark, M. P., Lundquist, J. D., and Raleigh, M. S.: Uncertainty in seasonal snow reconstruction: Relative impacts of model forcing and image availability, Adv. Water Resour., 55, 165–177,, 2013. 

Smith, C. D.: The relationship between snowfall catch efficiency and wind speed for the Geonor T-200B precipitation gauge utilizing various wind shield configurations, in: Proceedings of the 77th Annual Western Snow Conference, Canmore, AB, 115–121, 2009. 

Spence, C. and Mengistu, S.: Deployment of an unmanned aerial system to assist in mapping an intermittent stream, Hydrol. Proc., 30, 493–500,, 2016. 

Tangelder, H. and Fabri, A.: dD Spatial Searching, in: CGAL user and reference manual, version 4.10, available at:, last access: Januray 2018. 

Tarboton, D. G.: A new method for the determination of flow directions and upslope areas in grid digital elevation models, Water Resour. Res., 33, 309–319,, 1997. 

Thornton, P. E., Running, S. W., and White, M. A.: Generating surfaces of daily meteorological variables over large regions of complex terrain, Water Resour. Res., 190, 214–251,, 1997. 

Todini, E.: Rainfall-runoff modeling – Past, present and future, Water Resour. Res., 100, 341–352,, 1988. 

Tucker, G.: An object-oriented framework for distributed hydrologic and geomorphic modeling using triangulated irregular networks, Comput. Geosci., 27, 959–973,, 2001. 

Vaze, J., Post, D. A., Chiew, F. H. S., Perraud, J. M., Viney, N. R., and Teng, J.: Climate non-stationarity – Validity of calibrated rainfall–runoff models for use in climate change studies, Water Resour. Res., 394, 447–457,, 2010.  

Verseghy, D. L.: Class – A Canadian land surface scheme for GCMS, I. Soil model, Int. J. Climatol., 11, 111–133,, 1991. 

Verseghy, D. L., McFarlane, N. A., and Lazare, M.: Class – A Canadian land surface scheme for GCMS, II. Vegetation model and coupled runs, Int. J. Climatol., 13, 347–370,, 1993. 

Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le Moigne, P., Martin, E., and Willemet, J.-M.: The detailed snowpack scheme Crocus and its implementation in SURFEX v7.2, Geosci. Model Dev., 5, 773–791,, 2012. 

Viviroli, D., Dürr, H. H., Messerli, B., Meybeck, M., and Weingartner, R.: Mountains of the world, water towers for humanity: Typology, mapping, and global significance, Water Resour. Res., 43, W07447,, 2007. 

Vrugt, J. A., ter Braak, C. J. F., Clark, M. P., Hyman, J. M., and Robinson, B. A.: Treatment of input uncertainty in hydrologic modeling: Doing hydrology backward with Markov chain Monte Carlo simulation, Water Resour. Res., 44, W00B09,, 2008. 

Wagener, T. and Montanari, A.: Convergence of approaches toward reduci ng uncertainty in predictions in ungauged basins, Water Resour. Res., 47, W06301,, 2011. 

Walcek, C. J.: Cloud cover and its relationship to relative humidity during a springtime midlatitude cyclone, Mon. Weather Rev., 122, 1021–1035, 1994. 

Walvoord, M. A. and Kurylyk, B. L.: Hydrologic Impacts of Thawing Permafrost – A Review, Vadose Z. J., 15,, 2016. 

Wayand, N. E., Marsh, C. B., Shea, J. M., and Pomeroy, J. W.: Globally scalable alpine snow metrics, Remote Sens. Environ., 213, 61–72,, 2018. 

Wheater, H. S.: Water Security – science and management challenges, Proc. IAHS, 366, 23–30,, 2015. 

Winstral, A., Elder, K., and Davis, R. E.: Spatial snow modeling of wind-redistributed snow using terrain-based parameters, J. Hydrometeorol., 3, 524–538,<0524:ssmowr>;2, 2002. 

Wood, E. F., Roundy, J. K., Troy, T. J., Beek, L. P. H. van, Bierkens, M. F. P., Blyth, E., Roo, A. de, Döll, P., Ek, M., Famiglietti, J., Gochis, D., Giesen, N. van de, Houser, P., Jaffé, P. R., Kollet, S., Lehner, B., Lettenmaier, D. P., Peters-Lidard, C., Sivapalan, M., Sheffield, J., Wade, A., and Whitehead, P.: Hyperresolution global land surface modeling: Meeting a grand challenge for monitoring Earth's terrestrial water, Water Resour. Res., 47, W05301,, 2011. 

Zhao, L. and Gray, D. M.: Estimating snowmelt infiltration into frozen soils, Hydrol. Proc., 13, 1827–1842,<1827::aid-hyp896>;2-d, 1999. 

Short summary
The Canadian Hydrological Model (CHM) is a next-generation distributed model. Although designed to be applied generally, it has a focus for application where cold-region processes, such as snowpacks, play a role in hydrology. A key feature is that it uses a multi-scale surface representation, increasing efficiency. It also enables algorithm comparisons in a flexible structure. Model philosophy, design, and several cold-region-specific examples are described.