Hierarchical Sensitivity Analysis for Large Scale Process-based Hydrological Modeling with Application in an Amazonian Watershed

Sensitivity analysis is an effective tool for identifying important uncertainty sources and improving model 15 calibration and predictions, especially for integrated systems with heterogeneous parameter inputs and complex processes coevolution. In this work, an advanced hierarchical global sensitivity analysis framework, which integrates the concept of variance-based global sensitivity analysis with hierarchical uncertainty framework, was implemented to quantitatively analyze several uncertainties of a three-dimensional, process-based hydrologic model (PAWS). The uncertainty sources considered include model parameters, model structures (with/without overland flow module), and climate forcing. We apply 20 the approach in a ~9000 km Amazon catchment modeled at 1 km resolution to provide a demonstration of multiple uncertainty source quantification using a large-scale process-based hydrologic model. The sensitivity indices are assessed based on three important hydrologic outputs: evapotranspiration (ET), ground evaporation (EG), and groundwater contribution to streamflow (QG). It is found that, in general, model parameters (especially those within the streamside model grid cells) are the most important uncertainty contributor for all sensitivity indices. In addition, the overland flow module 25 significantly contributes to model predictive uncertainty. These results can assist model calibration and provide modelers a better understanding of the general sources of uncertainty in predictions of complex hydrological systems in Amazonia. We demonstrated a pilot example for comprehensive global sensitivity analysis of large-scale complex hydrological and environmental models in this research. The hierarchical sensitivity analysis methodology used is mathematically rigorous and capable of being implemented into a variety of large-scale hydrological models with various sources of uncertainty. 30 https://doi.org/10.5194/hess-2019-246 Preprint. Discussion started: 2 August 2019 c © Author(s) 2019. CC BY 4.0 License.

In this study, we applied an advanced hierarchical global sensitivity analysis approach to a PBHM (the process-based adaptive watershed simulator; PAWS) for a watershed in the Amazon considering multiple uncertainty sources. Hierarchical sensitivity analysis was first proposed by Dai and Ye (2015) and then applied to a groundwater modeling analysis for the Hanford 300 Area in Washington, U.S. (Dai et al., 2017a). The hierarchical sensitivity analysis method integrates the 70 concept of variance-based global sensitivity analysis with the hierarchical uncertainty quantification framework (Draper et al., 1999;Dai and Ye, 2015) to quantify the sensitivity of important processes to model parametric and structural uncertainty.
This advanced sensitivity analysis method considers three important sources of uncertainty (i.e., scenario, structural, and parametric) involved in hydrological models and dramatically decreases computational costs by combining uncertain inputs based on their characteristics and interdependencies. We also improved the hierarchical sensitivity analysis methodology by 75 introducing new parameter groups into the uncertainty framework and implementing new algorithms to make the assessment of global sensitivity analysis for large-scale PBHMs computationally affordable. This study is the first to implement a comprehensive hierarchical sensitivity analysis method in relation to a complex and large-scale PBHM. In this study, we revised the hierarchical sensitivity analysis method and defined a set of sensitivity indices to group and accurately quantify the importance of different uncertainty sources in PAWS when simulating hydrological processes in an Amazonian 80 watershed.
PAWS has been applied extensively in many watersheds, e.g., the large watersheds in Michigan state, U.S. (Niu et al., 2014(Niu et al., , 2017Ji et al., 2015;Shen et al., 2013Shen et al., , 2014Shen et al., , 2016Qiu et al., 2019), and the watershed in the Amazon basin (Niu et al., 2017), and demonstrates high efficiency and good performance. Here, we applied PAWS to the Amazon because this region includes more than half of the tropical rainforests globally (Morley, 2000) and plays an essential role in the world carbon 85 (Richey et al., 2002;Phillips et al., 2009;Lintner et al., 2017) and water (Fearnside, 2005;Phillips et al., 2009) cycles.
Amazonian forest hydrology can affect many processes, including nutrient budgets (Lesack and Melack 1996), the production and consumption of carbohydrates (Pegoraro et al., 2006), and root-zone moisture availability for plants (Oliveira et al., 2005;Tang et al., 2015). A large number of field (Lesack 1993;Leopoldo et al., 1995;Tomasella et al., 2008) and numerical modeling (Coe et al., 2008;Sorribas et al., 2016;Yamazaki et al., 2012) studies have been performed to improve 90 our understanding of hydrologic processes in the Amazon. As some studies have indicated that groundwater is the key controller of Amazon hydrology (Leopoldo et al., 1995;Hodnett et al., 1997a,b;Beighley et al., 2009;Pokhrel et al., 2014;Niu et al., 2017), some researchers have already implemented simulations involving subsurface water in the Amazon (Miguez-Macho and Fan, 2012a,b). The severe hydrological hazards that frequently occurred in the Amazon in the past, such as the droughts in 2005 and 2010 and the floods of 2009 (Tomasella et al., 2008;Marengo et al., 2010;Espinoza et al., 2011), 95 also increase the requirement of including climate scenario uncertainty and related model uncertainty in sensitivity analyses for the Amazon region. 4 The PAWS model applied in this work was originally used in Niu et al. (2017) to simulate the hydrologic cycle in an Amazonian watershed. Here, we build on that work by classifying model system uncertainty sources into three groups: uncertainty related to climate scenarios (six climate scenarios that differ in terms of radiation, precipitation, temperature, 100 humidity, and wind speed), model structure (different thicknesses to represent aquifers), and parameters, including soil saturated hydraulic conductivity, Ks (m day -1 ), the Van Genuchten equation parameters α (m -1 ) and N (unitless) (van Genuchten, 1980), unconfined aquifer hydraulic conductivity, K1 (m day -1 ), confined aquifer hydraulic conductivity, K2 (m day -1 ), and length of flow path for runoff contribution to overland flow domain, L (m). We consider the Van Genuchten parameters α and N here because the correlation between α and N can largely affect the soil water release and infiltration 105 processes in the vadose zone (Pan et al., 2011). A new set of subdivided parametric sensitivity indices was first defined to provide more detailed information for parametric sensitivities. Because of the high complexity and dimensionality of this model, the highly efficient parameter sampling method of Latin hypercube sampling (LHS) and a binning method were applied to estimate the sensitivity indices. We also investigated the effects of prior weights on the climate scenarios and numerical models. By implementing the hierarchical sensitivity analysis method, we aim to provide a pilot example of 110 comprehensive global sensitivity analysis for large-scale PBHMs considering all uncertainty sources instead of only parameters and investigate the most important source of uncertainty for modeling hydrological processes in the Amazon.
We introduce the study area and the numerical model in Section 2.1. Sections 2.2 and 2.3 present the hierarchical sensitivity analysis method and its algorithms in detail. Then, we describe the definition and generation of uncertainty sources based on the study site information in Section 2.4. We present and discuss the results in Section 3. Section 4 summarizes the key 115 findings of this research.

Study site and numerical model
The study site is located in northern Manaus, Brazil (Fig. 1), with a drainage area of ~9,000 km 2 . Within the central Amazon, the watershed is mostly covered by tropical forest, with ~12% cropland and ~5% wetland (based on Community Land Model 120 (CLM) land surface data; Niu et al., 2017). With the relatively high elevation (90 -210 m) of the upper landscape and the relatively low elevation (45 -55 m) of the swampy valleys, a dense drainage network was formed in the region. The watershed has 4 rivers: the Urubu, Preto da Eva, Tarumã-aç u, and Tarumã -mirim rivers. The average precipitation in this region has large seasonal variability. From December to May is the wet season, while from June to November is the dry season. 125 The model tool used in this study is the PAWS model (Niu et al., 2017), which implements an effective non-iterative scheme to couple hydrologic processes including both land surface and subsurface water. The details of the numerical implementation and the governing equations of PAWS can be found in Appendix A. Briefly, four flow domains are https://doi.org/10.5194/hess-2020-87 Preprint. Discussion started: 10 March 2020 c Author(s) 2020. CC BY 4.0 License. 5 simulated in PAWS, including the stream channel, overland flow, the vadose zone, and saturated groundwater. The structured grid-based finite-volume method is the main numerical scheme applied to discretize the governing equations of 130 the various hydrologic components. PAWS also simulates two land surface subdomains, i.e., infiltration and evaporation, which are depicted in the ponding subdomain, while overland flow occurs in the surface flow subdomain. PAWS considers the horizontal interaction of both surface runoff and groundwater flow between the model grids, which represents the actual hydrological processes and is often ignored by other regional and global hydrologic models. The 1-D diffusive wave equation is solved to simulate channel flow, and the 2-D version is used for overland flow. The leakance concept is the basis 135 applied to explicitly simulate the exchange between the channel and groundwater. PAWS has been coupled with the CLM (Shen et al., 2013), which calculates the surface energy balance and soil and plant carbon and nitrogen cycles. Canopy interception and ET demand (both transpiration and soil evaporation) are also computed in the CLM at each time step.
For the numerical model case in this study, a 1 km × 1 km grid is used for horizontal discretization, resulting in 118 × 122 grid cells for the study site. In this model, 20 vertical layers were defined to discretize the vadose zone, and for the fully 140 saturated groundwater, there are two layers: the unconfined aquifer on the top and the confined aquifer at the bottom.
In this study, 90 m resolution NASA Shuttle Radar Topography Mission (SRTM) (U.S. Geological Survey; http://eros.usgs.gov) data are applied as DEM input, but for the channel network and watershed boundary delineation, the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) provides the 30 m resolution Global Digital Elevation Model Version 2 (GDEM V2) instead. CLM land surface data are applied as land use and land cover (LULC) 145 inputs. Details regarding these data can be found in Niu et al. (2017). More information on the governing equations of PAWS can be found in Shen and Phanikumar (2010) and Niu et al. (2014).

Hierarchical sensitivity analysis method
The essential concept of the hierarchical sensitivity analysis method involves categorizing and quantifying different complex uncertainties of certain model systems while considering their dependence relationships. Different uncertainty sources (or 150 uncertain inputs) are placed at different layers of a hierarchical uncertainty framework, which is then integrated with the variance-based global sensitivity analysis method to form a new set of sensitivity indices to accurately quantify the importance of different uncertainty sources.
In this study, we implemented the hierarchical sensitivity analysis method to assess the sensitivity of a large-scale processbased hydrological model, PAWS, with uncertainty from climate scenarios, numerical models and parameters. Similar to the 155 work of Dai et al. (2017a), there are also three layers of uncertainties in the hierarchical uncertainty framework constructed for this research, including the upper layer of scenario uncertainty, the middle layer of model uncertainty, and the bottom layer of parameter uncertainty (Fig. 2). However, unlike the previous work, we further decomposed the layer of parametric uncertainty into three groups of uncertainties to investigate the importance of different uncertain parameters (Fig. 2). All of https://doi.org/10.5194/hess-2020-87 Preprint. Discussion started: 10 March 2020 c Author(s) 2020. CC BY 4.0 License. 6 these uncertainties were placed into the proper levels based on their dependence relationships. Climate scenario uncertainties 160 were placed on the first layer because the climate scenarios (differing in precipitation, radiation, temperature, humidity, and wind speed) are the driving forces of the model system, multiple models can be built under a single scenario, and the choice of models depends on the characteristics of the scenarios. The second layer includes the numerical model uncertainties, which represent different plausible mathematical conceptualizations describing the study site. In this study, they represent three models with different aquifer thicknesses. Under each climate scenario, these three models are likely to coexist. 165 Parametric uncertainty is included in the bottom layer because each model can contain a different set of parameters (Meyer et al., 2007). Note that in Fig. 2, the term PR represents the multiple parameter sets under a certain numerical model. In the following content, we will introduce the revised mathematical formulas involved in the hierarchical sensitivity analysis method used in this study.
For the model: where  is the model output and   1 ,..., m X X X = is a group of uncertainty 170 inputs, using the law of total variance, the total variance of  can be decomposed as (Dai and Ye, 2015): where the first term of partial variance on the right-hand side is the within-i X partial variance and represents the variance contribution by i X and X i represents all the inputs except i X . The second term on the right-hand side represents the variance contributed by the model inputs excluding i X as well as the interactions of all the inputs. Based on the definition of 175 the first-order sensitivity index (Saltelli et al., 1995), the percentage of uncertainty contributed by input i X can be accurately quantified.
Using the hierarchical framework in Fig. 2, the variance-based sensitivity analysis method enables the decomposition of the total variance into individual contributors as follows: where CS represents the set of alternative climate scenarios and NM represents the multiple plausible numerical models, PR, is the uncertain parameter. The term ~CS represents the uncertain inputs excluding climate scenarios, which are NM and PR in this study. The term NM, PR|CS represents the change in the combination of the model and parameters under one fixed climate scenario. The first term of partial variance on the right-hand side of this equation represents the variance caused by multiple climate scenarios. The second term on the right-hand side is the partial variance caused by other uncertain inputs 185 and can be further decomposed as: https://doi.org/10.5194/hess-2020-87 Preprint. Discussion started: 10 March 2020 c Author(s) 2020. CC BY 4.0 License.
where the first partial variance term on the right-hand side of this equation represents the uncertainty contributed by multiple plausible models. The subscripts NM|CS and PR|NM, CS refer to the change in the models under one climate scenario and the change in the parameters under one model and climate scenario, respectively. The second term represents the within-190 model partial variance caused by the uncertain parameters. By substituting Eq. (3) back into Eq. (2), we can obtain: The three terms on the right-hand side of Eq. (4) represent the partial variances contributed by the parameters, models and climate scenarios, respectively. The equation indicates that the total variance can be decomposed into the variances contributed by the alternative climate scenarios, CS, plausible numerical models, NM, and uncertain parameters, PR. Then,195 we can define the new sensitivity indices for PR, NM and CS following the first-order sensitivity index definition (Dai et al., 2017a): In this study, the partial variance of the parameters was further decomposed to assess the subdivided parametric sensitivities contributed by the vadose zone parameters (PRVDZ), groundwater parameters (PRGW) and overland flow parameter (PROVN).
Using the total variance decomposition method, the partial variance of the parameters can be further decomposed as: https://doi.org/10.5194/hess-2020-87 Preprint. Discussion started: 10 March 2020 c Author(s) 2020. CC BY 4.0 License.
The subscript PR VDZ |NM, CS represents the change in vadose zone parameters under a fixed model and climate scenarios. 205 The subscript PR~V DZ |PR VDZ , NM, CS refers to other uncertain parameter inputs excluding vadose zone parameters, which are groundwater parameters and the overland flow parameter. The first term of Eq. (8) on the right-hand side is the partial variance of PRVDZ, and the second term represents the partial variance of the other parameters, which are groundwater parameters and the overland flow parameter. Note that Eq. (8) is decomposed at the level of vadose zone parameters; when we decompose the partial variance of parameters at the level of groundwater parameters or the overland flow parameter, the 210 partial variance of the parameters can be expressed as Eq. (9) and Eq. (10): The first term of Eq. (9) and (10)

Sensitivity index estimation using the Latin hypercube sampling and binning method 220
In this study, the parameters were sampled with a sample size of 600 within the feasible range via Latin hypercube sampling (LHS) (Kanso et al., 2006;Zhang and Pinder, 2003). LHS is a type of constrained Monte Carlo sampling that can accurately reflect the function distribution of the input data. Compared with Monte Carlo sampling, LHS greatly reduces the demand for sample size and computation time and is widely used in modeling simulation and optimization calculations. For the , by using LHS, the range of i X , 1, 2,..., ik = can be divided into n non-225 overlapping intervals with equal probabilities. The n values obtained from 1 X are randomly paired with n values obtained from 2 X ; these n paired values are then combined with those n values from 3 X . We repeat this process until the new nk  matrix, X , is developed. The sample matrix X can be used to calculate the sensitivity index for the model output. More details regarding LHS are described in previous studies (McKay, 1979;Owen, 1998;Helton, 2003).
Using the variance definition, the partial variance of V(PR) can be expressed as: 230 In this study, there are l=6 alternative climate scenarios, and under each climate scenario, there are k=3 plausible models, and we have n=600 sampled parameter sets for each model and climate scenario. After applying the formula of expectation and the LHS method, the terms of ( ) can be expressed as: where n and j represent the parameter LHS sample number and index, respectively, . The values of the weights for alternative models or climate scenarios could be 240 selected using prior knowledge or objective criteria, e.g., posterior probabilities of the Bayesian theorem (Neumann, 2012;Schoniger et al., 2014).
To calculate the subdivided parametric sensitivity indices, i.e., the sensitivity indices for vadose zone parameters, groundwater parameters, and the overland flow parameter, a binning method was implemented in this study. The main concept of the binning method is to approximate the mean value for a single parameter by substituting the mean value for 245 this parameter in a bin. Taking the sensitivity index of vadose zone parameters as an example, since (11) can be expressed as: Using LHS and the binning method, the number of realizations is reduced to the size of the parameter sets obtained from the LHS method. Thus, the computation cost for estimating the subdivided parametric indices can be highly reduced. In this case, the sample size of PRVDZ is 600, and the number of possible combinations of PRGW and PROVN is 600×600 2 ⁄ =180,000; then, the number of realizations of regular Monte Carlo simulations is 600×180,000=1.08×10 8 . Using LHS and the binning method, the number of realizations is reduced to 180,000. Dai et al. (2017b) confirmed a similar accuracy of 36,000,000 265 Monte Carlo realization results with 16,000 realizations when applying only the binning method for a synthetic example.
The combination of the LHS method with the binning method makes it possible to further reduce the computational cost for such a large-scale, complex hydrologic model.

Uncertainty sources and the generation of uncertain inputs
In this study, climate scenarios, model structures, and parameters are treated as random uncertain inputs or uncertainty 270 sources to assess the sensitivity for two model outputs of interest: evapotranspiration (ET) and groundwater contribution to streamflow (QG). For the climate scenarios, we generated six typical and alternative scenarios based on NASA's Tropical Measuring Mission (TRMM) data (http://trmm.gsfc.nasa.gov/) and the default CLM CRU-NCEP (CRUNCEP) dataset (Piao et al., 2012) from 1998 to 2013. The precipitation data were obtained from the TRMM while the temperature, solar radiation, humidity and wind speed data are based on the CRUNCEP because the model fails to capture the peak stream discharges 275 using the CRUNCEP rainfall data (Niu et al., 2017). We first divided the full climate dataset into dry and wet seasons according to the precipitation values (six months for each season). Then, we sorted the different seasons according to their total precipitation values during the whole season. Next, we divided these wet and dry seasons into three different groups representing six climate scenarios from wet to dry (Fig. 3). The mean and standard deviation of the values of the different climate variables (e.g., precipitation, maximum temperature) for each group were calculated using the daily data (Table 1). 280 Finally, we generated random daily weather data for each climate scenario based on these mean and standard deviation data using a normal distribution. The mean and standard deviation for each climate scenario's daily data are listed in Table 1, and https://doi.org/10.5194/hess-2020-87 Preprint. Discussion started: 10 March 2020 c Author(s) 2020. CC BY 4.0 License. Fig. 3 displays a box plot of the precipitation data for the six climate scenarios (CS1, CS2, CS3, CS4, CS5, CS6). In this study, we assumed that the different scenarios have equal probability. However, we changed the probabilities for CS1 (the wettest climate scenario) and CS6 (the driest climate scenario) in Section 3.5 to investigate the influences of extreme climate 285 scenarios on the sensitivity analysis results.
According to Pelletier et al. (2016) assumed to have equal weights (probabilities), but changes in these probabilities were also investigated in Section 3.5.
We adjusted the values of six model parameters and classified them into three groups to explore the parametric sensitivities.  Table 2. To ensure that the sensitivity analysis was uncomplicated and computationally tractable, each grid was given the same value for a specific parameter.

Model predictions 305
The total number of PAWS+CLM simulations considering all possible combinations of the three uncertain factors is 6×3×600=10,800, which represents six climate scenarios, three model conceptualizations of aquifer thickness, and 600 sampled parameter sets. In this study, we first used equal weights for the climate scenarios and numerical models, i.e., P(CS l )= 1 6 and P(NM k |CS l )= 1 3 . In Section 3.5, we investigated the influences of different prior probabilities for the climate scenarios and numerical models. The simulation time for all the simulations was six months (180 days, 4320 hours), which is 310 the length of the dry or wet season in the central Amazon region. The results given by PAWS were represented in two forms: (1) space-averaged output values over the whole grid at each time step; (2) time-averaged output values over the whole simulated period for each grid. In this study, the time step is one hour. Figure 4 depicts the spatially averaged model https://doi.org/10.5194/hess-2020-87 Preprint. Discussion started: 10 March 2020 c Author(s) 2020. CC BY 4.0 License.

13
predictions at the simulation time of 2868 hours (wall-clock time: 12:00) for the two outputs of interest, ET and QG, with different inputs of scenarios, models, and parameter sets. The results demonstrate the large variabilities or uncertainties 315 among the predictions of the different model realizations, and these uncertainties were contributed by all three sources of uncertainty: scenarios, models, and parameters. Further accurate sensitivity analysis is capable of identifying the most important sources of uncertainty.

Sensitivity indices for evapotranspiration
We first calculated the sensitivity indices for the spatially averaged ET over the whole watershed at all time steps using Eqs. 320  Figure 5 indicates that the sensitivity to various factors is strongly temporally dependent. However, it should be noted that at the wall-clock time of 12:00-13:00, the climate scenarios are always the most important factors affecting the sensitivity of ET because ET is directly influenced by solar radiation values, and the radiation forcing used in this study reaches its maximum value at approximately 12:00. At the wall-clock time of 24:00-1:00, the sensitivity indices for the parameters (S PR ) show absolute dominance since precipitation and radiation forcing all decrease to 330 zero, leading to the differences in rainfall and radiation among the climate scenarios greatly decreasing.
Six time points (simulation times = 1428 hour, 1440 hour, 2868 hour, 2880 hour, 4308 hour, and 4320 hour) were chosen as examples to show the sensitivity indices (Fig. 6). Simulation times of 1428 hour, 2868 hour, and 4320 hour belonged to different days but all corresponded to 12:00 local time. At these time points, the climate scenario uncertainty (S CS ) is the most important contributor to the total ET prediction uncertainty, accounting for 54-77% of the total uncertainty, and 335 parameters (S PR ) contribute the second most to uncertainty. However, at different time points (1440 hour, 2880 hour, and 4320 hour, corresponding to 24:00 local time), the parameters are the dominant uncertainty contributor, with S PR ranging from 89 to 92%. We also calculated the sensitivity indices for every grid cell within the model domain using the timeaveraged ET predictions over all simulation periods (4320 hours). Figure 7 shows the spatial variability of the sensitivity indices for the temporal mean ET predictions. The maps demonstrate that for most grids, parameters are the most important 340 uncertainty contributor to ET predictions (S PR >0.50), and climate scenarios are the second most important contributor to uncertainty. However, for the stream grid cells, the importance of model uncertainty increases. The parameters and aquifer thickness are both important for the ET predictions in river grid cells. This sensitivity occurs because the aquifer thickness along the streams will affect the exchange between the groundwater cells and the river cells (the relevant process is shown in Appendix A). 345 https://doi.org/10.5194/hess-2020-87 Preprint. Discussion started: 10 March 2020 c Author(s) 2020. CC BY 4.0 License.

Sensitivity indices for groundwater contribution to streamflow
The same sensitivity analysis procedures were also conducted for the model predictions of QG. Groundwater is an important component of hydrologic processes in the central Amazon region; therefore, it is essential to conduct sensitivity analysis for QG. Figure 8(a) shows that regardless of the time step, parameters are always the dominant contributor to the total QG prediction uncertainty because soil parameters strongly affect the soil water redistribution process, including the infiltration 350 into groundwater. The sensitivity indices of the models ( S NM ) and climate scenarios ( S CS ) reach peak values at approximately 1:00 (Figs. 8(b)-(g)). This may be because the exchange between groundwater and river flow always occurs hours later than the rainfall process, and the value of the exchange process always reaches its peak at night, at approximately 1:00. Furthermore, the thickness of an aquifer will greatly influence the water redistribution process in the aquifer. It should be noted that, in Fig. 8(a), the values of the sensitivity indices for the climate scenarios are cumulative. In the beginning, the 355 importance of the climate scenarios is negligible, but over time, the importance of the climate scenarios increases. This pattern can be explained by the fact that the climate scenarios have no immediate direct influence on groundwater flow NM S , but over time, the amount of groundwater in aquifers will be affected, and this influence is long-term and cumulative.
Because groundwater exchange with stream flow occurs only at grid cells along the streams, the sensitivity indices only have valid values in those stream grid cells (Fig. 9). Our results indicate that for most stream grid cells, the parameters are the 360 most important contributor to the total uncertainty of time-averaged QG predictions. The second most important factor is aquifer thickness.

Sensitivity indices for subdivided parameters
Based on the sensitivity analysis for ET and QG, the results show that parameters are the most important uncertain input for both the spatially averaged and time-averaged uncertainties. In this study, we used Eqs. (11)-(13) to further calculate the 365 subdivided parametric sensitivity indices, which can provide more detailed sensitivity analysis for model simulation.
Through this investigation, the parametric sensitivity was subdivided into three groups: (1) the sensitivity for vadose zone parameters (PRVDZ), (2) the sensitivity for groundwater parameters (PRGW), and (3) the sensitivity for the overland flow parameter (PROVN). Using the binning method, we calculated the spatially averaged and time-averaged subdivided parametric sensitivity indices for ET and QG. We plotted frequency histograms of the subdivided parametric sensitivity 370 indices over 4320 hours in Fig. 10. Figure 10(a) depicts the results for ET. The value of S PR VDZ is concentrated in the range of 0.1-0.9, and S PR GW is concentrated in the range of 0.003-0.032. The value of S PR OVN is so small that the influence of the overland flow parameter can be ignored.
This indicates that vadose zone parameters (PRVDZ) dominate the total parametric uncertainties for ET. Figure 10(b) shows the frequency histogram of spatially averaged subdivided parametric sensitivity results for QG. S PR VDZ is still concentrated in the larger number range (0.2-0.53), and the value of S PR GW changes from 0.04 to 0.3. The number of S PR OVN is also the lowest, indicating that the overland flow parameter has little effect on QG.
We plotted the time-averaged subdivided parametric sensitivity indices for ET in Fig. 11(a) and for QG in Fig. 11(b).
Considering ET as our output, for most grids, the vadose zone parameters are the most important contributor to parametric uncertainties. Compared with that on other grids, the influence of groundwater parameters on the river grids is more 380 significant ( Fig. 11(a)). In terms of QG, for most grids, the vadose zone parameters dominate the parametric sensitivities ( Fig.   11(b)). For both time-averaged ET and QG, the impact of the overland flow parameter can be ignored.

Effects of prior weights on sensitivity indices
In this section, we changed the prior weights of the climate scenarios and numerical models to investigate their influences on the space-averaged sensitivity indices. Because the number of space-averaged results for ET and QG is too large to be well 385 exhibited, we chose one time step (4308 hour, 12:00 wall-clock time) to show the trend. We randomly changed the values of the prior weights for NM1 (the thickness of the unconfined aquifer is 50 m, and that of the confined aquifer is 250 m), CS1 (the wettest climate scenario), and CS6 (the driest climate scenario) to between 0 and 1. The values of these factors' prior weights were randomly changed. Figure 12(a) indicates that when we consider ET as our output, with the increase in the prior weight of NM1, the uncertainty of the climate scenarios will decrease to 50%, while the uncertainty of the parameters 390 will increase to 50%. Both parameters and climate scenarios have quite important effects on ET. Different from the results for ET, with the increase in the prior weights of NM1, the sensitivity index of the numerical models for QG decreases to 0 (because only one model exists under this condition), and the scenario uncertainty changes only slightly. Moreover, the uncertainty of parameters always dominates the total uncertainty for QG ( Fig. 12(b)) regardless of the prior weight value. In general, the different prior weight values for the numerical models only slightly change the sensitivity analysis results. 395 Figures 12(c)-(f) exhibit the influences of prior weights for the wettest and the driest climate scenarios on ET. These figures first demonstrate that changing the values of the prior weights of CS1 and CS6 have larger impacts on ET predictions than on QG predictions. This pattern coincides with the fact that the parameter uncertainty dominates the total predictive uncertainty of QG and that the scenario uncertainty is relatively small. Therefore, the selection of prior weight values for the scenarios does not have a significant effect on the sensitivity analysis results for the QG predictions, and the parameter sensitivity index 400 is always the largest (Fig. 12(d) and (f)). For the sensitivity analysis results pertaining to ET predictions, changing the values of the weights for CS1 and CS6 has different effects. The sensitivity index values of the climate scenarios for ET predictions monotonically decrease while the importance of parameters continues to increase as the prior weight of CS1 rises. However, the value of S CS for ET predictions first increases and then dramatically decreases after the prior weight of CS6 approaches 80%, and S PR shows the opposite trend ( Fig. 12(e)). The different effects on ET predictions of CS1 and CS6 might be caused 405 by their distinct temperatures and radiation intensities.

Conclusions
In this study, we implemented an advanced hierarchical sensitivity analysis method and a complex large-scale process-based hydrological model (PAWS) to identify important uncertain inputs for ET and QG predictions in an Amazonian watershed.
This sensitivity analysis method is capable of providing accurate measurements of the importance of uncertain inputs 410 through variance decomposition, and it can also categorize and combine different uncertain inputs considering their dependence relationships to decrease the high dimensionality induced by a complex and large-scale problem. Three groups of uncertainty sources or uncertain inputs were considered in this study, including six climate scenarios, three plausible aquifer models, and six uncertain parameters (i.e., soil saturated conductivity, Van Genuchten α and N, unconfined aquifer conductivity, confined aquifer conductivity, and the length of the flow path for runoff contribution to the overland flow 415 domain). A new set of subdivided parametric sensitivity indices was defined for three groups of parametric uncertainty sources (i.e., vadose zone, groundwater, and overland flow parameters). When researchers use this hierarchical sensitivity analysis method, and there are some processes at a certain layer that can be further separated, this study provides a feasible way to investigate more detailed information for the grouped uncertainties. The implementation of the Latin hypercube sampling method and the binning method reduced the high computational cost. 420 The sensitivity analysis results in this study demonstrate that parameter uncertainty is important in both time-averaged and space-averaged predictions regardless of whether the output is evapotranspiration or groundwater contribution to streamflow.
Among all the parameter uncertainties, the vadose zone parameters are the most important, and the parameter of overland flow is negligible. The climate scenarios are also important uncertainties in evapotranspiration predictions, especially at the wall-clock time of 12:00 noon. Along the river grid cells, the thickness of the aquifer has a non-ignorable influence on both 425 evapotranspiration and groundwater contribution to streamflow. On the basis of the results of this study, we suggest that when modelers use sophisticated hydrological simulators such as PAWS, they should pay attention to the weather variable values at approximately 12:00 noon (always the daily peak values), investigate the thickness of groundwater aquifers near rivers and adjust the parameters of the vadose zone.
This study represents a pilot example of comprehensive global sensitivity analysis considering all uncertainty sources in a 430 large-scale hydrological model. The sensitivity analysis results can provide key information on uncertainty sources for modelers and greatly improve the model calibration and uncertainty analysis processes. By categorizing multiple uncertainties into processes and placing them into a proper layer in a hierarchical framework, this advanced hierarchical sensitivity analysis method can largely reduce the computational cost associated with complex, large-scale hydrological models. Its combination with Latin hypercube sampling and the binning method can further decrease the computational cost. 435 The proposed method is mathematically rigorous and general and can be applied to extensive, large-scale hydrological or environmental models with more or different sources of uncertainty.
where h represents the soil water pressure head, z is the elevation (positive upward), K(h) represents the soil unsaturated conductivity and W(h) is the source or sink term, including the influence of evaporation, root extraction and lateral flow. The differential water capacity can be described as C(h)= ∂θ ∂h ⁄ , where h is the soil pressure head and θ is the water content. The 450 pressure head, h, is related to the unsaturated hydraulic conductivity, K(h). According to the Mualem-van Genuchten formula (Mualem, 1976;van Genuchten,1980), the soil saturated hydraulic conductivity, Ks, Van Genuchten α and N will influence the unsaturated conductivity, K(h): where S is the relative saturation, θs is the saturated water content, θr is the residual water content, N is related to the poresize distribution, α indicates the reciprocal of air suction and λ is a parameter measuring pore connectivity.
The aquifers in PAWS are depicted as a series of 2-D layers . In each layer, the 2-D groundwater equation where S is the storability; T is the transmissivity of the aquifer; T=Kb, where K is the aquifer conductivity and b is the saturated thickness of the aquifer; H is the aquifer hydraulic head; R is recharge or discharge; W is the source and sink term; and Dp is percolation into deeper aquifers.
PAWS applies one-dimensional diffusive wave equations to portray the channel flow model (Shen and Phanikumar, 2010;Shen et al., 2014). After calculating the channel flow, the exchange between groundwater and channel flow (QG) is 465 immediately computed. The calculation of QG is based on the leakance concept (Shen and Phanikumar, 2010 where w is the wetted perimeter. If the river width is greater than 10 m, w can be approximated as the river width. PAWS retains its own flow scheme, but the surface processes use the CLM 4.0 model. This enables the simulation of detailed surface processes such as surface heat flux, water vapor flux, surface radiation balance, crop growth, and plant 475 photosynthesis. The calculation of ET demand is performed in the CLM model based on the climate data, and then ET demand will be transferred to PAWS as a source term for the vadose zone. More details about the calculation of ET (both evaporation and transpiration information can be found in the technical note of CLM 4.0, http://www.cesm.ucar.edu/models/cesm1.1/clm/CLM4_Tech_Note.pdf). The coupling with the CLM makes PAWS a more comprehensive and robust surface-subsurface hydrological model. 480      https://doi.org/10.5194/hess-2020-87 Preprint. Discussion started: 10 March 2020 c Author(s) 2020. CC BY 4.0 License.