Abstract Space-based, global-extent digital elevation models (DEMs) are key inputs to many Earth sciences applications. However, many of these applications require the use of a ‘bare-Earth’ DEM versus a digital surface model (DSM), the latter of which may include systematic positive biases due to tree canopies in forested areas. Critical topographic features may be obscured by these biases. Vegetation-free datasets have been created by using statistical relationships and machine learning to train on local-scale datasets (e.g., lidar) to de-bias the global-extent datasets. Recent advances in satellite platforms coupled with increased availability of computational resources and lidar reference products has allowed for a new generation of vegetation- and urban-canopy removals. One of these is the Forest And Buildings removed Copernicus DEM (FABDEM), based on the most recent and most accurate global DSM Copernicus-30. Among the more challenging landscapes to quantify surface elevations are densely forested mountain catchments, where even airborne lidar applications struggle to capture surface returns. The increasing affordability and availability of UAV-based lidar platforms have resulted in new capacity to fly modest spatial extents with unrivalled point densities. These data allow an unprecedented ability to validate global sub-canopy DEMs against representative UAV-based lidar data. In this work, the FABDEM is validated against up-scaled lidar data in a steep and forested mountain catchment considering elevation, slope, and Terrain Position Index (TPI) metrics. Comparisons of FABDEM with SRTM, MERIT, and the Copernicus-30 dataset are made. It was found that the FABDEM had a 24% reduction in elevation RMSE and a 135% reduction in bias compared to the Copernicus-30 dataset. Overall, the FABDEM provides a clear improvement over existing deforested DEM products in complex mountain topography such as the MERIT DEM. This study supports the use of FABDEM in forested mountain catchments as the current best-in-class data product.

Estimates of near-surface wind speed and direction are key meteorological components for predicting many surface hydrometeorological processes that influence critical aspects of hydrological and biological systems. However, observations of near-surface wind are typically spatially sparse. The use of these sparse wind fields to force distributed models, such as hydrological models, is greatly complicated in complex terrain, such as mountain headwaters basins. In these regions, wind flows are heavily impacted by overlapping influences of terrain at different scales. This can have a great impact on calculations of evapotranspiration, snowmelt, and blowing snow transport and sublimation. The use of high-resolution atmospheric models allows for numerical weather prediction (NWP) model outputs to be dynamically downscaled. However, the computation burden for large spatial extents and long periods of time often precludes their use. Here, a wind-library approach is presented to aid in downscaling NWP outputs and terrain-correcting spatially interpolated observations. This approach preserves important spatial characteristics of the flow field at a fraction of the computational costs of even the simplest high-resolution atmospheric models. This approach improves on previous implementations by: scaling to large spatial extents O(1M km2); approximating lee-side effects; and fully automating the creation of the wind library. Overall, this approach was shown to have a third quartile RMSE of 1.8 and a third quartile RMSE of 58.2° versus a standalone diagnostic windflow model. The wind velocity estimates versus observations were better than existing empirical terrain-based estimates and computational savings were approximately 100-fold versus the diagnostic model.

• Snow, glaciers, wetlands, frozen ground and permafrost needed in hydrological models. • Water quality export by coupling biochemical transformations to cold regions processes. • Hydrological sensitivity to land use depends on cold regions processes. • Strong cold regions hydrological sensitivity to climate warming. Cold regions involve hydrological processes that are not often addressed appropriately in hydrological models. The Cold Regions Hydrological Modelling platform (CRHM) was initially developed in 1998 to assemble and explore the hydrological understanding developed from a series of research basins spanning Canada and international cold regions. Hydrological processes and basin response in cold regions are simulated in a flexible, modular, object-oriented, multiphysics platform. The CRHM platform allows for multiple representations of forcing data interpolation and extrapolation, hydrological model spatial and physical process structures, and parameter values. It is well suited for model falsification, algorithm intercomparison and benchmarking, and has been deployed for basin hydrology diagnosis, prediction, land use change and water quality analysis, climate impact analysis and flood forecasting around the world. This paper describes CRHM’s capabilities, and the insights derived by applying the model in concert with process hydrology research and using the combined information and understanding from research basins to predict hydrological variables, diagnose hydrological change and determine the appropriateness of model structure and parameterisations.

Despite the proliferation of computer-based research on hydrology and water resources, such research is typically poorly reproducible. Published studies have low reproducibility due to incomplete availability of data and computer code, and a lack of documentation of workflow processes. This leads to a lack of transparency and efficiency because existing code can neither be quality controlled nor reused. Given the commonalities between existing process-based hydrologic models in terms of their required input data and preprocessing steps, open sharing of code can lead to large efficiency gains for the modeling community. Here, we present a model configuration workflow that provides full reproducibility of the resulting model instantiations in a way that separates the model-agnostic preprocessing of specific data sets from the model-specific requirements that models impose on their input files. We use this workflow to create large-domain (global and continental) and local configurations of the Structure for Unifying Multiple Modeling Alternatives (SUMMA) hydrologic model connected to the mizuRoute routing model. These examples show how a relatively complex model setup over a large domain can be organized in a reproducible and structured way that has the potential to accelerate advances in hydrologic modeling for the community as a whole. We provide a tentative blueprint of how community modeling initiatives can be built on top of workflows such as this. We term our workflow the “Community Workflows to Advance Reproducibility in Hydrologic Modeling” (CWARHM; pronounced “swarm”).

Distributed hydrological models predict the spatial variability in processes that govern observed mass and energy fluxes. A challenge associated with the use of these models is the computational burden associated with representing the Earth's (sub)surface via millions of computational elements. This burden is exacerbated as more complex process representations are included because their parameterizations involve computationally intensive mathematical functions. Lookup tables (LUTs) approximate a mathematical function by interpolating precomputed values of the function. Highly accurate approximations are possible for substantially reduced computational costs. In this work, a general methodology using the C++ LUT library FunC is applied to identify and replace computationally intensive mathematical function evaluations in the Canadian Hydrological Model (CHM). The use of LUTs introduces a pointwise relative error below 10 − 8 and provides a reduction in run time of almost 20%. This work shows how LUTs can be implemented with relatively little pain and yield significant computational savings for distributed hydrological models. • The Canadian Hydrological Model (CHM) is profiled and expensive mathematical functions identified. • FunC was used to replace the expensive mathematical functions in CHM with lookup tables. • The run-time performance of CHM was improved by approximately 20% on two realistic simulations. • A general methodology for using FunC to replace expensive mathematical functions with lookup tables is given.

Abstract. The interaction of mountain terrain with meteorological processes causes substantial temporal and spatial variability in snow accumulation and ablation. Processes impacted by complex terrain include large-scale orographic enhancement of snowfall, small-scale processes such as gravitational and wind-induced transport of snow, and variability in the radiative balance such as through terrain shadowing. In this study, a multi-scale modelling approach is proposed to simulate the temporal and spatial evolution of high-mountain snowpacks. The multi-scale approach combines atmospheric data from a numerical weather prediction system at the kilometre scale with process-based downscaling techniques to drive the Canadian Hydrological Model (CHM) at spatial resolutions allowing for explicit snow redistribution modelling. CHM permits a variable spatial resolution by using the efficient terrain representation by unstructured triangular meshes. The model simulates processes such as radiation shadowing and irradiance to slopes, blowing-snow transport (saltation and suspension) and sublimation, avalanching, forest canopy interception and sublimation, and snowpack melt. Short-term, kilometre-scale atmospheric forecasts from Environment and Climate Change Canada's Global Environmental Multiscale Model through its High Resolution Deterministic Prediction System (HRDPS) drive CHM and are downscaled to the unstructured mesh scale. In particular, a new wind-downscaling strategy uses pre-computed wind fields from a mass-conserving wind model at 50 m resolution to perturb the mesoscale HRDPS wind and to account for the influence of topographic features on wind direction and speed. HRDPS-CHM was applied to simulate snow conditions down to 50 m resolution during winter 2017/2018 in a domain around the Kananaskis Valley (∼1000 km2) in the Canadian Rockies. Simulations were evaluated using high-resolution airborne light detection and ranging (lidar) snow depth data and snow persistence indexes derived from remotely sensed imagery. Results included model falsifications and showed that both wind-induced and gravitational snow redistribution need to be simulated to capture the snowpack variability and the evolution of snow depth and persistence with elevation across the region. Accumulation of windblown snow on leeward slopes and associated snow cover persistence were underestimated in a CHM simulation driven by wind fields that did not capture lee-side flow recirculation and associated wind speed decreases. A terrain-based metric helped to identify these lee-side areas and improved the wind field and the associated snow redistribution. An overestimation of snow redistribution from windward to leeward slopes and subsequent avalanching was still found. The results of this study highlight the need for further improvements of snowdrift-permitting models for large-scale applications, in particular the representation of subgrid topographic effects on snow transport.

Abstract. The interaction of mountain terrain with meteorological processes causes substantial temporal and spatial variability in snow accumulation and ablation. Processes impacted by complex terrain include large-scale orographic enhancement of snowfall, small-scale processes such as gravitational and wind-induced transport of snow, and variability in the radiative balance such as through terrain shadowing. In this study, a multi-scale modeling approach is proposed to simulate the temporal and spatial evolution of high mountain snowpacks using the Canadian Hydrological Model (CHM), a multi-scale, spatially distributed modelling framework. CHM permits a variable spatial resolution by using the efficient terrain representation by unstructured triangular meshes. The model simulates processes such as radiation shadowing and irradiance to slopes, blowing snow redistribution and sublimation, avalanching, forest canopy interception and sublimation and snowpack melt. Short-term, km-scale atmospheric forecasts from Environment and Climate Change Canada's Global Environmental Multiscale Model through its High Resolution Deterministic Prediction System (HRDPS) drive CHM, and were downscaled to the unstructured mesh scale using process-based procedures. In particular, a new wind downscaling strategy combines meso-scale HRDPS outputs and micro-scale pre-computed wind fields to allow for blowing snow calculations. HRDPS-CHM was applied to simulate snow conditions down to 50-m resolution during winter 2017/2018 in a domain around the Kananaskis Valley (~1000 km2) in the Canadian Rockies. Simulations were evaluated using high-resolution airborne Light Detection and Ranging (LiDAR) snow depth data and snow persistence indexes derived from remotely sensed imagery. Results included model falsifications and showed that both blowing snow and gravitational snow redistribution need to be simulated to capture the snowpack variability and the evolution of snow depth and persistence with elevation across the region. Accumulation of wind-blown snow on leeward slopes and associated snow-cover persistence were underestimated in a CHM simulation driven by wind fields that did not capture leeside flow recirculation and associated wind speed decreases. A terrain-based metric helped to identify these lee-side areas and improved the wind field and the associated snow redistribution. An overestimation of snow redistribution from windward to leeward slopes and subsequent avalanching was still found. The results of this study highlight the need for further improvements of snowdrift-permitting models for large-scale applications, in particular the representation of subgrid topographic effects on snow transport.

Blowing snow is ubiquitous in cold, windswept environments. In some regions, blowing snow sublimation losses can ablate a notable fraction of the seasonal snowfall. It is advantageous to predict alpine snow regimes at the spatial scale of snowdrifts (≈1 to 100 m) because of the role of snow redistribution in governing the duration and volume of snowmelt. However, blowing snow processes are often neglected due to computational costs. Here, a three‐dimensional blowing snow model is presented that is spatially discretized using a variable resolution unstructured mesh. This represents the heterogeneity of the surface explicitly yet, for the case study reported, gained a 62% reduction in computational elements versus a fixed‐resolution mesh and resulted in a 44% reduction in total runtime. The model was evaluated for a subarctic mountain basin using transects of measured snow water equivalent (SWE) in a tundra valley. Including blowing snow processes improved the prediction of SWE by capturing inner‐annual snowdrift formation, more than halved the total mean bias error, and increased the coefficient of variation of SWE from 0.04 to 0.31 better matching the observed CV (0.41). The use of a variable resolution mesh did not dramatically degrade the model performance. Comparison with a constant resolution mesh showed a similar CV and RMSE as the variable resolution mesh. The constant resolution mesh had a smaller mean bias error. A sensitivity analysis showed that snowdrift locations and immediate up‐wind sources of blowing snow are the most sensitive areas of the landscape to wind speed variations.

Abstract. 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.

Recent studies of water flow through dry porous media have shown progress in simulating preferential flow propagation. However, current methods applied to snowpacks have neglected the dynamic nature of the capillary pressure, such as conditions for capillary pressure overshoot, resulting in a rather limited representation of the water flow patterns through snowpacks observed in laboratory and field experiments. Indeed, previous snowmelt models using a water entry pressure to simulate preferential flow paths do not work for natural snowpack conditions where snow densities are less than 380 kg m−3. Because preferential flow in snowpacks greatly alters the flow velocity and the timing of delivery of meltwater to the base of a snowpack early in the melt season, a better understanding of this process would aid hydrological predictions. This study presents a 2‐D water flow through snow model that solves the non‐equilibrium Richards equation. This model, coupled with random perturbations of snow properties, can represent realistic preferential flow patterns. Using 1‐D laboratory data, two model parameters were linked to snow properties and model boundary conditions. Parameterizations of these model parameters were evaluated against 2‐D snowpack observations from a laboratory experiment, and the resulting model sensitivity to varying inputs and boundary conditions was calculated. The model advances both the physical understanding of and ability to simulate water flow through snowpacks and can be used in the future to parameterize 1‐D snowmelt models to incorporate flow variations due to preferential flow path formation.

Abstract Unstructured triangular meshes are an efficient and effective landscape representation that are suitable for use in distributed hydrological and land surface models. Their variable spatial resolution provides similar spatial performance to high-resolution structured grids while using only a fraction of the number of elements. Many existing triangulation methods either sacrifice triangle quality to introduce variable resolution or maintain well-formed uniform meshes at the expense of variable triangle resolution. They are also generally constructed to only fulfil topographic constraints. However, distributed hydrological and land surface models require triangles of varying resolution to provide landscape representations that accurately represent the spatial heterogeneity of driving meteorology, physical parameters and process operation in the simulation domain. As such, mesh generators need to constrain the unstructured mesh to not only topography but to other important surface and sub-surface features. This work presents novel multi-objective unstructured mesh generation software that allows mesh generation to be constrained to an arbitrary number of important features while maintaining a variable spatial resolution. Triangle quality is supported as well as a smooth gradation from small to large triangles. Including these additional constraints results in a better representation of spatial heterogeneity than from classic topography-only constraints.

Abstract Horizontal and altitudinal redistribution of snow by wind transport and avalanches can be important controls on small- and large-scale snow accumulation patterns that control meltwater supply in alpine environments. Redistribution processes control the spatial variability of snow accumulation, which not only controls meltwater supply, but also regulates snowmelt timing, duration, and rates, as well as snow-covered area depletion and the variable contributing area for meltwater runoff generation. However, most hydrological models and land surface schemes do not consider snow redistribution processes, and those that do are difficult to verify without spatially distributed snow depth measurements. These are rarely available in both high resolution and covering large scales. As an increased number of hydrological models include snow redistribution processes there is a need for additional snowcover metrics to verify snow redistribution schemes over large areas using readily available data. This study develops novel high-resolution (20 m), snowcover indices from remotely sensed imagery (Landsat-8 and Sentinel-2) to evaluate snow redistribution models over alpine areas without in-situ or airborne snow observations. A snowcover absence (SA) index, calculated from snow-free areas during the winter, identifies areas of wind erosion or avalanche source areas. A snowcover persistence (SP) index, calculated from snow-covered areas during the summer, identifies snow deposition in drifts and avalanche deposits. The snowcover indices captured the relative differences in surface observations of snow presence and absence between exposed and sheltered sites on an intensely instrumented ridge in the Canadian Rockies Hydrological Observatory. Within the Tuolumne River Basin in central California (1100 km2), the SP index captured roughly half of the spatial variability (R2 = 0.49 to 0.56) in peak SWE as estimated from airborne LiDAR-derived snow depths. At the individual mountain ridge scale (~800 m), variability in both ablation and snow redistribution controlled the SP patterns over 7979 ridges. Differences in shortwave irradiance explained 76% of the SP variance across ridges, but could not explain smaller-scale (~100 m) SP peaks that are associated with snowdrifts and avalanche deposits. The snowcover indices can be used to evaluate snow redistribution models of the finer scale impacts of snow redistribution by wind and gravity as long as the larger scale influences of spatially variable solar irradiance effects are also simulated.