HydroTrend is a numerical model that creates synthetic river discharge and sediment load time series as a function of climate trends and basin morphology.
In this example, we'll perform an experiment on HydroTrend using the stoch_collocation method to evaluate how changing two input parameters:
starting_mean_annual_temperature (T) andtotal_annual_precipitation (P)affects the median values of one output parameter, long-term suspended sediment load at the river mouth (Qs), over a 10-year run.
Before we start, make sure that you've installed Dakota, HydroTrend, and this package on your computer, using the instructions in the README file.
Start by importing the Dakota class.
In [ ]:
from dakotathon import Dakota
Create a Dakota instance to perform a study of HydroTrend using the stochastic collocation method and input parameters characterized by normal distributions.
In [ ]:
d = Dakota(method='stoch_collocation', variables='normal_uncertain', plugin='hydrotrend')
Define the HydroTrend input variables (T and P) to be used in the study. We assume they're random variables with normal distributions, and we can obtain their default mean and standard deviation values from the HydroTrend parameters file (discussed below).
In [ ]:
d.variables.descriptors = ['starting_mean_annual_temperature', 'total_annual_precipitation'] # T and P
d.variables.means = [14.26, 1.59]
d.variables.std_deviations = [0.55, 0.30]
Define the HydroTrend outputs to be used in the study, as well as the statistics to be calculated from them.
In [ ]:
d.responses.response_descriptors = 'Qs_median'
d.responses.response_files = 'HYDROASCII.QS'
d.responses.response_statistics = 'median'
In the stochastic collocation method,
a polynomial expansion is computed for the response function
from four Gauss points in both T and P,
for a total of 16 integration points.
From this expansion, 1000 samples are selected using Latin hypercube sampling to perform the analysis.
The method calculates the PDF and the CDF of Qs using the bins specified by the response_levels keyword.
Variance-based decomposition is turned on to provide additional sensitivity statistics,
such as Sobol' indices.
In [ ]:
d.method.quadrature_order = 4
d.method.sample_type = 'lhs'
d.method.samples = 1000
d.method.seed = 17
d.method.response_levels = [2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5]
d.method.variance_based_decomp = True
HydroTrend requires a set of files to run. They're included in the data directory of the repository containing this example. They can also be obtained directly from the HydroTrend GitHub repository. Set paths to these files with the following statements.
In [ ]:
import os
data_dir = os.path.join(os.getcwd(), 'data')
template_file = os.path.join(data_dir, 'hydrotrend.in.tmpl')
parameters_file = os.path.join(data_dir, 'parameters.yaml')
hypsometry_file = os.path.join(data_dir, 'HYDRO0.HYPS')
The template file provides the configuration file for HydroTrend, but with all parameter values replaced by variables in the form {parameter_name}. The parameters file provides descriptions, ranges, and default values for all of the parameters represented in the template file. The hypsometry file describes the change in elevation along the river's course from source to sea.
From the template and parameters files, we can create an input file that HydroTrend can run. Included in the CSDMS Dakota package is a routine that replaces the variables in the template file with default values from the parameters file. Import this routine and use it to create a HydroTrend input file.
In [ ]:
from dakotathon.plugins.base import write_dflt_file
default_input_file = write_dflt_file(template_file, parameters_file, run_duration=10*365)
print default_input_file
Next, we must replace the default values for the variables for starting_mean_annual_temperature and total_annual_precipitation with variable names for Dakota to substitute into. The CSDMS Dakota package also includes a routine to do this. Import this routine and use it to create a Dakota template file.
In [ ]:
from dakotathon.plugins.base import write_dtmpl_file
dakota_template_file = write_dtmpl_file(template_file, default_input_file, d.variables.descriptors)
print dakota_template_file
Associate the Dakota template file and the hypsometry file with the Dakota instance.
In [ ]:
d.template_file = dakota_template_file
d.auxiliary_files = hypsometry_file
Call the setup method to create files needed by Dakota, then run the experiment.
In [ ]:
d.setup()
d.run()
Check the output. First, the dakota.dat file. It shows the quadrature points, and the median values of Qs calculated at those points.
In [ ]:
%cat dakota.dat
Next, the dakota.out file. At the end of this file, statistics calculated by Dakota are reported.
In [ ]:
%cat dakota.out
The coefficient of variation, $c_v = \sigma / \mu$, can be used as a measure of the uncertainty in the inputs, T and P. For T, $c_v = 1.59/14.26 = 3.9\%$, while for P, $c_v = 0.30/1.59 = 18.9\%$. These uncertainties are propagated through the model by Dakota, which calculates the first four moments of the response, Qs. From these moments, we can also calculate a coefficient of variation for Qs, $c_v = 0.346/4.34 = 8.0\%$.
Sobol' indices are calculated, which provide a measure of sensitivity of Qs to the model inputs. For this study, 75 percent of the variance in Qs is caused by variance in T, while only about 25 percent is caused by P, and with marginal interactions between the two.