Search

In this Reply, we address the concerns of Knoben and Clark (2023, https://doi.org/10.1029/2022WR034294) or KC23 that “the assumptions needed to effectively use difference-based variance estimation methods are not always met by hourly streamflow records.” There should be little doubt that the assumptions of our difference-based estimator will sometimes be violated in hourly streamflow records but the results from de Oliveira and Vrugt (2022, https://doi.org/10.1029/2022wr032263) and confirmed by the findings in our Reply show that such violations are sporadic enough not to affect much the error variance estimates. Snowmelt as pointed out by KC23 (https://doi.org/10.1029/2022WR034294) may not have received sufficient attention in our paper, yet their 365-day record is simply not long enough to demonstrate bias of our discharge error variance estimates (and their dependence on flow level). This would require analysis of a much longer, multi-year, streamflow record of a snowmelt-driven watershed. The snowmelt catchment analyzed in dOV22 (https://doi.org/10.1029/2022wr032263) did not demonstrate bias in the discharge error variance estimates. We also provide additional clarification on the interpretation of the variance estimates obtained with the nonparametric estimator, and discuss the main issues in the test case presented in Knoben and Clark (2023, https://doi.org/10.1029/2022WR034294)).

Variance-based analysis has emerged as method of choice for quantifying the sensitivity of the output, y, of a scalar-valued square-integrable function, f ∈ L2((Formula presented.)), to its d ≥ 1 input variables, x = {x1, …, xd}, with support (Formula presented.). The prototype of this approach, Sobol's method is a generalization of the analysis of variance (ANOVA) to d > 2 independent input variables and decomposes y, as sum of elementary functions of zeroth-, first-, second-, up to dth-order. This independence assumption is mathematically convenient but may not be borne out of the causal or correlational relationships between the x's. This paper is concerned with variance-based sensitivity analysis (SA) for correlated input variables, for example, multivariate dependencies in a posterior parameter distribution. We use high-dimensional model representation (HDMR) of Li et al. (2010, https://doi.org/10.1021/jp9096919), Li and Rabitz (2012, https://doi.org/10.1007/s10910-011-9898-0) and replace Sobol's elementary functions with so-called component functions with unknown expansion coefficients to disentangle the structural, correlative and total contribution of input factors. We contrast the default HDMR methodology with cubic B-splines and sequential coefficient estimation against its successor, HDMRext of Li and Rabitz (2012, https://doi.org/10.1007/s10910-011-9898-0), which uses polynomial component functions with an extended orthonormalized basis. Benchmark experiments confirm that HDMR and HDMRext parse out the structural and correlative contributions of input factors to the model output and infer an optimal experimental design with parameter correlation. Our last study applies HDMRext to probabilistic SA of a watershed model. The multivariate posterior parameter distribution supports model emulation and yields sensitivity indices that pertain to measured discharge data.

In this paper, we introduce a relatively simple data-driven method for the representation of the uncertainty in daily discharge records. The proposed method relies only on hourly discharge data and takes advantage of a nonparametric difference-based estimator in the characterization of random errors in discharge time series. We illustrate with corrupted streamflow data that the nonparametric estimator provides an accurate characterization of the nature (homoscedastic or heteroscedastic) and magnitude of these errors. In addition, we demonstrate the practical usefulness of the estimator using discharge time series of 500+ watersheds of the Catchment Attributes and MEteorology for Large-sample Studies data set. This analysis reveals that the magnitude of errors of aleatory nature in the investigated discharge records is rather small (less than 3% for 80% of the records). We then combine the effect of random errors and measurement frequency into a daily variance estimate, which serves as input to a streamflow generation approach. This procedure produces replicates of the discharge record which portray accurately the assigned streamflow uncertainty, preserve key statistical properties of the discharge record and are hydrologically realistic. The proposed method facilitates Bayesian analysis and supports tasks such as model diagnostics, data assimilation, uncertainty quantification and regionalization.

Nonstationary extreme value analysis (NEVA) can improve the statistical representation of observed flood peak distributions compared to stationary (ST) analysis, but management of flood risk relies on predictions of out-of-sample distributions for which NEVA has not been comprehensively evaluated. In this study, we apply split-sample testing to 1250 annual maximum discharge records in the United States and compare the predictive capabilities of NEVA relative to ST extreme value analysis using a log-Pearson Type III (LPIII) distribution. The parameters of the LPIII distribution in the ST and nonstationary (NS) models are estimated from the first half of each record using Bayesian inference. The second half of each record is reserved to evaluate the predictions under the ST and NS models. The NS model is applied for prediction by (1) extrapolating the trend of the NS model parameters throughout the evaluation period and (2) using the NS model parameter values at the end of the fitting period to predict with an updated ST model (uST). Our analysis shows that the ST predictions are preferred, overall. NS model parameter extrapolation is rarely preferred. However, if fitting period discharges are influenced by physical changes in the watershed, for example from anthropogenic activity, the uST model is strongly preferred relative to ST and NS predictions. The uST model is therefore recommended for evaluation of current flood risk in watersheds that have undergone physical changes. Supporting information includes a MATLAB® program that estimates the (ST/NS/uST) LPIII parameters from annual peak discharge data through Bayesian inference.

What is the “best” model? The answer to this question lies in part in the eyes of the beholder, nevertheless a good model must blend rigorous theory with redeeming qualities such as parsimony and quality of fit. Model selection is used to make inferences, via weighted averaging, from a set of K candidate models, (Formula presented.), and help identify which model is most supported by the observed data, (Formula presented.). Here, we introduce a new and robust estimator of the model evidence, (Formula presented.), which acts as normalizing constant in the denominator of Bayes’ theorem and provides a single quantitative measure of relative support for each hypothesis that integrates model accuracy, uncertainty, and complexity. However, (Formula presented.) is analytically intractable for most practical modeling problems. Our method, coined GAussian Mixture importancE (GAME) sampling, uses bridge sampling of a mixture distribution fitted to samples of the posterior model parameter distribution derived from MCMC simulation. We benchmark the accuracy and reliability of GAME sampling by application to a diverse set of multivariate target distributions (up to 100 dimensions) with known values of (Formula presented.) and to hypothesis testing using numerical modeling of the rainfall-runoff transformation of the Leaf River watershed in Mississippi, USA. These case studies demonstrate that GAME sampling provides robust and unbiased estimates of the evidence at a relatively small computational cost outperforming commonly used estimators. The GAME sampler is implemented in the MATLAB package of DREAM and simplifies considerably scientific inquiry through hypothesis testing and model selection.

Distribution forecasts P over future quantities or events are routinely made in hydrology but usually traded for a (likelihood‐weighted) mean or median prediction to accommodate error measures or scoring functions such as the mean absolute error or mean squared error. Case in point is the so‐called KG efficiency (KGE) of Gupta et al. (2009, https://doi.org/10.1016/j.jhydrol.2009.08.003) and improvements thereof (Lamontagne et al., 2020, https://doi.org/10.1029/2020wr027101), which have rapidly gained popularity among hydrologists as alternative scoring functions to the commonly used Nash and Sutcliffe (1970, https://doi.org/10.1016/0022‐1694(70)90255‐6) efficiency, but are equally exclusive in how they quantify model performance using only single‐valued output of the quantities of interest. This point‐valued mapping necessarily implies a loss of information about model performance. This paper advocates the use of probabilistic watershed model training, evaluation and diagnostics. Distribution evaluation opens a mature literature on scoring rules whose strong statistical underpinning provides, as we will demonstrate, the theory, context and guidelines necessary for the development of robust information‐theoretically principled metrics for watershed signatures. These so‐called hydrograph functionals are scalar‐valued mappings of major behavioral watershed functions embodied in a strictly proper scoring rule. We discuss past developments that led to the current state‐of‐the‐art of distribution evaluation in hydrology and review scoring rules for dichotomous and categorical events, quantiles (intervals) and density forecasts. We are particularly concerned with elicitable functionals and scoring rule propriety, discuss the decomposition of scoring rules into a sharpness, reliability and entropy term and present diagnostically appealing strictly proper divergence scores of hydrograph functionals for flood frequency analysis, flow duration and recession curves. The usefulness and power of distribution‐based model evaluation and diagnostics by means of scoring rules is demonstrated on simple illustrative problems and discharge distributions simulated with watershed models using random sampling and Bayesian model averaging. The presented theory (a) enables a more complete evaluation of distribution forecasts, (b) offers a statistically principled means for watershed model training, evaluation, diagnostics and selection using hydrograph functionals and/or extreme events and (c) provides a universal framework for metric development of watershed signatures, promoting metric standardization and reproducibility. Plain Language Summary: The past decades have witnessed an unbridled growth in goodness‐of‐fit metrics of hydrologic models. These metrics may satisfy the needs of hydrologists but lack conforming theory and principles. This state of affairs (a) elicits improper model training and evaluation, (b) provokes and supports misguided inferences, (c) impedes statistically‐principled uncertainty quantification, metric standardization and development of universal model benchmarks and (d) obfuscates determination of whether the model has finished learning. What is more, most hydrologic model evaluation metrics in use today are rather exclusive in how they quantify model performance using only single‐valued simulated output of the quantities of interest. Predictive distributions derived from (quasi)‐Bayesian methods or ensembles are usually traded for a (likelihood‐weighted) mean or median prediction to accommodate error measures (scoring functions) such as the mean absolute error. This implies a large loss of information. This paper develops a distribution‐based approach to hydrologic model evaluation and diagnostics. Distribution evaluation opens the necessary theory and guidelines for development of robust information‐theoretically principled metrics of watershed signatures. These so‐called hydrograph functionals are scalar‐valued mappings of major behavioral watershed functions embodied in a strictly proper scoring rule. The hydrograph functionals offer a statistically principled means for hydrologic model evaluation, diagnostics and selection. Key Points: Scoring rules of hydrograph functionals provide an information‐theoretically principled means for watershed model training, evaluation, and diagnosticsWe present strictly proper (divergence) scores for flood frequency analysis, flow duration, and recession curvesPropriety and elicitability offer useful working paradigms for metric development of hydrograph functionals [ABSTRACT FROM AUTHOR]

The depth to bedrock controls a myriad of processes by influencing subsurface flow paths, erosion rates, soil moisture, and water uptake by plant roots. As hillslope interiors are very difficult and costly to illuminate and access, the topography of the bedrock surface is largely unknown. This essay is concerned with the prediction of spatial patterns in the depth to bedrock (DTB) using high-resolution topographic data, numerical modeling, and Bayesian analysis. Our DTB model builds on the bottom-up control on fresh-bedrock topography hypothesis of Rempe and Dietrich (2014) and includes a mass movement and bedrock-valley morphology term to extent the usefulness and general applicability of the model. We reconcile the DTB model with field observations using Bayesian analysis with the DREAM algorithm. We investigate explicitly the benefits of using spatially distributed parameter values to account implicitly, and in a relatively simple way, for rock mass heterogeneities that are very difficult, if not impossible, to characterize adequately in the field. We illustrate our method using an artificial data set of bedrock depth observations and then evaluate our DTB model with real-world data collected at the Papagaio river basin in Rio de Janeiro, Brazil. Our results demonstrate that the DTB model predicts accurately the observed bedrock depth data. The posterior mean DTB simulation is shown to be in good agreement with the measured data. The posterior prediction uncertainty of the DTB model can be propagated forward through hydromechanical models to derive probabilistic estimates of factors of safety.

Many watershed models used within the hydrologic research community assume (by default) stationary conditions, that is, the key watershed properties that control water flow are considered to be time invariant. This assumption is rather convenient and pragmatic and opens up the wide arsenal of (multivariate) statistical and nonlinear optimization methods for inference of the (temporally fixed) model parameters. Several contributions to the hydrologic literature have brought into question the continued usefulness of this stationary paradigm for hydrologic modeling. This paper builds on the likelihood-free diagnostics approach of Vrugt and Sadegh () and uses a diverse set of hydrologic summary metrics to test the stationary hypothesis and detect changes in the watersheds response to hydroclimatic forcing. Models with fixed parameter values cannot simulate adequately temporal variations in the summary statistics of the observed catchment data, and consequently, the DREAM(ABC) algorithm cannot find solutions that sufficiently honor the observed metrics. We demonstrate that the presented methodology is able to differentiate successfully between watersheds that are classified as stationary and those that have undergone significant changes in land use, urbanization, and/or hydroclimatic conditions, and thus are deemed nonstationary.

We present a Bayesian inversion method for the joint inference of high-dimensional multi-Gaussian hydraulic conductivity fields and associated geostatistical parameters from indirect hydrological data. We combine Gaussian process generation via circulant embedding to decouple the variogram from grid cell specific values, with dimensionality reduction by interpolation to enable Markov chain Monte Carlo (MCMC) simulation. Using the Matérn variogram model, this formulation allows inferring the conductivity values simultaneously with the field smoothness (also called Matérn shape parameter) and other geostatistical parameters such as the mean, sill, integral scales and anisotropy direction(s) and ratio(s). The proposed dimensionality reduction method systematically honors the underlying variogram and is demonstrated to achieve better performance than the Karhunen-Loève expansion. We illustrate our inversion approach using synthetic (error corrupted) data from a tracer experiment in a fairly heterogeneous 10,000-dimensional 2-D conductivity field. A 40-times reduction of the size of the parameter space did not prevent the posterior simulations to appropriately fit the measurement data and the posterior parameter distributions to include the true geostatistical parameter values. Overall, the posterior field realizations covered a wide range of geostatistical models, questioning the common practice of assuming a fixed variogram prior to inference of the hydraulic conductivity values. Our method is shown to be more efficient than sequential Gibbs sampling (SGS) for the considered case study, particularly when implemented on a distributed computing cluster. It is also found to outperform the method of anchored distributions (MAD) for the same computational budget. Key Points: Joint Bayesian inference of Gaussian conductivity fields and their variograms A dimensionality reduction that systematically honors the underlying variogram Distributed multiprocessor implementation is straightforward

The quest for a more powerful method for model evaluation has inspired Vrugt and Sadegh (2013) to introduce ‘‘likelihood-free’’ inference as vehicle for diagnostic model evaluation. This class of methods is also referred to as Approximate Bayesian Computation (ABC) and relaxes the need for a residual based likelihood function in favor of one or multiple different summary statistics that exhibit superior diagnostic power. Here we propose several methodological improvements over commonly used ABC sampling methods to permit inference of complex system models. Our methodology entitled DREAM(ABC) uses the DiffeRential Evolution Adaptive Metropolis algorithm as its main building block and takes advantage of a continuous fitness function to efficiently explore the behavioral model space. Three case studies demonstrate that DREAM(ABC) is at least an order of magnitude more efficient than commonly used ABC sampling methods for more complex models. DREAM(ABC) is also more amenable to distributed, multi-processor, implementation, a prerequisite to diagnostic inference of CPU-intensive system models.

The ever increasing pace of computational power, along with continued advances in measurement technologies and improvements in process understanding has stimulated the development of increasingly complex hydrologic models that simulate soil moisture flow, groundwater recharge, surface runoff, root water uptake, and river discharge at different spatial and temporal scales. Reconciling these high-order system models with perpetually larger volumes of field data is becoming more and more difficult, particularly because classical likelihood-based fitting methods lack the power to detect and pinpoint deficiencies in the model structure. Gupta et al. (2008) has recently proposed steps (amongst others) toward the development of a more robust and powerful method of model evaluation. Their diagnostic approach uses signature behaviors and patterns observed in the input-output data to illuminate to what degree a representation of the real world has been adequately achieved and how the model should be improved for the purpose of learning and scientific discovery. In this paper, we introduce approximate Bayesian computation (ABC) as a vehicle for diagnostic model evaluation. This statistical methodology relaxes the need for an explicit likelihood function in favor of one or multiple different summary statistics rooted in hydrologic theory that together have a clearer and more compelling diagnostic power than some average measure of the size of the error residuals. Two illustrative case studies are used to demonstrate that ABC is relatively easy to implement, and readily employs signature based indices to analyze and pinpoint which part of the model is malfunctioning and in need of further improvement. ©2013. American Geophysical Union. All Rights Reserved.

Simultaneous scaling of soil water retention and hydraulic conductivity functions provides an effective means to characterize the heterogeneity and spatial variability of soil hydraulic properties in a given study area. The statistical significance of this approach largely depends on the number of soil samples collected. Unfortunately, direct measurement of the soil hydraulic functions is tedious, expensive and time consuming. Here we present a simple and cost-effective hybrid scaling approach that combines the use of ancillary information (e.g., particle-size distribution and soil bulk density) with direct measurements of saturated soil water content and saturated hydraulic conductivity. Our results demonstrate that the presented approach requires far fewer laboratory measurements than conventional scaling methods to adequately capture the spatial variability of soil hydraulic properties. © 2013. American Geophysical Union. All Rights Reserved.

Markov Chain Monte Carlo (MCMC) methods have become increasingly popular for estimating the posterior probability distribution of parameters in hydrologic models. However, MCMC methods require the a priori definition of a proposal or sampling distribution, which determines the explorative capabilities and efficiency of the sampler and therefore the statistical properties of the Markov Chain and its rate of convergence. In this paper we present an MCMC sampler entitled the Shuffled Complex Evolution Metropolis algorithm (SCEM-UA), which is well suited to infer the posterior distribution of hydrologic model parameters. The SCEM-UA algorithm is a modified version of the original SCE-UA global optimization algorithm developed by Duan et al. [1992]. The SCEM-UA algorithm operates by merging the strengths of the Metropolis algorithm, controlled random search, competitive evolution, and complex shuffling in order to continuously update the proposal distribution and evolve the sampler to the posterior target distribution. Three case studies demonstrate that the adaptive capability of the SCEM-UA algorithm significantly reduces the number of model simulations needed to infer the posterior distribution of the parameters when compared with the traditional Metropolis-Hastings samplers.

Practical experience with the calibration of hydrologic models suggests that any single-objective function, no matter how carefully chosen, is often inadequate to properly measure all of the characteristics of the observed data deemed to be important. One strategy to circumvent this problem is to define several optimization criteria (objective functions) that measure different (complementary) aspects of the system behavior and to use multicriteria optimization to identify the set of nondominated, efficient, or Pareto optimal solutions. In this paper, we present an efficient and effective Markov Chain Monte Carlo sampler, entitled the Multiobjective Shuffled Complex Evolution Metropolis (MOSCEM) algorithm, which is capable of solving the multiobjective optimization problem for hydrologic models. MOSCEM is an improvement over the Shuffled Complex Evolution Metropolis (SCEM-UA) global optimization algorithm, using the concept of Pareto dominance (rather than direct single-objective function evaluation) to evolve the initial population of points toward a set of solutions stemming from a stable distribution (Pareto set). The efficacy of the MOSCEM-UA algorithm is compared with the original MOCOM-UA algorithm for three hydrologic modeling case studies of increasing complexity.

We have developed a sequential optimization methodology, entitled the parameter identification method based on the localization of information (PIMLI) that increases information retrieval from the data by inferring the location and type of measurements that are most informative for the model parameters. The PIMLI approach merges the strengths of the generalized sensitivity analysis (GSA) method [Spear and Hornberger, 1980], the Bayesian recursive estimation (BARE) algorithm [Thiemann et al., 2001], and the Metropolis algorithm [Metropolis et al., 1953]. Three case studies with increasing complexity are used to illustrate the usefulness and applicability of the PIMLI methodology. The first two case studies consider the identification of soil hydraulic parameters using soil water retention data and a transient multistep outflow experiment (MSO), whereas the third study involves the calibration of a conceptual rainfall-runoff model.

Bayesian analysis is widely used in science and engineering for real-time forecasting, decision making, and to help unravel the processes that explain the observed data. These data are some deterministic and/or stochastic transformations of the underlying parameters. A key task is then to summarize the posterior distribution of these parameters. When models become too difficult to analyze analytically, Monte Carlo methods can be used to approximate the target distribution. Of these, Markov chain Monte Carlo (MCMC) methods are particularly powerful. Such methods generate a random walk through the parameter space and, under strict conditions of reversibility and ergodicity, will successively visit solutions with frequency proportional to the underlying target density. This requires a proposal distribution that generates candidate solutions starting from an arbitrary initial state. The speed of the sampled chains converging to the target distribution deteriorates rapidly, however, with increasing parameter dimensionality. In this paper, we introduce a new proposal distribution that enhances significantly the efficiency of MCMC simulation for highly parameterized models. This proposal distribution exploits the cross covariance of model parameters, measurements, and model outputs and generates candidate states much alike the analysis step in the Kalman filter. We embed the Kalman-inspired proposal distribution in the DiffeRential Evolution Adaptive Metropolis algorithm during burn-in and present several numerical experiments with complex, high-dimensional, or multimodal target distributions. Results demonstrate that this new proposal distribution can greatly improve simulation efficiency of MCMC. Specifically, we observe a speedup on the order of 10-30 times for groundwater models with more than 100 parameters. [ABSTRACT FROM AUTHOR]

Many watershed models used within the hydrologic research community assume (by default) stationary conditions, that is, the key watershed properties that control water flow are considered to be time invariant. This assumption is rather convenient and pragmatic and opens up the wide arsenal of (multivariate) statistical and nonlinear optimization methods for inference of the (temporally fixed) model parameters. Several contributions to the hydrologic literature have brought into question the continued usefulness of this stationary paradigm for hydrologic modeling. This paper builds on the likelihood-free diagnostics approach of Vrugt and Sadegh (2013) and uses a diverse set of hydrologic summary metrics to test the stationary hypothesis and detect changes in the watersheds response to hydroclimatic forcing. Models with fixed parameter values cannot simulate adequately temporal variations in the summary statistics of the observed catchment data, and consequently, the DREAM(ABC) algorithm cannot find solutions that sufficiently honor the observed metrics. We demonstrate that the presented methodology is able to differentiate successfully between watersheds that are classified as stationary and those that have undergone significant changes in land use, urbanization, and/or hydroclimatic conditions, and thus are deemed nonstationary. [ABSTRACT FROM AUTHOR]

The quest for a more powerful method for model evaluation has inspired Vrugt and Sadegh (2013) to introduce 'likelihood-free' inference as vehicle for diagnostic model evaluation. This class of methods is also referred to as Approximate Bayesian Computation (ABC) and relaxes the need for a residual-based likelihood function in favor of one or multiple different summary statistics that exhibit superior diagnostic power. Here we propose several methodological improvements over commonly used ABC sampling methods to permit inference of complex system models. Our methodology entitled DREAM(ABC) uses the DiffeRential Evolution Adaptive Metropolis algorithm as its main building block and takes advantage of a continuous fitness function to efficiently explore the behavioral model space. Three case studies demonstrate that DREAM(ABC) is at least an order of magnitude more efficient than commonly used ABC sampling methods for more complex models. DREAM(ABC) is also more amenable to distributed, multi-processor, implementation, a prerequisite to diagnostic inference of CPU-intensive system models. [ABSTRACT FROM AUTHOR]

Spatially distributed hydrologic models are increasingly being used to study and predict soil moisture flow, groundwater recharge, surface runoff, and river discharge. The usefulness and applicability of such complex models is increasingly held back by the potentially many hundreds (thousands) of parameters that require calibration against some historical record of data. The current generation of search and optimization algorithms is typically not powerful enough to deal with a very large number of variables and summarize parameter and predictive uncertainty. We have previously presented a general-purpose Markov chain Monte Carlo (MCMC) algorithm for Bayesian inference of the posterior probability density function of hydrologic model parameters. This method, entitled differential evolution adaptive Metropolis (DREAM), runs multiple different Markov chains in parallel and uses a discrete proposal distribution to evolve the sampler to the posterior distribution. The DREAM approach maintains detailed balance and shows excellent performance on complex, multimodal search problems. Here we present our latest algorithmic developments and introduce MT-DREAM(ZS), which combines the strengths of multiple-try sampling, snooker updating, and sampling from an archive of past states. This new code is especially designed to solve high-dimensional search problems and receives particularly spectacular performance improvement over other adaptive MCMC approaches when using distributed computing. Four different case studies with increasing dimensionality up to 241 parameters are used to illustrate the advantages of MT-DREAM(ZS). [ABSTRACT FROM AUTHOR]

The current study evaluates the impacts of various sources of uncertainty involved in hydrologic modeling on parameter behavior and regionalization utilizing different Bayesian likelihood functions and the Differential Evolution Adaptive Metropolis (DREAM) algorithm. The developed likelihood functions differ in their underlying assumptions and treatment of error sources. We apply the developed method to a snow accumulation and ablation model (National Weather Service SNOW17) and generate parameter ensembles to predict snow water equivalent (SWE). Observational data include precipitation and air temperature forcing along with SWE measurements from 24 sites with diverse hydroclimatic characteristics. A multiple linear regression model is used to construct regionalization relationships between model parameters and site characteristics. Results indicate that model structural uncertainty has the largest influence on SNOW17 parameter behavior. Precipitation uncertainty is the second largest source of uncertainty, showing greater impact at wetter sites. Measurement uncertainty in SWE tends to have little impact on the final model parameters and resulting SWE predictions. Considering all sources of uncertainty, parameters related to air temperature and snowfall fraction exhibit the strongest correlations to site characteristics. Parameters related to the length of the melting period also show high correlation to site characteristics. Finally, model structural uncertainty and precipitation uncertainty dramatically alter parameter regionalization relationships in comparison to cases where only uncertainty in model parameters or output measurements is considered. Our results demonstrate that accurate treatment of forcing, parameter, model structural, and calibration data errors is critical for deriving robust regionalization relationships. [ABSTRACT FROM AUTHOR]