Land models are often applied at a regularly spaced grid. Land models are typically set up at a range of spatial configurations, ranging from grid sizes of 0.02 to 2^∘ (approximately 2 to 200 km) and applied at sub-daily temporal resolutions for the simulation of energy fluxes. A priori specification of the grid size of the land models is often derived from forcing resolutions, modeling objectives, available geospatial data, and computational resources and is usually based on modeling convenience. Figure 1e–h illustrate the typical land model configuration – here the modeler selects a cell size, and then the soil, vegetation, and forcing files are all aggregated or disaggregated to the target cell size. Original data resolution and spatial distribution of soil, land cover, and forcing data are smeared while upscaled to the resolution of interest. Any change in the modeling resolution will require upscaling or downscaling of the geophysical data set once again.

In this study, we configure the land models using non-regular shapes. Figure 1a–d present an example of non-regular shapes created through spatial intersections of the land covers and soil type shapes. These vector-based configurations of the geospatial data are then forced at the original meteorological forcing resolution or its upscaled or downscaled values, resulting in computational units. Therefore, each computational unit has unique geospatial data such as soil, vegetation, slope, and aspect and is forced with a unique forcing. In this configuration, changing meteorological forcing resolution does not affect the decisions needed to upscale the geospatial data such as soil type and land cover to the grid resolution.

Figure 1Panels (a–d) indicate vector-based configuration of a land model. (a) Meteorological forcing at its original resolution or upscaled and downscaled resolutions, (b) land covers, (c) soil types with their spatial extent, and (d) vector-based configuration with 28 computational units each with unique forcing, soil type, and land cover type. The bottom row indicates typical grid-based configuration of a land model; (e) a priori resolution should be decided; (f) meteorological forcing should be upscaled or downscaled to the grid resolution; (g) land cover percentage should be calculated for each modeling grid, or a dominate land cover should be selected to represent that grid; and finally (h) soil characteristics for each modeling grid should be identified.

Download

The benefits of vector-based configuration of land models can be summarized as follows:

1. There is no need for a priori assumption on modeling grid size. In traditional land model implementation, the modeler selects a grid resolution (which is often a regular latitude/longitude grid). The soil parameters and forcing data from any resolution must be aggregated, disaggregated, resampled, or interpolated for every grid size. The land cover data are often only considered as a percentage for every grid, and spatial location of the land cover is lost. However, in the vector-based setup, these decisions are only based on the input and forcing data that are chosen to be used in the modeling practice, and no upscaling or downscaling to grid size is needed. Furthermore, the size of computational units can vary across the modeling domain depending on the variability of the meteorological forcing and geospatial heterogeneity. For example, the spatial density of computational units can be higher in mountainous areas, where temperature and precipitation gradients are larger while avoiding an unnecessarily high number of computational units in areas with a lower gradient in meteorological forcing.

2. There is a reasonable relation between available meteorological forcing and geospatial data resolution and the number of computational units. Computational units that are the result of available geophysical data sets forced with the original forcing data logically represent the maximum number of computational units that can be hydrologically unique. A higher number of computational units than the proposed setup will arguably provide an unnecessary computational burden due to identical forcing data and geospatial information.

3. Direct simplification of geospatial data is possible. The vector-based implementation facilitates easier aggregation of computational units. It is easier to aggregate similar soil types or similar forested areas into unified shapes with a basic GIS function (dissolving for example) than this would be if all data had to be upscaled or downscaled into a different grid size.

4. Direct specification of physical parameters is possible, and unrealistic combinations of land cover, soil, and other geophysical information are avoided. As each computational unit has a specific type of land cover, soil type, and other physical characteristics, it is straightforward to specify parameter values based on lookup tables (i.e., no averaging; upscaling is needed). This is favorable because the modeler does not need to make decisions about methods used for upscaling of geophysical data at the grid level. Also, this might avoid the unrealistic combination of parameter sets that might be considered by the model at a grid scale, such as equiprobable combination of land cover on soil type which may not exist in reality, which will be increasing the fidelity of the model representation of the processes (we will elaborate this further in the context of the VIC model in the Discussion). We emphasize that the ease of parameter allocation for vector-based implementation of land models does not address the challenge of finding the right parameter sets for each computational unit.

5. It has the ability to compare and constrain the parameter values for computational units and their simulations. The impact of land cover, soil type and elevation zone can be evaluated separately. For example, the vector-based implementation makes it easier to test if forested areas generate less surface runoff than grasslands. This might be more challenging at the gird-based configuration, in which there are a combination of different land cover types at grid scale. Similarly, the vector-based implementation may simplify regularization efforts across large geographical domains. These relative constraints can be utilized to translate often patchy expert knowledge into a sophisticated land model.

6. It is possible to incorporate additional data. If needed, additional data, such as slope and aspect, for example, can be incorporated in building the computational units, accounting for changes in shortwave radiation or lapse rates for temperature. The changes can be implemented outside of the model in the forcing files. Computational units can be built also based on variation of the leaf area index (LAI), giving an additional layer of information in addition to the land cover type. The additional information can be easily ingested into the model without extra effort in contrast to changing the model parameter files at the grid scale.

7. Comparison of model simulations and in situ point-scale observation and visualization are easier. The vector-based implementation of land models makes it easier to compare the point measurement to model simulation as the model simulations preserve the extent of geospatial features.

8. The selection of models is modular and controlled. The vector-based implementation identifies the characteristics and spatial boundary of geospatial domains. A model might not be suitable for processes of some of the geospatial domains. Alternatively, processes of a computational unit that is beyond the capacity of one model can be replaced with an alternative model. For example, computational units that are glaciered can be replaced with more suitable models, while the spatial configuration and forcings remain identical. Consequently, the effect of features such as glaciers can be better studied as more expert models can be applied to glaciers, while the rest of the computational units can be simulated with a model that includes general processes.

In this study, we proposed a vector-based configuration for land models and applied this setup to the VIC model. We used a vector-based routing scheme, mizuRoute, which was forced using output from the land model (one-way coupling). Unlike the grid-based approach, there is no upscaling of land cover percentage or soil characteristics to a new grid size. This enables us to separate the effects of changes in forcing resolution from changes in the spatial configurations. As mentioned earlier in Sect. 2, the vector-based configuration of land models may help to avoid unrealistic configuration of soil type, land cover, or elevation zones that may happen in traditional grid-based implementation and hence increase the model fidelity. As an example, VIC configuration at grid scale assumes equal distribution of land cover over different elevation zones. Figure 7b illustrates how the traditional VIC configuration at grid scale wrongly considers forested land cover above the treeline. This issue is avoided in the vector-based configuration as the setup will only include two computational units of forested area below the treeline and bare soil above the treeline (Fig. 7a). The vector-based setup provides more flexibility in comparing the model simulations across computational units (as an example, refer to Fig. 5) and also comparing model simulations with point measurements, such as snow water equivalent. Moreover, the vector-based routing results in complete decoupling of the land model computational units' spatial extent from routing sub-basins. For the grid-based configuration of land models, it is often the case that in land, model grids and routing grids are identical, which result in further decisions on upscaling of the routing direction to the land model grid scale.

Figure 7(a) The realistic configuration of a natural system with land cover consisting of 50 % bare soil and 50 % forest within a grid located in two different elevation zones above and below the treeline which is preserved with vector-based configurations and (b) the traditional VIC configurations for the given system at the grid for the two elevation zones and two land cover types which results in an unrealistic combination of forested land cover above the treeline and bare soil below the treeline.

Download

Our results illustrate that various vector-based spatial configuration of the VIC model generates similar large-scale simulations of streamflow when the setups are calibrated by maximizing the Nash–Sutcliffe score at the basin outlet. Similarly, we have shown that often behavioral parameter sets yield similar E_NS values and can be significantly uncertain (Fig. 3a) or have significant differences for their internal behavior, which may be very well masked by aggregation of the result at the grid scale or basin scale (Figs. 3b and 5). Generally, parameter, state and flux uncertainties are not often evaluated or reported on in land models (Demaria et al., 2007) or are ignored by tying parameters, linking specific hydraulic conductivity to the slope of the water retention curve, for example, so that the possible combinations of parameters are reduced. Moreover, the behavior of the K_slow parameter can be revealing of the VIC model structural deficiencies, which are not often explored for land models. The recession coefficient obtained from recession analysis on the observed hydrograph is approximately 0.01 d⁻¹, while the calibrated K_slow has much higher values of around 0.90 d⁻¹. This can be due to a dampened response from the two top soil layers and a lack of macropore water movement to the baseflow component. Similarly, and due to a lack of macropore water movement in the VIC model, and land models in general, it is impossible to infer the K_slow based on recession analysis on the observed hydrograph (for further reading on this and also recession analysis refer to Gharari et al., 2019). This finding can be generalized to the five-parameter VIC baseflow, highlighting the need to properly evaluate the often not observable but calibrated baseflow parameters for the VIC model and if it is possible to identify five parameters based on the recession limps of a hydrograph.

Land models are often applied at large spatial scales. The results clearly show that the deviation of streamflow is much lower in river segments with a larger upstream area (Figs. 4 and 6a). It is often the case that the model parameters and associated processes are inferred through calibration on the streamflow at the basin outlet or over a large contributing area. We argue that this may not be a valid strategy for process understanding at the smaller scale (read computational units), given the large uncertainty exhibited by the parameters. Therefore, hyper-resolution modeling efforts (Wood et al., 2011) may suffer from poor process representation and parameter identification at the scale of interest (Beven et al., 2015). What is needed instead of efficiency metrics that aggregate model behavior across both space (e.g., at the outlet of the larger catchment) and time (e.g., expressing the mismatch between observations and simulations across the entire observation period as a single number) is diagnostic evaluation of the model's process fidelity at the scale at which simulations are generated in the case of available observations (e.g., Gupta et al., 2008; Clark et al., 2016).

One might argue that the spatial discretization is important for the realism of model fluxes and states. Moving to a significantly high number of computational units may result in computational units that are similar in their forcing and geospatial fabric (such as soil and land cover types). Based on the result of this study for snow water equivalent (Fig. 5), we can argue that the snow patterns are fairly similar for the configurations that have elevation zones and finer resolution of forcing (Cases 3 and 4 and forcing resolution less than 0.125^∘). ($m\left(\stackrel{\mathrm{‾}}{x}\mathrm{|}\stackrel{\mathrm{‾}}{\mathit{\theta }}\right)\sim \stackrel{\mathrm{‾}}{m\left(x\mathrm{|}\mathit{\theta }\right)}$, in which m is the model, x and θ are model forcing and the model parameter set, and $\stackrel{\mathrm{‾}}{x}$ and $\stackrel{\mathrm{‾}}{\mathit{\theta }}$ are upscaled forcing and parameter values at coarser spatial representation.)

The analysis of the accuracy–efficiency trade-off presented in this study, Fig. 6, can be used in model analysis such as sensitivity and uncertainty. One can assume that a configuration with fewer computational units can be a surrogate for a model with more computational units, under the condition that both models are known to behave similarly for a given parameter set. The calibration can be done on the model configuration with fewer computational units, and the parameters can be transferred directly to the model with more computational units or can be used as an initial point for optimization algorithm to speed up the calibration process. Similarly, the sensitivity analyses can be done primarily on the model with fewer computational units.

In this study, and following the concept hydrological similarity, we assume the parameters of computational units are identical for computational units with similar soil and land cover. The degree of validity of hydrological similarity concepts is debatable. For example, at the catchment scale, Oudin et al. (2010) have shown that the overlap between catchments with similar physiographic attributes and catchments with similar model performance for a given parameter set is only 60 %. Thus, physiographic similarity (in our case expressed through GRUs) does not necessarily imply similarity of hydrologic behavior, even though this is the critical assumption underlying GRUs. The VIC parameters can be linked to many more characteristics such as slope, height above nearest drainage (HAND; Renno et al., 2008; Gharari et al., 2011), or topographic wetness index (Beven and Kirkby, 1979), as has been done by Mizukami et al. (2016) and Chaney et al. (2018). Techniques such as multiscale parameter regionalization (MPR; Samaniego et al., 2010) can be used to scale parameter values for different model configurations. However, the functions that are used to link computational units and physical attributes to model parameters remain mostly based on inference, (i.e., calibration), and the reproducibility of those relationships are not very well explored. However, application of these techniques, such as in this case that has significant parameter and process uncertainty and significant accuracy–efficiency trade-off, should be put through rigorous tests (Merz et al., 2020; Liu et al., 2016).

A key outstanding challenge is for models to provide the right results for the right reasons (Kirchner, 2006). Thoughtful strategies to formulate parameter and process constraints based on expert knowledge can reduce the plausible range of behavioral parameter sets. In this study, we imposed a simple parameter constraint that the root zone moisture storage of forested area should be larger than the non-forested area (Table 1). Additional process constraints, if available, can be increasingly difficult to satisfy. More rigorous parameter estimation methods that satisfy the fidelity constraints based on expert knowledge are required (e.g., Gharari et al., 2014).

In this study, the vector-based routing configuration does not include lakes and reservoirs. This is often a neglected element of land modeling efforts and has only attracted limited attention compared to the its impact on terrestrial water cycle (Haddeland et al., 2006; Yassin et al., 2019). The presence of lakes and reservoirs and their interconnections reduces the, already limited, ability of inference of land model parameters based on calibration on the observed streamflow (Dang et al., 2020).

Although not primarily the result of this study, however, the nivo-glacial regime of the Bow River basins is mostly dominated by snowmelt that contributes mostly to streamflow through baseflow (slow component of the hydrograph). The high Nash–Sutcliffe efficiency, E_NS, is partly due to the fact that it is rather easy for the land model to capture the yearly cycle of the streamflow only with snow processes (see, e.g., Knoben et al., 2020, demonstrating this for the Kling–Gupta efficiency), while rapid subsurface water movement, such as macropore movement, is largely missing in the land models (Gharari et al., 2019). Therefore, more caution is needed for the calibration of land model parameters for flood forecasting (Vionnet et al., 2020) for the Bow region and all the nivo-glacial river systems in western Canada, McKenzie, Yukon, and Colombia River basins.

All the data used in this study are available publicly (refer to references).

SG, AP, and MPC developed the concept. SG did the model setup and simulations. SG and MC wrote the paper. WJMK, NM, and JSW contributed to the writing and editing of the paper, scientific comments, interpretation of results, and preparation of figures.

The authors declare that they have no conflict of interest.

This paper was edited by Niko Wanders and reviewed by two anonymous referees.