VISCOUS: A Variance‐Based Sensitivity Analysis Using Copulas for Efficient Identification of Dominant Hydrological Processes

Global sensitivity analysis (GSA) has long been recognized as an indispensable tool for model analysis. GSA has been extensively used for model simplification, identifiability analysis, and diagnostic tests. Nevertheless, computationally efficient methodologies are needed for GSA, not only to reduce the computational overhead, but also to improve the quality and robustness of the results. This is especially the case for process‐based hydrologic models, as their simulation time typically exceeds the computational resources available for a comprehensive GSA. To overcome this computational barrier, we propose a data‐driven method called VISCOUS, variance‐based sensitivity analysis using copulas. VISCOUS uses Gaussian mixture copulas to approximate the joint probability density function of a given set of input‐output pairs for estimating the variance‐based sensitivity indices. Our method identifies dominant hydrologic factors by recycling existing input‐output data, and thus can deal with arbitrary sample sets drawn from the input‐output space. We used two hydrologic models of increasing complexity (HBV and VIC) to assess the performance of VISCOUS. Our results confirm that VISCOUS and the conventional variance‐based method can detect similar important and unimportant factors. Furthermore, the VISCOUS method can substantially reduce the computational cost required for sensitivity analysis. Our proposed method is particularly useful for process‐based models with many uncertain parameters, large domain size, and high spatial and temporal resolution.

need to represent nature, leading to a higher number of physical processes included with greater degree of sophistication (Clark et al., 2017;Fatichi et al., 2016). Some observers suggest a sociological element where sophisticated process representations are linked to modelers' hubris . The increased complexity ultimately results in modifying the model structure through including new model components, for example, by adding extra parameters and new feedbacks, which in turn requires carefully evaluating interacting decisions in the model building process (Gharari et al., 2021).
This high level of complexity in hydrological models can inevitably cause a high level of uncertainty, which needs to be quantified to glean useful information about the system's behavior and make robust decisions (Clark et al., 2008;Gupta et al., 2008;Hall, 2007;Refsgaard et al., 2007). A comprehensive analysis of the behavior of such models can shed light on how an output of interest represents a distinct aspect of the hydrologic cycle and can provide an informative view of the full spectrum of the system's behavior under different conditions. This is of vital importance, particularly when coping with a range of decision-making problems, such as investigating the system's response to a diverse range of plausible scenarios (e.g., future hydroclimatic conditions).
In this regard, global sensitivity analysis (GSA) (Saltelli et al., 1993(Saltelli et al., , 2008 is a powerful tool for model analysis. GSA has been broadly applied to a wide range of modeling problems, including model calibration, identifiability analysis, model simplification, and diagnostic testing (e.g., Guse et al., 2016;Markstrom et al., 2016;Rakovec et al., 2014). GSA is well recognized to be an essential means of solving some of the outstanding challenges associated with complexity of the hydrologic models, including but not limited to (Razavi et al., 2021): (a) determining importance of the input factors, which is helpful for factor prioritization and data acquisition; (b) detecting non-influential factors, whose variations have no impact on the behavior of the model and can be fixed to reduce model complexity; and (c) diagnosing differences between models, hypotheses, and constraints.
A comprehensive GSA of a hydrologic model typically requires many model runs, and hence considerable computational resources. This is mainly due to the high number of interacting, and uncertain input factors causing the "curse of dimensionality" problem. When using a sampling-based strategy to perform GSA (e.g., varying all parameters except one parameter at multiple locations in the parameter space), the curse of dimensionality along with the highly nonlinear nature of the response surface in hydrologic models become two major obstacles for efficient GSA. In other words, adequate exploration of the factor space and characterizing nonlinearity of the response surface necessitates a large sample size in traditional sampling-based GSA, which results in a considerable number of model evaluations.
The problem of high computational demand incurred by the sampling-based GSA manifests itself in the fact that many of the recent GSA studies in environmental modeling and hydrology (approximately 70%) have used relatively low-dimensional models . In contrast, the number of factors in complex, state-of-the-art hydrologic models is relatively high. More importantly, due to the high number of model evaluations required by existing GSA techniques and the computationally expensive nature of the hydrologic models, analysts usually tend to conduct GSA without evaluating its stability and convergence (for more discussion see, e.g., Cosenza et al., 2013;Ferretti et al., 2016;Mai & Tolson, 2019;Nossent et al., 2011;Sarrazin et al., 2016). It is therefore common to choose the sample size (the number of model runs) only based on the available computational budget, which in turn can result in lack of robustness, and consequently contaminate the assessment of the sensitivities (Hill et al., 2016). Therefore, devising strategies to minimize the computational cost of the GSA techniques is one of the major challenges in the quest to better understand the behavior of hydrologic models.
In the next subsection, we critically review three cost-effective strategies that have been adopted in the literature to accelerate GSA of the computationally expensive models.

Improved Sampling Strategies for GSA
Sampling from the input/factor space is a cornerstone of the sampling-based model analysis methods. In sampling-based GSA methods, an initial set of sample points needs to be drawn from the factor space in SHEIKHOLESLAMI ET AL. 10.1029/2020WR028435 2 of 24 order to extract the sensitivity-related information. However, many properly distributed points are usually necessary to sample the high dimensional factor space for analyzing the nonlinear response surfaces. Moreover, the arrangement of sample points in parameter space are often method-dependent (i.e., ad-hoc experimental designs), meaning that the parameter samples used in one GSA method cannot be applied in another GSA method. This issue can make sampling-based GSA very time-consuming, if not unattainable, particularly when computationally expensive models are used. However, the number of model evaluations needed for a comprehensive GSA can be effectively reduced if the utilized sampling strategy conveys the maximum amount of information from the output space with minimum sample size (Andres, 1997). A number of studies have shown that utilizing an appropriate sampling plan for GSA can decrease the number of model runs several orders of magnitude (Castaings et al., 2012;Gan et al., 2014;Gong et al., 2016).
The efficiency of any sampling-based GSA can be improved by using a sampling strategy that limits the computational load by optimizing some basic sample quality merits (Crombecq et al., 2011). This is also known as the "optimal design" of GSA methods. An optimal design should contain sample points that are uniformly distributed over the entire factor space such that all regions of the factor space are equally explored (this criterion is referred as space-fillingness) (Pronzato & Müller, 2012;. As shown by Janouchová and Kučerová (2013), a sample set with poor space-filling properties can cause large errors when estimating sensitivity measures. Additionally, there are alternative sampling methods that were only designed for specific GSA methods to reduce the computational cost (see, e.g., Chan et al., 2000;Saltelli, 2002;Campolongo et al., 2007;Ruano et al., 2012;Razavi & Gupta, 2016). In a comprehensive study, Piano et al. (2021) investigated the computational efficiency of the several sample-based estimators for variance-based sensitivity indices. One of the main drawbacks of the aforementioned sampling strategies is specifying the sample size a priori, and thereby increasing the risk of unnecessary computational burden (Sheikholeslami & Razavi, 2017).

Emulation Techniques for GSA
Constructing a cheaper-to-run emulator to replace the original computationally expensive model is another common approach to lower the computational cost in GSA (Ratto et al., 2012). Given a set of pre-existing model runs as a "training data set," this strategy first fits a statistical model (called emulator or surrogate model) to approximate the relationship between input factors and model output. Then, sensitivity measures are estimated using the learned emulator via either Monte Carlo-based methods or analytical approaches. Examples of the Monte Carlo-based emulators for GSA include neural networks (Marseguerra et al., 2003) and polynomial regression (Iooss et al., 2006) (for a comprehensive comparison of various Monte Carlo-based emulation techniques see Storlie et al. (2009) andVerrelst et al. (2016)). On the contrary, there are emulators that allow to estimate sensitivity measures analytically from the fitted surrogate model, including Gaussian radial basis functions (Chen et al., 2005), polynomial chaos expansions (Sudret, 2008), Gaussian process model (Bounceur et al., 2015), complex linear model (Jourdan, 2012), and sparse polynomial chaos (Wu et al., 2020). Since variance-based sensitivity analysis methods are an established good GSA practice, many of these emulation-based algorithms have been only developed for estimating variance-based sensitivity indices. Note that unlike the above-mentioned emulators which fit the (static) model output, dynamic emulation modeling provides a low-order emulator that preserves the dynamical nature of the original model. Monte Carlo-based approach is often used to identify dominant modes of dynamic behavior in dynamic emulation modeling (Young, 1999).
Because emulators are only an approximation of the original model, the emulation-based GSA results are typically prone to two types of uncertainties (O'Hagan, 2006): (a) the code uncertainty and (b) the uncertainty due to finite sample size. The former is associated with the emulator's degree of accuracy, while the latter is due to not knowing the output of the model at unsampled regions outside of the training data set. Therefore, quantifying uncertainty of the emulation-based GSA results is necessary and can function as a quality check. Another eminent challenge is that these emulation techniques usually work well when the problem at hand has a relatively low number of factors. But when the dimensionality of the problem is reasonably high, as in more complex models, emulators can become practically unsuitable for finding influential factors due to the curse of dimensionality and over-fitting issues (Becker et al., 2018). Emulation-based GSA suffers from two more drawbacks: (a) choosing the best emulator may not be easy due to the existence SHEIKHOLESLAMI ET AL.

10.1029/2020WR028435
3 of 24 of a plethora of different emulators, and (b) some emulators also need ad-hoc experimental designs and specific arrangement of sample points for the sake of efficiency.

A Given-Data Approach to GSA
An important disadvantage of the GSA algorithms cited so far is their dependency on a specific sampling design. This means that available data from previous model runs (generated for other modeling purposes, e.g., for uncertainty propagation, model emulation or calibration) cannot be re-used for GSA. To overcome this issue, the given-data approach (otherwise known as data-driven method) has been introduced. The given-data principle states that the sensitivity-related information can be extracted from a finite (generic) sample of an existing data set, independent of having a model and/or without re-running the model (Borgonovo et al., 2017). As argued by Sheikholeslami and Razavi (2020), the given-data approach can potentially advance GSA paradigm on two fronts. First, it proposes a new research area, where GSA can be applied to any data set, even when the underlying relationships/mechanisms are unknown. Second, it makes GSA tractable for computationally expensive models, as it can produce robust GSA results with minimal computational cost.
In this context, Plischke (2010); Strong et al. (2012); Strong and Oakley (2013); Wainwright et al. (2014); Li and Mahadevan (2016) presented computationally efficient algorithms for calculating variance-based sensitivity indices. Apart from variance-based GSA methods, moment-independent methods, which do not rely on any specific moment, can also work with a given sample of input-output pairs (see, e.g., Pianosi & Wagener, 2018;Plischke et al., 2013). Basically, these methods (also referred to as distribution-based methods) quantify the difference between conditional and unconditional density functions using distance-based metrics to measure the importance of input factors. Another widely applied method is the regional sensitivity analysis of Hornberger and Spear (1981) which does not require a specific sampling design for GSA. This method measures the sensitivity of outputs with respect to input factors by dividing model outputs into several classes, and then comparing the conditional distribution of each input factor between classes. More recently, Sheikholeslami and Razavi (2020) introduced a new given-data estimator for GSA based on theory of variogram analysis of response surfaces.
The rationale for given-data approach is as follows: If the sample of input-output data contains enough information to determine the dominant factors that exert strong influence on system's behavior, then the given-data estimator can be utilized as a post-processing GSA module for the pre-computed model evaluations. Using this approach, given a sample of input-output data representing the behavior of the system, which might be obtained from observations, laboratory, or field experiments, we can directly calculate the sensitivity measures from a generic sample set, without re-running the model. Note that due to a large factor space of the high-dimensional models, generating "representative" input-output data for GSA may require many properly distributed sample points (or observations) to explore a full spectrum of the system's behavior. Given this, regardless of the chosen method for GSA from given data, it is beneficial to spend some time up front to find an optimal sample set by prudently using sampling algorithms before performing GSA.
Note that compared to aforementioned cost-effective techniques (i.e., improved sampling strategies and emulators that have been utilized to enable computationally tractable GSA) the given-data approach has several advantages. The given-data approach mainly focuses on situations where analysts want to conduct GSA retrospectively, that is, using data from previous experiments, such as calibration or uncertainty propagation. The computationally expensive models will still need to be run, but not using specific sampling strategies. Moreover, the given-data approach does not rely on emulators. As mentioned earlier, emulation-based GSA techniques have major shortcomings that render them less useful compared to given-data approach.

Objectives and Scope
The primary goal of this paper is to develop a data-driven variance-based GSA method based on copula functions. The method is mainly designed to alleviate the high-computational demand associated with GSA of process-based hydrologic models. At present, unlike the sampling-based GSA methods, most of the given-data methods only focus on estimating the variance-based first-order sensitivity index (i.e., main effect) and cannot estimate the total sensitivity and interaction effects (see, e.g., Sparkman Strong et al., 2012). Hence, the secondary objective of this paper is to develop an efficient data-driven method to simultaneously compute first-order and total effect indices.
Our proposed GSA method (hereafter called VISCOUS, A VarIance-based Sensitivity analysis using COp-UlaS) utilizes a class of functions known as Gaussian mixture copulas (Tewari et al., 2011) to model the joint probability distribution of the inputs and the output, which might have multimodal densities and nonlinear dependencies. By approximating the conditional distributions of the output, when the inputs are conditioned to particular values, VISCOUS extracts sensitivity-related information directly from the given data, without augmenting the computational cost. This avoids the need for additional model runs and will help modelers to efficiently identify dominant factors, particularly when dealing with computationally expensive models. Furthermore, VISCOUS can be effectively used to rank the strength of relationships in multivariate data sets (not shown in this paper). In this paper, we use two real-world case studies, when only a fixed amount of input-output data is available, to demonstrate how well VISCOUS approximates variance-based sensitivity indices. The case studies include a conceptual rainfall-runoff model (HBV) and a physically based hydrologic model (VIC) configured for simulating the streamflow values of a headwater basin in the Canadian Rockies.
The remainder of this paper is organized as follows: In Section 2 we briefly review the conventional variance-based GSA algorithm, summarize the concepts of copulas and mixture models, present our copula-based data-driven GSA method (VISCOUS), and discuss its implementation details. In Section 3, we describe two hydrologic modeling case studies followed by explaining the setup for computational experiments. Section 4 applies the proposed GSA method to these case studies and reports numerical results and analyses. Finally, we conclude the paper in Section 5 and provide recommendations for future work.

Original Variance-Based GSA Method
We denote the output of a model    Y g X as a function of an m-dimensional input vector m is an uncertain factor of the model, such as model parameters, initial conditions, or boundary conditions. Here, we also assume   . g is a deterministic, scalar function of the input variables. In the context of variance-based GSA, uncertainty of the input vector X is often modeled in a probabilistic framework, that is X is deemed to be a random vector and uncertainty of the X propagates through   g X . Therefore, the model output uncertainty can be expressed in terms of its variance. Basically, in variance-based GSA the main goal is to apportion the model output variance into fractions attributed to each input factor and their interactions. Thus, we can quantify how much the variability of an individual (uncertain) input factor i X (and its interaction with other factors) contributes to the variability of the output. To do this, Sobol method (Sobol, 2001) uses a well-known variance decomposition formula (known as generalized Hoeffding-Sobol decomposition (Hoeffding, 1992)), according to Equation 1: In Equation 1, variance of the   g X is decomposed into the partial variances, where the first-order partial variance, i V , represents the contribution of an individual input factor i X to the output variance y V , that is: In other words, if we fix i X , the expected reduction in output variance will be i V . Moreover, ij V describes the portion of the output variance explained by the interaction between i X and j X , and so on for higher order interactions (  , ijk V ).
Variance-based GSA uses a set of indices to measure the importance of each input factor.  (Saltelli et al., 2019). The main and total effect variance-based sensitivity indices are often estimated using a sampling-based strategy, wherein the Monte Carlo integration technique is applied to Equations 2 and 4. Several numerical schemes for sampling-based strategy have been developed to compute these indices (Kucherenko & Song, 2017;Tarantola et al., 2006). In this study, a widely used and computationally affordable algorithm introduced by Saltelli (2002) and Saltelli et al. (2010) is utilized. We briefly describe this algorithm in the following, as it is the basis for the comparison of our method with the conventional variance-based method.
The sampling-based algorithm used in the present study starts by randomly generating two independent sample matrices: A and B, each of size  N m. Then, for calculating the i-th sensitivity index, an additional of size  N m is constructed by combining the columns of A and B, such that   i C contains all columns of B except the i-th column which is taken from A. Using the Jansen (1999)'s estimator, the sensitivity indices can be obtained from: where k A denotes the k-th row of the matrix A.

The Proposed Copula-Based Method for GSA From Given Data
Our proposed given-data approach (VISCOUS) is a probabilistic framework that estimates variance-based GSA indices from a given input-output sample. VISCOUS is a non-parametric approach and does not require generating new sample set for GSA, instead it recycles the existing data for determining factor importance. An essential building block of our proposed method is to infer the joint probability density function of the given data and to characterize its dependency structure using copula models. The theoretical background and implementation details of VISCOUS are provided in the following subsections.

Linking Copula Theory to Variance-Based GSA
Estimating the expected value of the model output conditioned on a factor,   E | i Y X , is the cornerstone of the variance-based method (see Equations 2 and 4). For a given set of model output values obtained from running a simulation model using a random sample of input factors, assume the marginal probability density function (PDF) along the i-th direction is y x denote the joint PDF of the model output Y and the individual input factor i X . Then, the conditional expectation can be written as (Beckman & McKay, 1987): y x accurately describes the joint PDF of the input-output sample, we can simply compute various sensitivity indices based on Equation 7. This has been motivated a number of studies to investigate several strategies for characterizing the joint PDF of the input-output data for the purpose of variance-based GSA (see, e.g., Cheng et al., 2015;Jia & Taflanidis, 2016;Sparkman et al., 2016;Wei et al., 2015).
In this paper, we rely on the copula theory which bridges one-dimensional marginal distributions to multi-dimensional joint distribution. Copula theory provides a powerful tool to model dependence structure between multivariate random variables from their associated marginals. Based on Sklar's (1959) theorem, the joint cumulative distribution function (CDF) of a bivariate random vector   , Y X can be expressed as a unique function of the marginal distributions: where  is a copula function, and X F and Y F are marginal CDFs of X and Y .
Hence, if the copula function and the marginal distributions are differentiable, the joint PDF of the input-output data for m-dimensional distribution function can be obtained as product of the individual marginal densities and the copula density: where c is the PDF of the copula function, There are several methods for density estimation (Silverman, 1986). The kernel density estimation is used in this paper to transform the input-output sample into   , v u . A crucial step when using kernel estimator is the choice of optimal bandwidth value. Very small values of the bandwidth lead to spiky estimates (not much smoothing) and overfits to the data whereas larger values can cause oversmoothing. In many cases finding an optimal bandwidth requires a good number of realizations.
If the forms of the marginal PDFs and the copula density are known, Equation 9 is utilized for estimating the unknown multivariate joint PDF. On the other hand, in case of the known multivariate joint PDF, the copula density can be derived by rearranging Equation 9 as: For example, in case of standard normal multivariate Gaussian data, Equation 10 can be solved analytically to obtain Gaussian copula density. The Gaussian copula has been previously used in conjunction with various GSA algorithms, including method of Morris (Ţene et al., 2018), variance-based method (Sainte-Marie & Cournède, 2019), and moment-independent method (Wei et al., 2014). However, the usefulness of the formal Gaussian copula can be severely limited when the response surface has multiple modes with nonlinear relationships. This issue is crucial for GSA of the hydrologic models, as it is not uncommon to encounter a highly complex and nonlinear response surfaces with multimodality, sharp discontinuities, and small-scale roughness in hydrologic models (Duan et al., 1992;Kavetski & Clark, 2010).
Alternatively, it has been shown that other types of copula functions such as empirical checkerboard copula (Deheuvels, 1979) and Gaussian mixture copula model (Li et al., 2011;Tewari et al., 2011) are more effective in modeling nonlinear relationships, particularly of non-Gaussian data (see, e.g., Bilgrau et al., 2016;Fan et al., 2016;Genest et al., 2017;Kasa et al., 2020). To capture the complex behavior of the hydrologic models, we developed an estimator based on Gaussian mixture copula model for constructing the joint distributions as discussed in the following sub-sections. Note that our proposed framework is not limited to the Gaussian mixture copula model and can be extended to any types of copula functions.

Gaussian Mixture Copula-Based Estimator for Sensitivity Indices
Gaussian mixture copula model (GMCM) combines advantages of the copula theory with capability of the Gaussian mixture models (GMM) in approximating an arbitrary PDF, wherein GMM is defined as a mixture statistical model of a finite number of Gaussian distributions. In fact, GMCM infers the dependence structure of the multivariate data through coupling GMM into copulas. Using Equation 10, GMCM estimates the PDF of the copula function by (Tewari et al., 2011): where  1 denotes the inverse function of the standard normal distribution, and f Z is the joint PDF of GMM, which can be expressed by the weighted sum of -component Gaussian densities,  , that is, ,  k is the weight for the k-th GMM component (note that all elements of the  k sum to unity),   k μ is the mean vector,   k Σ represents the covariance matrix for the k-th component, and comprises all the weights, mean vectors, and covariance After building c g using the GMCM scheme, the conditional expectation can be calculated by combining Equations 9 and 11 with Equation 7. Importantly, the integration in Equation 7 can be performed without re-running the original computer model (or physical model), which makes VISCOUS a computationally frugal method. We can use a Monte Carlo-based approximation to evaluate the integrals, according to Equations 13 and 14: where  c X x is the given random input factor of interest for which the Sobol indices should be estimated, MC N is the number of Monte Carlo samples for integration, and   1 y F v is the inverse CDF of Y which can be estimated using kernel method.

Steps for Performing VISCOUS
There are six steps involved in applying the proposed VISCOUS framework. In the following, we outline the algorithmic implementation of these steps in detail: Step 1 From the given data extract the c-th column of the input sample matrix (and the corresponding output values Y ) and set  c X x for which the variance-based sensitivity indices should be estimated.
Step 2 Use the kernel density function to convert the given data matrix into marginal CDFs, and set Step 3 Transform marginal CDFs of Step 2 into data of standard normal distribution and estimate the inverse of the standard normal distribution: ,..., .
Step 4 Build a mixture model f Z (i.e., GMM) by fitting  Gaussian density components to the data obtained from Step 3. If p is the number of indices that GSA needs to be performed, we build p copula models for calculating each index. For example, estimating   1 2 12 , , S S S requires constructing  3 p copula models. A critical step here is choosing optimal number of the components, , in Equation 12. To find the optimal value of , we utilized a likelihood-based goodness-of-fit statistic known as Bayes information criterion (BIC). BIC compares multiple models with varying values of  and chooses the optimal number of components associated with the lowest BIC value (Steele & Raftery, 2010).
Step 5 Fit a copula c g (Equation 11) using the data of Step 3 and f Z obtained from Step 4. We employed the expectation-maximization (EM) algorithm of Tewari et al. (2011) for estimating GMCM's parameters,   . Starting from initial values for   , the EM algorithm iteratively attempts to maximize the likelihood function. To avoid the risk of trapping into local optima, we used several starting points within an iterative loop for global optimum investigation. To do this, we applied k-means clustering strategy for initialization of the EM algorithm (for more details see McLachlan and Peel (2000)). A similar strategy was implemented in the R-package of Bilgrau et al. (2016) and fitgmdist function of the MATLAB for parameter estimation.
Step 6 Compute the variance-based sensitivity indices using Equations 13 and 14 based on the GMCM model learned from Step 5.
It is noteworthy that there are two sources of error in VISCOUS. First, error resulting from the use of numerical integration in Step 6. To reduce this error, MC N should be a large number (we arbitrarily set  4 MC 10 N ). Second, the kernel density estimator in Step 2 introduces an additional layer of error when transforming the original data into the data of CDF values. In Appendix, we exemplify the use of VISCOUS and demonstrate its efficiency and effectiveness by focusing on an analytical test function adopted from Kucherenko et al. (2011).
We conclude this section by highlighting the key distinctions between our proposed given-data approach and emulation techniques for GSA. First, the latter is mainly designed to approximate an output value corresponding to an individual realization of the uncertain input factors (see Section 1.2.2), while the former approximates a distribution of the output given probability distributions of the input factors. Particularly, for a given set of input-output data, VISCOUS trains GMCM in the probability space, but emulators are always constructed in the output space. Second, to perform GSA and assess the output uncertainty associated with uncertainty of a given input factor, emulators typically need to be evaluated many times for multiple realizations of that uncertain factor. On the other hand, one run of our approach yields the conditional probability distributions of the output without repeated runs of the model (or surrogate model).
Another major difference is that the dimensionality of the joint PDF obtained by GMCM in VISCOUS is independent of the dimensionality of the original model. This makes VISCOUS a suitable candidate in case SHEIKHOLESLAMI ET AL.
10.1029/2020WR028435 9 of 24 of high-dimensional models for which fitting an emulator often suffers from the curse of dimensionality (see Section 1.2.2). In VISCOUS, the dimensionality of the joint PDF is  1 s , where s is the order of the variance-based GSA indices. By way of example, the joint PDF is only two-dimensional when computing the first-order indices.

Study Site
In this study, we set up two hydrologic models to simulate streamflow for the Bow River at Banff, Alberta, Canada ( Figure 1a) with a basin area of approximately 2,210 km 2 . The headwaters of the Bow River are in the eastern Front Range of the Canadian Rocky Mountains. The temperature in this part of the Rockies can dip below −30°C in winter on mountain tops and rise up to +30°C in summer in Bow Valley (basin average temperatures are depicted in Figure 2). Most parts of the basin receive a large proportion of their precipitation as snowfall (∼50%), which results in substantial water storage as snowpack. Snowmelt (beginning in mid-April) is the dominant process that strongly controls the dynamics and generation of the runoff in the Bow River Basin (Nivo-glacial regime). The average basin elevation is 2,130 m ranging from 1,380 m above mean sea level at the stream gauge (located in the town of Banff) to above 3,000 m at the highest SHEIKHOLESLAMI ET AL.   (2017)). The predominant land cover (Figure 1c) is conifer forest in the Bow Valley, and rocks and gravels for mountain peaks above the tree line.

A Conceptual Rainfall-Runoff Model: HBV
We first illustrate the utility of the proposed GSA algorithm using the HBV model, which is a lumped, conceptual hydrologic model (Bergström, 1995;Seibert, 1997). HBV simulates daily streamflow time series using input data of mean temperature, precipitation (daily values) and potential evapotranspiration (monthly estimates). The model consists of three main components: (a) snow module, (b) soil moisture module, and (c) flow routing module. The latter is composed of two reservoirs for quick flow (non-linear reservoir) and slow flow generations (linear reservoir). HBV transforms the simulated outflow using a triangular weighting function for channel routing. In this study, we implemented a slightly modified version of the HBV model developed by Razavi et al. (2019). The total number of calibration parameters was set to 10. Table 1 SHEIKHOLESLAMI ET AL.

Parameter [unit]
Range Definition summarizes parameter description, units of measurements, and their feasible ranges of variation. Note that the feasible ranges of these parameters were determined using a combination of expert knowledge and previous studies.

A Semi-Distributed Physically Based Hydrologic Model: VIC
As the second case study we utilized a more computationally expensive, process-based model, known as variable infiltration capacity (VIC). The VIC model is developed as a simple land model (Liang et al., 1994) which simulates both mass and energy fluxes. The model setup used in this study has 18 parameters and is based on the concept of flexible vector-based configuration for land models (for more details on the model set up refer to configuration case-2-4 km in Gharari et al., 2020). Description of the model parameters are presented in Table 2. Parameter ranges were specified using a combination of expert knowledge and previous studies, and recommendations by the developers.
In our implementation, the VIC model considers soil type and vegetation, and is forced with 7 meteorological variables from WRF simulations with spatial resolution of 4 km (Rasmussen & Liu, 2017). In the VIC model, the saturation excess mechanism generates all fast response runoff and a gravity-driven process controls the drainage between the three soil layers of each grid cell (Liang et al., 1994). The soil column is modeled using three layers, wherein the upper two layers are responsible for the surface flow generation, called runoff, and the third layer generates the slow response, called baseflow. In this study, the traditional 4-parameter baseflow layer is replaced with a simplified 1-parameter linear layer (for further implementation details refer to Gharari et al., 2019).
Furthermore, the streamflow was simulated using the vector-based routing tool mizuRoute  that uses a Gamma distribution-based unit hydrograph to route runoff from hillslope to the river channel (also known as in-basin routing). Moreover, mizuRoute has an in-channel impulse response function routing scheme based on diffusive wave along the river segment. Note that four mizuRoute parameters (i.e., parameters 15-18 in Table 2) were also considered for perturbation in our experiment.

Selecting Output Function
Residual-based summary statistics, such as Nash-Sutcliffe efficiency criterion ( NS e ), are common objective functions in hydrology and have been widely used for sensitivity analysis to investigate which parameters significantly influence the model predictions. This provides valuable information for parameter estimation and model calibration methods. Therefore, we calculated the NS e in our experiments and considered it as the response function for GSA: where T is the number of time steps, obs t Q is the observed streamflow values at time t, sim t Q is the simulated value of the streamflow at time t, and Q is the mean value of the observed streamflow values. Here, we compared the simulated flows with the observed streamflow data of the Water Survey of Canada (gauge ID of 05BB001).

Generating Synthetic Data Sets for VISCOUS Runs
To evaluate the effectiveness of the proposed VISCOUS method and ensure a meaningful comparison, the GSA results obtained by the sampling-based algorithm were deemed to be the "true" sensitivities. Therefore, we first ran the conventional sampling-based GSA algorithm for both case studies using the method described in Section 2.1. For the first case study, the HBV model was run using a large sample size to ensure the stability and convergence of the sampling-based algorithm, resulting in 1,000,008 evaluations. The computational efficiency of the HBV model made it feasible to select this large sample size, and consequently enabled us to benchmark our results against the true sensitivities. For the second case study, we focused on the GSA of the VIC model using a data set consisting of 21,000 realizations of the input parameters and the corresponding NS e values.
In the next step, we produced synthetic input-output data sets by randomly sampling (with progressively increasing sample size) from the available model evaluations (1,000,008 and 21,000 input-output pairs obtained from HBV and VIC runs, respectively). Each input-output sample set was considered as a set of "given data." Then, we applied our VISCOUS method to all the given sample sets and compared the results with true sensitivities.

Accounting for Randomness in Sampling Variability
Because of the randomness in our experiments, the performance of the GSA methods can be quite sensitive to the choice of the sample points, that is, the results can be afflicted by the randomness in the selection of sample points drawn from the input-output space. This necessitates the use of several trials with different random seeds for each experiment to obtain a more reliable measure of the algorithm performance. Therefore, we replicated all GSA experiments 50 times in this study, each replication using different random seed. This allowed us to see a range of possible performances for the proposed algorithm, and thus evaluate its robustness. Note that the progressive Latin hypercube sampling strategy (Sheikholeslami & Razavi, 2017) was used to randomly generate parameter sets from uniform distributions. In practice, however, the bootstrap method can be used to avoid a large number of simulations when assessing the robustness and convergence of the GSA results in computationally expensive models. Here, a high-  Figures 3 and 4 show the main effect indices ( i S ) for the 10-parameter HBV model and the 18-parameter VIC model estimated by the conventional variance-based and VISCOUS methods. For the conventional variance-based method, the sampling-based algorithm described in Section 2.1 was applied; for the VISCOUS method, Equations 13 and 14 were applied to the given set of input-output data. In Figure 3, calculation of the (final) true sensitivity indices (the dashed horizontal lines) required about 1 million model evaluations using the sampling-based method, while the proposed VISCOUS method has been applied by increasing the size of input-output data incrementally up to 50,000 (i.e., 1,000, 2,000, …, 50,000). For the VIC model, the VISCOUS-based estimates of the main effect indices depicted in Figure 4 were also obtained by randomly sampling with progressively increasing size (i.e., 1,000, 1,400, 1,800, …, 21,000) from the available pre-computed model evaluations. In addition, Figures 3 and 4 show the 90% intervals of the estimated sensitivity SHEIKHOLESLAMI ET AL.  indices obtained from 50 replicates of the experiments. Overall, these results suggest that both VISCOUS and the conventional variance-based methods yielded similar main effect indices across the parameter space, thereby confirming the accuracy of our proposed method in estimating the sensitivity indices for HBV and VIC.

Results for the Main Effects
To assess the effect of sampling uncertainty and evaluate the robustness of VISCOUS to sampling variability, the 5th and 95th percentiles (90% intervals) of the 50 realizations of each experiment are presented in Figures 3 and 4. For HBV (see Figure 3), the results show that our proposed method converged to the true sensitivity indices within the 90% interval when the samples size was larger than about 8,000. In addition, by analyzing Figure 3 we see that sensitivity indices associated with three important parameters of HBV (FRAC, C0, and beta in Figure 3) show a high variability compared to other parameters when the sample size was smaller than 8,000. For the VIC model (see Figure 4), the 90% intervals indicate that most of the sensitivity indices reached a reasonable convergence when the sample size was larger than about 13,000. Note that width of the 90% intervals of the sensitivity indices in the case of VIC are quite narrow for the sample seizes greater than 17,000.
It is worth mentioning that three HBV parameters with lower sensitivity (alpha, EFT, and K1 in Figure 3) have not converged to true indices (slightly overestimated) even when the sample size was 50,000. The main SHEIKHOLESLAMI ET AL.  effect values associated with these less important parameters are very small (close to zero), and therefore overestimation of the sensitivity indices is mainly because of the numerical errors in Equations 13 and 14. Like HBV, some of the less important parameters of VIC exhibit poor convergence behavior in Figure 4 (i.e., d 1rnf , S stomatal , and V). If the goal of conducting GSA is finding the accurate estimates of these sensitivity indices, a larger sample size should be used to achieve a higher degree of accuracy. However, when dealing with computationally expensive models, estimating accurate values of the sensitivity indices is often not of interest. In practice, one might focus on determining if a parameter belongs to a high, medium, or low-influence group rather than on its exact ranking, particularly considering that parameter ranking usually converges faster than sensitivity indices (Sarrazin et al., 2016;Sheikholeslami, Razavi, Gupta, et al., 2019).

Results for the Interaction Effects
One straightforward approach to quantify the interaction effects in variance-based GSA is to calculate the sum of main effects (Borgonovo et al., 2017). The sum of main effects explains to what extent model parameters are individually important, while the remainder to one implies the degree to which the interaction among parameters is influential, that is, the remainder      1 1 m i i S is a fraction of the output variance because of the interactions among all model parameters (see Equation 3). Based on our results for HBV, the sum of main effect indices is about 0.80 (80% of the variance), suggesting that interactions (expressed by total effect index, Ti S in Equation 4) have very limited relevance for this case study. Hence, the HBV response (i.e., NS e values) can be approximately considered as an additive function of the parameters over their feasible ranges.
However, for the VIC model the sum of main effect indices is about 0.50 (50% of the variance), indicating that interactions play a significant role in the variability of the VIC model output. As a result, in this section, we opted for presenting the total effect results only for VIC. To investigate the interaction effects, we calculated total effect sensitivity indices, Ti S , using the proposed VISCOUS algorithm. Figure 5 compares the total effect and main effect indices of an arbitrarily chosen realization out of 50 independent realizations of this experiment for the VIC model. As can be seen, most of the parameters show a significant total effect compared to their main effect values.
Figures 5a and 5b highlight a clear difference between total effect and main effect indices. We closely examine this difference in Figure 5c. As shown, most of the parameters lie further toward the top-left corner of Figure 5c, confirming the strong interactions between parameters of the VIC. This is an important sign of a possible non-identifiability issue in our model. In fact, it has been asserted that parameters associated with low main effects but having large total effects are typically less identifiable (see, e.g., Borgonovo et al. (2017);Hill et al. (2016); Ratto et al. (2007)). For example, Hill and Tiedeman (2007) reported a same observation for their groundwater problem with dominant interaction effects. On the other hand, if parameters lie close to the ideal (1:1) line (red line in Figure 5c), they are likely identifiable. Note that having small main effects but high total effects does not necessarily indicate that parameters are non-identifiable (i.e., insensitivity is not a sufficient condition for non-identifiability).
Figure 5c may also be an indication of a badly defined calibration problem, which typically corresponds to a highly over-parameterized model. In such models, non-identifiability problem can be mainly attributed to the negligible influence of most of the parameters and their associated physical processes on the NS e values (parameters with low sensitivity in Figure 5a). The non-identifiability can also occur due to the compensation among parameters, wherein changes in one parameter may be compensated for by changes in other parameters as evident by strong interaction between VIC parameters in Figure 5b. Overall, our results follow the identifiability analysis of Gharari et al. (2020). For example, they have also reported that two soil parameters (K sat and E exp ) were non-identifiable for the utilized configuration of VIC. As shown in Figure 5c, these parameters have very small main effects but high total effects, implying lack of identification.

High Importance Hydrological Processes Identified by the VISCOUS Method
To characterize dominant processes affecting model responses, we implemented the recently developed grouping-based importance ranking approach of Sheikholeslami, Razavi, Gupta, et al. (2019). This approach uses agglomerative hierarchical clustering to categorize the parameters into distinct groups based SHEIKHOLESLAMI ET AL.

10.1029/2020WR028435
16 of 24 on similarities between their sensitivity indices, and then ranks parameters according to importance group rather than individually. Figure 6 shows parameter grouping results of the HBV and the VIC models obtained from VISCOUS-based estimations of the main effect indices from 50 independent trials of the experiments.
As shown in Figure 6, we can categorize the HBV parameters into three importance groups according to their sensitivities: {FRAC, C0, K2} is the strongly influential group; {FC, TT, beta} is the moderately influential group; and {alpha, ETF, K1, LP} with the least influence on NS e . For VIC, parameters can be allocated into three groups as well: strongly influential parameters: {K slow , d 2f }; moderately influential parameters: {b inf , G scale , D diff }; and the rest of the parameters were ranked as the least influential parameters.
Factor importance ranking and grouping results in Figure 6 point at dominant hydrological processes that govern the streamflow generation mechanisms in the Bow River basin. Figure 6 reveals that slow reservoir coefficient in both models (K slow in VIC and K2 in HBV) is among the most important parameters influencing streamflow simulations. Note that we observe a strong yearly cycle for the Bow River due to the presence of snow and glaciers (see time series in Figure 2). In cold or mountainous regions, such as the Bow River where its headwaters are in the glaciated eastern slopes of the Canadian Rockies, snow melting and release of water from glacier storage greatly contribute to streamflow generation as opposed to rainfall. For example, Hopkinson and Young (1998)  dry years can be responsible for up to 50% of the late-summer flow in the Bow River at Banff. A significant contribution of these parameters to simulation of the streamflow can be attributed to the fact that K slow in VIC and K2 in HBV are responsible for the basin storage and timing of the release for getting the recession limb of the hydrograph correctly.
Soil moisture parameters in both models (b inf in VIC and beta in HBV) also play a key role in streamflow simulation (see moderately influential parameters in Figure 6). These parameters have strong impact on the partitioning of the precipitation into infiltration and runoff. An increase in b inf results in a lower infiltration capacity, and consequently yields higher runoff values simulated by the VIC model (Gharari et al., 2019). For HBV, a higher beta reduces runoff from soil, thereby increasing evaporation and reducing flow. Previous studies conducted by Gou et al. (2020); Melsen et al. (2016); Demaria et al. (2007) have also identified b inf as one of the important parameters of the VIC model.
Overall, based on our parameter grouping in Figure 6, it can be concluded that dominant factors of the HBV model mainly control the soil (FRAC and FC) and snowmelt (C0) processes. FRAC determines the fraction of soil release entering fast reservoir, as a result a higher FRAC value allows high flow of water from soil to the fast reservoir. FC accounts for partitioning of the precipitation into runoff and soil moisture. Hence, both FRAC and FC have a direct role on runoff generation in the basin. The high importance of C0 can be justified because the hydrological regime of the Bow River is highly influenced by the snow accumulation and melt, and the amount of snowmelt is linearly proportional to temperature via C0.
In the VIC model, d 2f is one of the most important parameters. An explanation might be the fact that d 2f together with S LAI are controlling the transpiration from plants. Also, notice the high total effect value associated with S LAI in Figure 5b. The transpiration process in VIC is based on the Jarvis's (1976) formulation, which is very sensitive to the leaf area index, S LAI , assigned to each vegetation type. This parameter together with d 2f can significantly affect the simulated evaporation from the forested areas (dominant vegetation type), and therefore adjust the transpiration as one of the major water balance components to get a closer streamflow to observation. Lastly, from Figure 6 we observe that two mizuRout routing parameters (G scale and D diff ) are highly important parameters of VIC. These parameters are mainly responsible for peak time and flashiness of the hydrograph.

Conclusions and Future Research
An outstanding challenge in GSA is to produce robust and stable results with low computational cost. In particular, the high computational cost incurred by the sampling-based GSA methods limits their application to complex process-based hydrologic models (Clark et al., 2017). To tackle this problem, we developed an effective data-driven algorithm for variance-based sensitivity analysis using copulas (VISCOUS). Our proposed strategy, VISCOUS, makes no assumptions on the underlying structure of the input-output data SHEIKHOLESLAMI ET AL.

10.1029/2020WR028435
18 of 24 wherein the probability distributions, correlations, and interactions between input factors and model output are learned only from available data. An intriguing possibility is that VISCOUS could provide more meaningful results than traditional applications of the computationally intensive, sampling-based Sobol method because the copula functions can effectively "smooth out" some of the numerical artefacts that are common in hydrological models (Kavetski & Clark, 2010).
Our primary motivations to develop VISCOUS include: 1. VISCOUS is a given-data approach. This means that it makes use of any collected data sets or pre-existing model runs, and accordingly does not require additional computational burden of new model simulations to satisfy the specific sampling strategies required by many existing GSA methods. Hence, VISCOUS enables modelers to efficiently conduct GSA for cases in which properties of the input-output distributions and of the underlying response surface are unknown and only a (sub-)sample of the input-output space is available. 2. VISCOUS opens up the possibility to incorporate information about the dependence structure between model parameters into the variance-based GSA. The dependencies can be specified by modeling their joint distribution. VISCOUS uses GMCM to approximate the joint probability density function of the input-output sample for estimating several variance-based sensitivity indices.
In this paper, we presented the numerical implementation and algorithmic details of our proposed given-data estimator, and examined its performance using two popular hydrologic models of increasing complexity (HBV and VIC). VISCOUS efficiently identified dominant parameters (and processes) of the VIC and HBV models that significantly influence prediction of the streamflow values in the Bow River Basin at Banff, Alberta, Canada. Overall, the numerical experiments indicate that: 1. Parameter sensitivities obtained by VISCOUS are similar to the conventional sampling-based method of Sobol, while our method requires lower computational demand. 2. VISCOUS is robust to sampling variability and randomness. The level of accuracy achieved by our proposed method is promising for practical purposes, such as for factor prioritization and factor ranking. 3. Parameter identifiability of the physically based distributed hydrologic models (VIC in the present study) can be troublesome and efficient GSA methods such as VISCOUS can serve as a computationally tractable tool for diagnosing lack of identification. 4. VISCOUS provides physically sound results regarding the sensitivity behavior of the parameters, even in cases where the sample size (or number of simulations) is limited. Our sensitivity analysis results correspond well with the expected behavior and dynamics of runoff generation in cold regions.
Future studies should include an extension of the proposed method to multi-criteria sensitivity analysis of the hydrologic models. Our data-driven approach, in principle, can be used in the setting of multi-criteria sensitivity analysis to efficiently estimate the sensitivity indices for models with multiple dependent outputs. One possible option is to combine VISCOUS with multiple-response emulators of Liu et al. (2019), which integrates the copula functions with multiple-response Gaussian process emulator. In addition, since the GMCM used in VISCOUS has a set of design parameters, a "sensitivity analysis of sensitivity analysis" (Puy et al., 2020) should be tried in future studies to investigate VISCOUS sensitivity to its own design parameters.
Our proposed method has a great potential to be applied to GSA of the new generation of complex hydrological models, such as the SUMMA framework (Clark et al., 2015a(Clark et al., , 2015b. VISCOUS can facilitate fast and effective GSA for computationally demanding models and can help achieve accurate and robust results. Finally, testing our method on real-world data sets obtained from field observations, remote sensing, or laboratory experiments is an interesting avenue for further research. This can help gain an in-depth insight into the underlying system's behavior, thereby better supporting model development and understanding.

Appendix A: Illustration Using an Analytical Test Function
We consider a type C benchmark problem for GSA introduced by Kucherenko et al. (2011), which has the following mathematical form: SHEIKHOLESLAMI ET AL. This function has dominant higher-order interaction terms with no unimportant subset of factors. Kucherenko et al. (2011) mathematically proved that for the function defined in Equation A1: where j S and Tj S are the first order and total order sensitivity indices for the j-th factor     1,2, , j d, respectively. The non-linearity and non-monotonicity of this function along with the availability of analytical results (Equation A2) make it a suitable problem for the study of GSA methods. As can be seen from Equation A1, this function is a product of contributions from each input factor, it is non-additive, and features interactions of all orders.
In this study, for illustrative purposes, the number of factors was set to  10 d , and the true values of 1 1 / T S S and  j S were used to evaluate the performance of our proposed approach. To implement VISCOUS, 100,000 points were randomly sampled from     10 0,1 using progressive Latin hypercube sampling strategy. In the next step, we produced synthetic input-output data sets by randomly sampling with progressively increasing sample size (i.e., 1,000, 2,000, …, 100,000) from the available 100,000 evaluations of the function. Each input-output sample set was considered to be a set of "given data." Then, we applied our VISCOUS method to all the given sample sets and compared the results with analytical values. In addition, we replicated this experiment 50 times, each replication using different random seed.  ing the accuracy of our proposed method in estimating sensitivity indices. Figure A1 also depicts the 5th and 95th percentiles (90% intervals) of the 50 realizations of each experiment. Note that width of the 90% intervals become narrower when the sample size increases.

Data Availability Statement
The WRF simulations are available from Rasmussen and Liu (2017). The HBV and VIC models are available through Razavi et al. (2019) and Gharari et al. (2020), respectively. The MATLAB codes for the proposed VISCOUS method will be available on GitHub upon publication of the manuscript (https://github.com/ Razi-Sheikh/VISCOUS). SHEIKHOLESLAMI ET AL.
10.1029/2020WR028435 21 of 24 Acknowledgments This paper was written during the COV-ID-19 pandemic which has transformed the world. We would like to take this opportunity to acknowledge a deep sense of gratitude to the workers on the frontline against COVID-19 around the world-healthcare workers, including doctors, nurses, paramedics, and others. We pay tribute to their courage and sacrifice and wish them rest and quiet from the storm. The authors would also like to thank editor D. Scott Mackay and associate editor Thorsten Wagener for handling this paper. We thank reviewers, Andrea Saltelli and three anonymous reviewers, for their insightful feedback on this work. The first author gratefully acknowledges Oxford University for covering open access publication costs.