Robust global sensitivity analysis of a river management model to assess nonlinear and interaction effects

The simulation of routing and distribution of water through a regulated river system with a river management model will quickly result in complex and nonlinear model behaviour. A robust sensitivity analysis increases the transparency of the model and provides both the modeller and the system manager with a better understanding and insight on how the model simulates reality and management operations. In this study, a robust, density-based sensitivity analysis, developed byPlischke et al.(2013), is applied to an eWater Source river management model. This sensitivity analysis methodology is extended to not only account for main effects but also for interaction effects. The combination of sensitivity indices and scatter plots enables the identification of major linear effects as well as subtle minor and nonlinear effects. The case study is an idealized river management model representing typical conditions of the southern Murray– Darling Basin in Australia for which the sensitivity of a variety of model outcomes to variations in the driving forces, inflow to the system, rainfall and potential evapotranspiration, is examined. The model outcomes are most sensitive to the inflow to the system, but the sensitivity analysis identified minor effects of potential evapotranspiration and nonlinear interaction effects between inflow and potential evapotranspiration.


Introduction
Water managers rely heavily on models to predict future water availability, optimize water use and evaluate water man-30 agement strategies in order to find a balance between environmental, social and economic demands on the system.It is therefore crucial to be aware of the ability of a model to capture the dynamics of the hydrological cycle relevant to the water management question.In recent decades, address-35 ing this issue has been the focus of much research in hydrological model calibration and predictive uncertainty analysis (Gupta et al., 2012).For a modeler, to arrive at a 'well'-calibrated model or to produce sensible and robust prediction intervals, it is essen-40 tial to have a thorough understanding of how the hydrological system works and how this system is represented in the model; how a variation in parameters, boundary conditions or driving forces will affect the prediction of interest.The knowledge gained from such sensitivity analysis is 45 not only of relevance during model development, it provides added value to the model as it can focus management and monitoring to those aspects of the system and model that are most important to the management of water resources (Saltelli et al., 2008).Additionally, discussing model sensi-50 tivities with stakeholders will remove the notion of the model being a 'black box' and can provide stakeholders with a better appreciation of the accuracy of the model, which has proven to be a key aspect of adoption of model results in management (Patt, 2009;Bark et al., 2013).

55
River management models such as eWater Source (Welsh et al., 2013) are increasingly used, especially in Australia, in the development of basin-wide water allocation plans.As these plans directly affect the livelihood of people and the health of ecosystems, it is essential that the models underpinning these plans have wide support and are robust.It is therefore essential that practitioners have a set of tools for sensitivity analysis available, tailored to the needs of water allocation modelling.The most straight forward sensitivity analysis technique is One-At-a-Time (OAT) sensitivity analysis in which one model aspect is changed while the others are fixed.The sensitivity of the model output to variation of the tested parameter is proportional to the gradient of the response surface.This is formalized in gradient-based calibration routines, such as Levenberg-Marquardt optimization.Examples of such OAT sensitivity analysis are Doherty and Hunt (2009), Foglia et al. (2009), Castaings et al. (2009) and Peeters et al. (2011).This methodology is attractive as it requires a very limited number of model runs, about 2 to 3 model runs per parameter evaluated, and, as long as the model behaves linearly, parameter interaction effects can be explored (Hill and Tiedeman, 2007).Saltelli and Annoni (2010) highlight that OAT sensitivity analysis only provides reliable and robust results if it can be shown that the model behavior is linear.This condition is seldom satisfied for hydrological models or even known before a sensitivity analysis.The Elementary Effects method (Campolongo et al., 2007) is more robust against non-linearity in the model behavior, whilst still being frugal in the number of model runs.Global sensitivity analysis techniques however do not require the model behaviour to be linear (Saltelli et al., 2008).The most straightforward global sensitivity analysis is either random or density based sampling of parameter space and visualizing scatter plots of the parameter value against the prediction of interest (Wagener and Kollat, 2007;Peeters et al., 2013).Variance based methods, such as Sobol' sensitivity analysis (Saltelli and Annoni, 2010;Nossent et al., 2011), use a scheme of structured resampling of a random base sampling to decompose the variance of the metric of interest into the main effects of a parameter and interaction effects of other parameters.The main drawback of variance based methods is that it assumes that the entire effect of a parameter can be summarized by the variance (Borgonovo, 2007;Borgonovo et al., 2011).Variance based sensitivity indices will therefore be less reliable if the response to a parameter has a skewed or multi-modal distribution.Density-based sensitivity analysis techniques attempt to account for this by incorporating the entire distribution of the response of a prediction of interest in the metric in a way that does not require any assumptions on the shape of the distribution.The methodology suggested by Plischke et al. (2013) implements such a density-based sensitivity analysis technique which is independent of the parameter sampling scheme.This has the added benefit that as no model runs need to be devoted to the resampling of a base sampling, more computing resources can be directed to exploration of parameter space.The goal of this study is to apply a density-based sensitivity analysis in a river management modelling context to assess its capability to identify and quantify non-linear effects and 115 to extend the methodology to account for interaction effects.An idealised, hypothetical river management model implemented in the eWater Source platform (Welsh et al., 2013) serves as testing platform to assess the ability of the sensitivity analysis methodology to quantify the influence of a 120 small number of forcing variables upon a variety of model outcomes.The next section presents the theoretical background and numerical implementation of the Plischke et al. (2013) global sensitivity analysis method.The river management model is 125 briefly introduced before presenting the results of the sensitivity analysis and summarizing the findings in the discussion and conclusion sections.

Methods
The sensitivity analysis introduced in Plischke et al. (2013) 130 provides a robust, global density-based sensitivity analysis, independent of sampling strategy.This section provides a short summary of this methodology, for a detailed overview the interested reader is referred to Plischke et al. (2013).Consider X and Y the set of variables that comprise the in-135 put and output respectively of a river system model.Fixing X to a single realisation, the parameter combination x, results in a conditional cumulative distribution of Y equal to F Y |X=x (y) and an equivalent density function f Y |X=x (y).The importance of fixing X to x can be quantified by the sep-140 aration between the unconditional F Y (y) and the conditional F Y |X=x (y) or, similarly, the separation between f Y (y) and f Y |X=x (y).Using the L1-norm, the separation between the two density functions can be written as: The importance of factor X on outcome Y can then be defined as: The sensitivity index δ(X, Y ) varies between 0 and 1 and it can be shown that this index is zero when X and Y are completely independent (Plischke et al., 2013).
To compute δ(X, Y ) the integrals in eq. 2 need to be ap-155 proximated numerically.This can be achieved by taking n samples of the parameter space X and compute the corresponding values for Y .The method does not impose any restrictions on the sampling strategy of the parameter space.dom sampling, quasi-random sampling (e.g.Latin Hypercube Sampling or Sobol' sequences) or Markov Chain Monte Carlo simulation.The resulting dataset is partitioned into M classes C m with m = 1, ..., M .For each class C m , the density function can be approximated with a kernel smoothing function with kernel K(.) and bandwidth α (Devroye and Gyorfi, 1985): where n m is the number of samples in class C m and α m the corresponding bandwidth for the kernel smoothing function.The next step is to approximate the L1 norm between the two distributions for each class.Using a predefined number of quadrature points {ỹ j , j = 1, ..., l}, the separation can be computed as: 180 The sensitivity index δ can then be approximated by: To avoid bias in the sensitivity index and to assess the robustness of the sensitivity index estimate, it is recommended to perform a bootstrap of the sensitivity index (Efron, 1977) and to adjust δ with the mean of the bootstrap δ * : δ provides the sensitivity index of the main effect of a variable.Plischke et al. (2013) however does not provide a method to explore second order effects, i.e. the interaction 190 between two variables.To estimate second order effects between variables X 1 and X 2 , the samples are subdivided into n groups of equal intervals for X 1 .The sensitivity index δ for X 2 , δX2 , is computed for each interval.If there is no interaction effect between X 1 and X 2 , then δX2 will not vary with the level of X 1 .To quantify this, the variance of δX2 is computed over all n levels of X 1 .Small variances indicate small interaction effects and vice versa.

Model Description and Setup
The case study is a hypothetical river system model (Fig. 200 1), based on a simplified version of the Murrumbidgee River Model in New South Wales, Australia (Dutta et al., 2012;Podger et al., 2014).Using the full version of the Murrumbidgee River Model was not warranted, not only because of the complexity of the system and the manage-205 ment rules, but, more importantly, because of legal issues with regards to model licensing and confidentiality.The idealised, hypothetical model retains most of the relevant complexity practitioners encounter when creating water allocation models, which is more than sufficient to illustrate 210 the sensitivity analysis methodology.In the model, water is routed from a storage reservoir through three river reaches.Routing starts in reach 1 at the storage reservoir with hydropower generators that receive water from a single tributary inflow.In Reach 1, water is 215 taken from the system for town water supply and irrigation and water is received from unregulated rain-fed tributaries.From the Upper Gauge at the end of Reach 1, water is routed through Reach 2. In this reach, interaction with groundwater is taken into account by an exchange flux.As in reach 1, 220 water is received from unregulated, rain-fed tributaries and water is taken out for irrigation and town water supply.In addition to these offtakes, water is diverted into an off-river wetland system.Reach 3 starts at the middle gauge and is similar to reach 2. It also has offtake for town water supply, 225 irrigation and off-river wetlands and recieves inflow from rainfed tributaries.Groundwater-surface water interaction is not taken into account in this reach.Each reach has a term representing unaccounted losses.The loss relationships are taken from the more complex model.The total travel time 230 from headwater to end-of-system is 18 days (3 days reach 1, 6 days reach 2 and 9 days reach 3).These values, together with the other parameters influencing routing of water are also taken and aggregated from the more complex model.types.There are environmental demands for the wetlands in reach 2 and 3, which are designed to establish and maintain favorable habitat conditions for indigenous fauna and flora (Janssen, 2012).Two aspects of water management are considered: 347m 3 /s 250 order constraint on storage releases, i.e. the maximum flow that can be requested by water users in the system of the

Results
In the sensitivity analysis, the three main forcing vari-295 ables are considered; the system inflow (Inf low), the precipitation (Rain) and the potential evapotranspiration (P ET ).The latter two affect the inflow into the reaches and the irrigation demand.Inspired by the work of Leblanc et al. (2012), the forcing variables are changed through a 300 multiplier to the corresponding input time series with the range of the multiplier for each variable is to be between 0.5 and 1.5.This range encompasses both historical variation in hydrological input and output, as well as the expected change under various climate change models and scenario's.While 305 elaborate schemes are available to perturb hydrological time series, this is not warranted in this study as the focus is on metrics that integrate the entire flow time series.As such the emphasis of this research is on changes in total flow in or out the model, rather than in changes of the timing of flow.

310
Using Sobol' sequences (Sobol, 1976), 100, 000 quasirandom samples of the three input variables are generated.For each of these samples a range of output time series is calculated (Pickett et al., 2013).Table 1 lists the names of the output series and a short description.
The choice of this metric is motivated by the fact that, since the case study is an idealised, hypothetical model, it is not possible to directly compare the results with observations.In addition to this, and more importantly, the variety of model outcomes examined in this study are more than likely to be affected by different aspects of the hydrograph.Similar to choosing an objective function in traditional calibration or a likelihood function in uncertainty analysis, such metric needs to be tailored to be able to capture the relevant aspects of the hydrograph.Choosing an ill-suited metric can have huge consequences for the sensitivity analysis, calibration or uncertainty analysis, as pointed out in Montanari and Koutsoyiannis (2012) and Nearing (2014).The metric presented in Eq. 7 is designed to provide an as general and robust as possible measure of the difference between two time series as not to bias the interpretation of the sensitivity analysis.

Main Effects 345
Fig. 2 shows the scatter plots of sensitivity metric M for all combinations of forcing data and output variables.It is clear that the dominant influencing driving variable is Inf low as a strong response is noticeable for variations in this driving variable for all output variables with the exception 350 of HydroP ower.The effects of Rain and P ET are less pronounced.A very striking feature are the many nonlinearities in the response surface of the hypothetical model.This is mostly due to a number of threshold values used in the management rules of the river management system.For 355 instance, generation of hydro-power is only possible when the storage level in the dam exceeds a predefined threshold related to the height of the water intake point for the turbines.Inf low and the same is true for the output variables related to monetary value ($AlgalBloom, $Stor and $T otalAg).HydroP ower is least influenced by Inf low, which, from Fig. 2, is clearly related to the threshold-induced non-linear behavior.

370
The methodology is also able to quantify the often small and non-linear effects of the other forcing variables.This is especially noticeable for P ET .There is a clear but highly nonlinear effect of P ET on $Stor, which is reflected in a higher δ.The output variable HydroP ower has a bimodal distri-375 bution where the majority of simulations have an M close to zero.Nevertheless, the global sensitivity method is able to distinguish and quantify the subtle trends in the non-zero values for the different input variables.

Interaction Effects
The previous section established the importance of Inf low as the main driving variable.It is however from both a management and modeling perspective interesting to have an understanding of how the interaction between variables affects the model outcome.
Fig. 4 shows plots with the factor values on the x and y-axis, with a color scale to visualise M for the three combinations of interaction of the driving forces (Inf low-Rain, Inf low-P ET and Rain-P ET ) for all 8 model outcomes.
The first column shows that the effect of Inf low on most of the model outputs does not vary with the value of Rain.
There is however a clear interaction between Inf low and P ET for most of the model outputs; while the Inf low response is the dominant feature in the plots, the shape of this response depends on the value of P ET .HydroP ower is a noted exception as it displays very little structure in the scatterplots.This is because hydropower is generated by release of water from the reservoir in function of the demand and the water level in the reservoir.These management rules create a buffer to immediate impact from rainfall and inflow and also result in non-linear, threshold related behaviour.Very little structure is noticeable in the third column of Fig. 4, which shows the interaction between Rain and P ET , reflects the limited influence both driving forces have as a main effect.To quantify the interaction effect for each interaction combination in Fig. 4, the variance of the δ of the variable on the y-axis is computed for 100 equal intervals of the variable on the x-axis.By using Sobol' sequences to generate the 100, 000 samples of the parameter space, each equal interval of the x-axis variable has approximately 1, 000 samples to compute the δ.The most dominant interaction effects are between Inf low and P ET for $T otalAg and U pperF low, followed by $AlgalBloom, $Stor and M iddleF low.

435
The sensitivity analysis of the hypothetical river management model highlights inflow as a crucial variable of the model and how this affects the economic, environmental and sociological functions of the river.This emphasizes the importance of an accurate characterization of the flow rates of upstream 440 areas when modeling flow routing in regulated systems comparable to the case study, i.e. the regulated river systems of the Murray Darling Basin in Australia.An accurate characterization of flow rates not only entails maintaining a dense river gauge network, it also means adequately describing the 445 measurement uncertainty in the flow rates, not in the least the uncertainty introduced by the rating curve that describes the stage-discharge relationship (Tomkins, 2012).The work of Hughes et al. (2014) illustrates this as they identify the inflow from ungauged catchment as crucial in the calibration 450 of river management models.Direct precipitation in the storage, wetlands and irrigation areas has a very minor influence on the model outcomes.This is mostly due to the small volume of rainfall (0.633km 3 /yr) compared to the inflow volume (4.4km 3 /yr) and the corre-455 lation between the inflow volume and rainfall.Any effect of rainfall will therefore be dwarfed by the effect of inflow to the system.The interaction effect of Inf low and P ET is mostly due to the feedback mechanism as irrigation requirements increase with increasing potential evapotranspiration.

460
Such parameter interaction is well-known in other areas of hydrological modelling, such as in rainfall-runoff modelling (Gallagher and Doherty, 2007;Zhang et al., 2013;Peeters et al., 2013) and in groundwater modelling (Doherty and Hunt, 2009), although it has not received much attention 465 in river system modelling.Letcher et al. (2007) discuss the importance of interacting effects in water allocation models, without however providing a rigorous quantitative framework to evaluate the effects.The sensitivity analysis in this study was limited to multiply-470 ing factors on three driving forces.It would be very insightful to include other model parameters in the sensitivity analysis, especially those controlling storage volumes and irrigation requirements.Along the same lines, including the parameters of the management rules, e.g.rules on allocations, in 475 the sensitivity analysis can yield additional understanding of the operational management of the river system, as shown by Micevski et al. (2011).

Conclusions
The density-based sensitivity analysis of Plischke et al.

485
The extended sensitivity analysis method presented in this paper provides a quantitative measure of sensitivity of main and interaction effects and, through a combination with qualitative visual inspection of scatter plots, proved to be able to identify not only major effects but also subtle interactions, 490 even in the presence of strong non-linearities.Due to the small dimensionality of the case study, it was possible to visualise all main effects and their interactions through scatter plots for all model outcomes.Although this will be challenging for higher dimensional problems, the visual inspection of scatter plots is an invaluable complement to the sensitivity indices.Understanding the dynamics of river system models is often not intuitive, especially in larger or basin-scale models (Johnston and Smakhtin, 2014).A robust and comprehensive sensitivity analysis is an invaluable step in model development to elucidate the often intricate interactions between driving forces, management rules and parameters.Increased understanding of the model will not only lead to improvements in calibration and prediction, it also has enormous potential in 505 establishing credibility and understanding of models.

Figure 1 .
Figure 1.a) Map showing the extent (indicated by pink shading) of the idealised river system model within the Murray-Darling Basin and b) schematic structure of the river management model

Fig. 3
Fig. 3 shows a barplot of the sensitivity indices δ for all

Figure 2 .
Figure 2. Scatter plots of M , the difference between kernel density estimates for each simulation and the kernel density estimate of the reference simulation for all forcing data and model output variables for the eWater Source hypothetical river management model.

Figure 3 .
Figure 3. Sensitivity indices, δ, for all forcing data and model output variables for the eWater Source hypothetical river management model.

Fig. 5
Fig.5illustrates this for the interaction effects of Inf low, Rain and P ET on $Stor.The sensitivity index values for Rain are low and hardly vary for different levels of Inf low, which is an indication of very limited interaction between Rain and Inf low, as confirmed by the scatterplot (Fig.4).The δ values for P ET do vary markedly with the level of Inf low.This sensitivity index reaches a minimum for Inf low values close to 1, while reaching peaks close to values of 0.75 and 1.1.This is reflected in the variance of the δ values which is 4.5x10 −4 for the Inf low-Rain couple and 3.5x10 −3 for Inf low-P ET .Fig.6shows the variance of the sensitivity indices for all interaction pairs for all model outcomes.The values for Hydropower are much higher than for the other model outcomes due to the nonlinear behaviour.They are omitted from Fig.6as they distorted the visualization.
480(2013) has been applied to a river management model rep-

Figure 4 .
Figure 4. Scatterplots of interaction of the driving forces.The intensity of the color scale is proportional to the model outcome value, where dark red colors indicate high values and light red colors indicate low values

Figure 5 .
Figure 5. Sensitivity index δ of the effect of Rain (blue) and P ET (red) on $Stor for 100 equal intervals of Inf low

Table 1 .
Output variables of the Source river system model Each of the output variables in Table1is a daily time series.The metric for the sensitivity for different forcing data ( M ) is the difference between the kernel density estimate of the daily times series of a randomly selected reference simulation ( fY ref (y)) and the kernel density estimate of the daily time series for the changed forcing data ( fY sim (y)): Saulnier, G.-M.: Sensitivity analysis and parameter estimation for distributed hydrological modeling: potential of variational methods, Hydrology and Earth System Sciences, 13, 503-517, http://www.hydrol-earth-syst-sci.net/13/503/2009/,2009.Crase, L. and Gillespie, R.: The impact of water quality and water level on the recreation values of Lake Hume, Australasian Journal of Environmental Management, 15, 21-29, 2008.DRET: Destination visitor survey: strategic regional research report -New South Wales, Victoria and South Australia.Impact of the drought on tourism in the Murray River Region., Tech.rep., Department of Resources, Energy and Tourism, Tourism Research Australia.57pp., 2010.