Parameter uncertainty estimation of transient storage model using Bayesian inference with formal likelihood based on breakthrough curve segmentation

Journal Pre-proof Parameter uncertainty estimation of transient storage model using Bayesian inference with formal likelihood based on breakthrough curve segmentation Soo Yeon Choi, Il Won Seo, Young-Oh Kim PII:

S1364-8152(19)30734-0

DOI:

https://doi.org/10.1016/j.envsoft.2019.104558

Reference:

ENSO 104558

To appear in:

Environmental Modelling and Software

Received Date: 1 August 2019 Accepted Date: 12 October 2019

Please cite this article as: Choi, S.Y., Seo, I.W., Kim, Y.-O., Parameter uncertainty estimation of transient storage model using Bayesian inference with formal likelihood based on breakthrough curve segmentation, Environmental Modelling and Software (2019), doi: https://doi.org/10.1016/ j.envsoft.2019.104558. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Parameter uncertainty estimation of transient storage model using Bayesian inference with formal likelihood based on breakthrough curve segmentation

1 2 3

Soo Yeon Choi1, Il Won Seo1, and Young-Oh Kim1

4 5 6

1

Dept. of Civil and Environmental Engineering, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul, 08826, Korea.

7 8

Corresponding author: Il Won Seo ([email protected])

9 10

Highlights

11

•

The uncertainty estimation framework of transient storage model is suggested

12

•

Heteroscedasticity of errors is shown in the breakthrough curve

13

•

Formal likelihood is suggested based on mixture of segmented the breakthrough curve

14

•

The formal likelihood yields the prediction interval that captures observations

•

The formal likelihood identifies the behavioral parameters well

15 16 17

Abstract

18

27

The transient storage model (TSM) for the analysis of pollutant mixing in rivers has been hampered by parameter uncertainties due to the equifinality problem. The generalized uncertainty estimation method, which was frequently used to quantify the parameter uncertainty of TSM, has been criticized because this method uses informal likelihood which can cause overestimation of the uncertainty. Thus, in this study, we suggest a Bayesian inference method using a segment mixture (SM) likelihood, which is a formal likelihood based on the mixture distribution of segmented breakthrough curve. The uncertainty estimation was conducted using three synthetic data and the real tracer test data achieved from Uvas Creek in the USA. The results show that the SM likelihood estimated uncertainties of TSM appropriately by correctly representing the error distribution of the TSM and identifying the behavioral parameters.

28 29 30

Keyword: Transient storage model, Parameter uncertainty estimation, Bayesian inference, Formal likelihood, GLUE, Breakthrough curve segmentation

19 20 21 22 23 24 25 26

31

32

1 Introduction

33

In river mixing studies, the one-dimensional Transient Storage Model (TSM) has been frequently used to simulate behaviors of solute transport since it was proposed by Hays in 1966 (Hays, 1966; Fischer et al., 1979; Nordin and Troutman, 1980; Bencala and Walters, 1983; Bencala, 1984; Seo and Maxwell, 1992; Runkel and Chapra, 1993; Cheong and Seo, 2003; Cheong et al., 2007; Camacho and Gonzalez, 2008; Kelleher et al., 2013; Ward et al., 2017). The reason for the popularity of the TSM is that it can accurately reproduce the observed breakthrough curve (BTC), which is characterized by its long tail, via explaining river mixing processes with three important mixing mechanisms: advection, shear dispersion, and storage effect. However, apart from its successful representation, it has been pointed out that the parameters of TSM are highly uncertain to estimate, which would lead to difficulty in characterizing river mixing processes and the application of TSM for prediction (Harvey and Bencala, 1993; Harvey et al., 1996; Wagner and Harvey, 1997; Wagener et al., 2002; Kelleher et al., 2013; Ward et al., 2017).

34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63

In TSM, there are four parameters that represent different river mixing mechanisms such as advection, shear dispersion, and storage effect as shown in Fig. 1. These parameters should be estimated indirectly, since their values are hard to observe physically in natural rivers (Wagener and Gupta, 2005; Kelleher et al., 2013). In general, the parameter combination is determined using inverse modeling in which a parameter combination producing the best-fit result with the observed BTC is selected as the optimal parameter combination. In other words, the parameter combination is usually estimated based on an optimization technique with a certain objective function, which can be formulated with various goodness-of-fit measures (Dennis et al., 1981; Wagner and Gorelick, 1986; Runkel and Broshears, 1991). However, in this parameter estimation process, more than one parameter combination can generally show similarly good levels of fitting to the observed BTC, which has been referred to as the equifinality problem (Beven and Binley, 1992; Kelleher et al., 2013). It is known that this phenomenon results in difficulty in identifying a unique parameter combination of TSM, and consequently results in the parameter uncertainty. In particular, the parameters standing for the storage effect have been known for their higher uncertainties than the other parameters related to the advection and shear dispersion mechanisms in the free-flowing water zone (Wagner and Harvey, 1997; Wagener et al., 2002; Kelleher et al., 2013).

64 65 66

Fig. 1. Schematic sketch of river mixing. (a) Advection, dispersion, and storage effect in natural river. (b) Conceptual model of river mixing defined in TSM.

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93

When estimating parameters of TSM with the inverse modeling, the parameter uncertainty has usually been estimated in local space of the optimized value by calculating the confidence interval of parameter based on the assumption that the parameter estimate is normally distributed (Donaldson and Tryon, 1990; Runkel, 1998; Ward et al., 2017). However, this method of assessing uncertainty can be extremely incorrect for a highly nonlinear model, even though its application is relatively easy (Donaldson and Tryon, 1990). Furthermore, this method cannot report generalized information on uncertainty across the whole parameter space (Ward et al. 2017). In order to overcome this locality problem, Wagener et al. (2002), Camacho and Gonzalez (2008), and Ward et al. (2017) estimated the parameter uncertainty of TSM using the generalized likelihood uncertainty estimation (GLUE) framework, which has been widely applied in the field of rainfall-runoff modeling since Beven and Binley proposed the GLUE in 1992. In principle, the GLUE was developed based on the Bayesian inference in which a posterior distribution of parameter is achieved by combining the prior knowledge of the parameter distribution and the likelihood of describing the newly observed information. However, the GLUE has been regarded as a pseudo-Bayesian inference, because it uses the informal likelihood based on some goodness-of-fit measure, rather than the formal likelihood based on classical probability distribution (Freni and Manniba, 2012). Use of the informal likelihood in GLUE in the estimation of parameter uncertainty has recently been warned against by several previous studies of Mantovan and Todini (2006), Stedinger et al. (2008), and Clark et al. (2011). The previous studies suggested that the formal likelihood based on the probability distribution representing the error structures should be used in their rainfall-runoff modeling examples, since the GLUE with the informal likelihood would overestimate the parameter uncertainty. For this reason, the recent applications have often employed the formal likelihood, rather than the informal likelihood, to estimate the parameter uncertainty of various hydrologic models such as the rainfall-runoff model (McMillan and Clark, 2009; Schaefli et al, 2007; Vrugt et al., 2009; Schoups and Vrugt, 2010; Smith et al, 2010), water

94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116

distribution model (Hutton et al., 2014), and TMDL model (Hantush and Chaudhary, 2014; Chaudhary and Hantush, 2017). The aforementioned criticisms have brought about the reasonable doubt that the parameter uncertainty in TSM can also be unsoundly estimated with informal likelihood in the previous studies (Wagener et al., 2002; Camacho and Gonzalez, 2008; Ward et al, 2017) even though the informal likelihood is relatively simple to implement. In this context, there has been a need to apply the formal likelihood for parameter uncertainty estimation of TSM; however, as of yet, the use of formal likelihood has never been employed in the study of TSM. Meanwhile, the formal likelihood should be carefully formulated to appropriately represent the error structure of a given model, and to prevent the biased uncertainty estimation of parameters (Stedinger et al., 2008; Schoups and Vrugt, 2010). In other words, an ill-formulated formal likelihood might result in a worse estimation of uncertainty than the informal likelihood (Yang et al., 2008; Freni and Mannina, 2012). The previous studies on rainfall-runoff modeling proposed different types of formal likelihood that can capture key statistical characteristics of error structures such as normality, autocorrelation, and heteroscedasticity. For example, Schaefli et al. (2007) and Smith et al. (2010) used the formal likelihood function of mixture distributions to handle the varying variance of the residual errors in high and low discharge observations separately. This heteroscedasticity condition that appeared in the discharge curve can be frequently found in the observed BTC of pollutant studies in natural streams. The BTC usually consists of segments of rising limb, peak, falling limb, and tail, as shown in Fig. 2 (a). Among them, the tail part of BTC tends to show distinct error structure from the other segments of BTC, due to its long extension with low concentration, as shown in Fig. 2 (b). Therefore, capturing these characteristics of BTCs is important in formulating the formal likelihood of TSM.

117

118 119 120 121

Fig. 2. The characteristics of the breakthrough curve. (a) The segmentation of the breakthrough curve that consists of rising limb, peak, falling limb, and tail segments. (b) The example of the residual errors between observed breakthrough curve and the best fitting simulation of TSM.

122 123 124 125 126 127 128

This study, therefore, proposed a Bayesian framework with the formal likelihood that can capture the distinct statistical characteristics of the BTCs. The BTC segmentation was a particular interest of this study, in which the heteroscedasticity in residual errors, as well as the biasedness and autocorrelation, were effectively reflected by adopting the likelihood function based on the finite mixture distribution for segmented BTC. This study also aimed to demonstrate the superiority of the suggested formal likelihood, by comparing it with the informal likelihood that has been applied in the uncertainty estimation of TSM.

129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152

2 Materials and Methods 2.1 Transient storage model TSM was proposed to represent the observed BTC that cannot be reproduced properly by the conventional one-dimensional advection-dispersion (1D-AD) model of Taylor (1954). The observed BTC is usually more skewed than the curve simulated by the 1D-AD model. This is because the 1D-AD model simulates the behavior of solutes without considering the storage effect that exists along the river channel. As shown in Fig. 1, the storage zone is a relatively immobile zone compared to the main free-flowing water zone where the flow is constantly moving. The storage zone is usually made due to irregular channel boundaries such as pools, vegetation, hydraulic structures, and gravel beds. When solutes pass by in the river channel, some portion of solutes is captured by the storage zone, and slowly released back into the main free-flowing water zone after the main cloud has already moved downstream, which leads to the long-tailed and the skewed BTC. This storage effect was taken into account in TSM by adding the equation of the storage zone to the 1D-AD model (Seo and Maxwell, 1992; Cheong and Seo, 2003). In TSM, the river channel is divided into two conceptual areas: the main free-flowing water zone, and the storage zone, as shown in Fig. 1 (b). In the main free-flowing water zone, advection and shear dispersion are the dominant transport mechanisms, while the storage zone represents the portion of the stream that contributes to the storage of pollutants (Runkel, 1998). TSM consists of two governing equations that model the different mixing processes in each zone. Based on the assumptions of steady uniform flow, conservative solute, and completely mixed storage zone, the governing equations can be established as follows (Bencala and Walters, 1983; Runkel and Broshears, 1991; Runkel, 1998):

153 154

155 156

∂C ∂C ∂  ∂C  + UF =  KF  + α (CS − C ) ∂t ∂x ∂x  ∂x  ∂CS A = α F (C − CS ) ∂t AS

(1) (2)

157

where C , U F , K F , and AF are the solute concentration, mean flow velocity, longitudinal

158

dispersion coefficient, and area of the main free-flowing water zone, respectively; CS and AS

159

are the solute concentration and area of the storage zone, respectively; and α is the mass exchange rate between the main free-flowing water zone and the storage zone.

160 161

In the TSM governing equations, as aforementioned, K F , AF , AS , and α are usually

162

treated as the parameters that should be determined by the inverse modeling using the observed

163

BTC, as explained earlier. Among the four parameters, K F and AF stand for the shear dispersion

164

and advection mechanisms, respectively, while AS and α represent the storage effect. In some

165

cases, U F in Eq. (1) is also treated as a model parameter to be estimated. However, in this study,

166

U F was considered as an input variable that could be measured through field study (Seo et al.,

167

2016). Also, although the governing equations of the TSM above can be extended to include additional mechanisms such as lateral flow, sorption/desorption, and chemical reaction, this study took account of only the main mixing mechanisms such as advection, shear dispersion, and storage effect, which results in four parameters to be estimated. In this study, the numerical model of the TSM based on Eqs. (1) and (2) was constructed in MATLAB code using the finite volume method and the Crank-Nicolson method to approximate the spatial and temporal derivatives, respectively (Runkel, 1998).

168 169 170 171 172 173 174 175

2.2 Uncertainty estimation

176

2.2.1 Bayesian inference framework

177 178 179 180

The Bayesian inference is a useful framework for obtaining a probability distribution of a model parameter that can be quantified to a degree of uncertainty. In the Bayesian inference, a posterior probability distribution of the model parameter is updated by multiplying a prior distribution of the model parameter and a likelihood, as shown in the following equation:

181 182

% ) = cf (C % θ) f (θ) f (θ C

(3)

183 184

where θ = {θ1 ,..., θ d } is the parameter vector, i.e. in this study, θ = {K F , AF , AS , α } and

185

d = 4 ; C% = {C%1 ,..., C% T } is the vector of observed concentrations in the main free-flowing water

186

% ) is the posterior probability density function for the model zone at times t = {1,..., T } ; f (θ C

187 188

% parameter vector θ given a set of observed concentration in the main free-flowing water zone C % θ) is the likelihood function to represent the probability of ; c is a normalization constant; f (C

189

% for a given set of model parameters θ ; and f (θ) is the the observed concentration vector C

190

prior probability density function for the model parameter vector θ .

191

Since the likelihood plays an important role in obtaining the posterior distribution using Eq. (3), the selection of an appropriate likelihood is critical in uncertainty estimation (Smith et al., 2010). As aforementioned, types of likelihood can be classified into formal and informal likelihoods. Depending on which type of likelihood is used, the uncertainty estimation method based on the Bayesian inference can be classified into the formal Bayesian inference, and the GLUE. In the GLUE, the informal likelihood is formulated based on a goodness-of-fit measure such as root-mean-square-error (RMSE), residual sum of squares (RSS), and Nash Sutcliffe efficiency (NSE) (Beven and Binley, 1992; Stedinger et al., 2008; Beven and Binley, 2014). In the previous application of the GLUE for the uncertainty estimation of the TSM (Wagener et al., 2002; Camacho and Gonzalez, 2008; Ward et al., 2017), the RMSE-based informal likelihood has been widely used, as follows:

192 193 194 195 196 197 198 199 200 201 202

203

(

 T C% − C t t % f (C θ) = ∑  t =1 T 

)

2

   

− N /2

(4)

204 205 206 207 208 209 210 211

where C% t and Ct are the observed and simulated concentration in the main free-flowing water zone at time t , respectively; T is the number of observed concentration data; and N is a shaping factor that is used to put more weight to the parameters that give a better goodness-of-fit. The formal likelihood can be rigorously formulated based on the probabilistic structure of the errors which are the differences between simulated and observed concentrations. When we assume that the errors are normally distributed, unbiased, independent, and homoscedastic, the formal likelihood can be given as follows:

212

213

% θ) = f (C

 T 2  ∑ εt 1 exp  − t =1 2 2  2σ ε 2πσ ε  

  T % 2    ∑ (Ct − Ct )  1 =  exp  − t =1 2 2σ ε2    2πσ ε      

(5)

214 215

where ε t is the error between the simulated and observed concentrations in the main free-

216

flowing water zone at time t , i.e. ε t = C%t − Ct ; and σε2 is the variance of the errors. However, in

217 218

many cases including TSM, the assumption of unbiasedness, independence, and homoscedasticity of the errors can be easily violated.

219 220 221

In order to solve the problem of biasedness and autocorrelation, the first-order autoregressive model (AR(1)) can be adopted as follows (Sorooshian and Dracup, 1980; Hantush and Chaudhary, 2014).

222 223

ε t − µ = ρ ( ε t −1 − µ ) + ωt , ωt ~ N (0, σ ω2 )

(6)

224 225

where ρ is the first-order autocorrelation coefficient of errors; µ is the mean of errors; ωt is

226

the residual error at time t , which is normally, and independently distributed with zero mean and

227

variance σω2 . By combining Eqs. (5) and (6), the formal likelihood reflecting the biasedness and

228

autocorrelation of errors can be derived as follows:

229

230

% θ) = f (C

2   T ε t − µ − ρ ( ε t −1 − µ ) )  ( ∑  1  exp  − t =1 2σ ω2   2πσ ω2    

(7)

231 232 233 234 235 236 237 238 239 240 241 242 243

In order to rectify the heteroscedasticity of the model residual errors, data transformation, such as the Box-Cox transformation, has been employed in previous studies (Stedinger et al., 2008; Vrugt et al., 2009; Hantush and Chaudhary, 2014). However, this solution can be problematic, since the parameter uncertainty can be overestimated by the data transformation (Schaefli et al, 2007; Freni and Mannina, 2012). Another solution is to employ the finite mixture distribution, which expresses a complex distribution with the weighted sum of normal distributions of different means and variances (Schaefli et al., 2007; Smith et al., 2010). In this study, since the residual error corresponding to the tail segment of BTC tended to show distinct variance from that corresponding to the body segment (i.e. the rising limb, the peak, and the falling limb) of BTC, the segment mixture likelihood, hereafter named the SM likelihood, in which two different distributions corresponding to the body and the tail segments of BTC were combined, was suggested, as follows:

244

% θ) = ( 2πσ 2 )−T1 /2 ( 2πσ 2 )−T2 /2 f (C ω ,1 ω ,2 245

e2  e1 2 2  ε µ ρ ε µ − − − ( ) ( ) (ε t − µ2 − ρ2 (ε t −1 − µ2 ) )  ∑  ∑ t 1 1 t −1 1 t =s t =s  × exp  − 1 − 2 2   2σ ω ,1 2σ ω2 ,2    

(8)

246 247

where subscripts j = 1 and 2 stand for the body and tail segments of BTC, respectively; T j is the

248

duration of the segment j ; s j and e j are the starting and ending times of the segment j , i.e.

249

s1 = 1 , e1 = T1 , s2 = T1 + 1 and e2 = T ; µ j and ρ j are the mean and first-order autocorrelation

250

coefficient of errors of the segment j , respectively; and σ ω2 , j is the variance of the residual errors

251

of the segment j .

252 253 254 255 256 257 258 259 260 261 262 263

2.2.2 Parameter uncertainty estimation In order to perform the parameter uncertainty estimation of TSM based on the Bayesian inference as shown in Eq. (3), TSM should be implemented with a large number of the parameter combinations that were sampled by the sampling technique, since Eq. (3) cannot be implemented in an analytical way. This study employed two sampling techniques according to the type of likelihood. When applying the formal likelihood, a differential evolution adaptive Metropolis (DREAM) sampling technique was used to sample the parameter combinations based on multichain Markov chain Monte Carlo (MCMC) simulation algorithm (Vrugt et al., 2009). Meanwhile, in the case of informal likelihood, a Latin hypercube sampling technique (LHS) was used. In LHS, parameters were sampled from the equally divided partitions of cumulative probability distributions of the parameters (McKay et al., 1979).

265

The flow charts in Fig. 3 explain the simplified procedures of the parameter uncertainty estimation using the DREAM and the LHS. In DREAM, the first step was to determine the prior

266

distributions of the parameters, and then the initial parameter combinations θmk=1 , m = {1,..., M }

267 269

were sampled for M chains from the prior distributions. After simulating the TSM with the initial parameter combinations, the likelihoods of the initial parameter combinations were calculated based on the SM likelihood, Eq. (8). Afterwards, the evolution of chains began by

270

generating the candidate parameter combinations θmcand via the evolutionary algorithm, and by

271

determining to accept the candidate parameter combinations based on the Metropolis acceptance probability pacc . The Metropolis acceptance probability could be calculated as a ratio of the

264

268

272

275

likelihoods of the previous and the present parameter combinations. If the Metropolis acceptance probability was larger than the random number drawn from the uniform distribution on the interval (0,1), then the candidate parameter combination was accepted as the k th parameter

276

m combination of the m th chain, i.e. θ mk = θcand . Otherwise, the previous parameter combination

277

was added to the chain as the k th parameter combination, i.e. θmk = θmk −1 . This evolution process

273 274

278 279 280 281 282 283 284

was performed iteratively until the convergence condition was met, i.e. the Gelman-Rubin  was less than 1.2. The parameter combinations that were achieved convergence diagnostic R after the convergence were used to construct the posterior distribution of parameters. In addition to what is described here, detailed explanation of the DREAM can be found in the studies of Vrugt et al. (2009), and Vrugt (2016). When using the LHS in the case of informal likelihoods, the prior distributions of the parameters were determined, and then n parameter combinations were sampled all at once from

285 286 287 288 289 290 291 292 293 294 295 296

the prior distributions. After the likelihood for each parameter combination was calculated based on the informal likelihood, Eq. (4), a cutoff value was introduced to distinguish behavioral parameter combinations from non-behavioral parameter combinations, in which a behavioral parameter combination is defined as a parameter combination that provides a simulation series that is well fitted to the observation with a high degree of accuracy. Non-behavioral parameter combinations were discarded, and only the behavioral parameter combinations were used to obtain the posterior distributions of parameters. Though the parameters of the error distribution, such as µ j , σ ω2 , j and ρ j in Eq. (8), should be estimated by using the sampling technique as in the model parameters, this study assumed that these values were close to the maximum likelihood estimates (MLE). This study employed the following equations suggested by Chaudhary and Hantush (2017) in order to reduce the dimensionality of the parameters to be sampled.

297 ej

∑ (ε

298

ρˆ j =

t =s j

t −1

− ε t −1 )( ε t − ε t )

ej

∑ (ε

t =s j

t −1

− ε t −1 )

(9) 2

e

j 1 ∑ (ε t − ρ jε t −1 ) T j (1 − ρ j ) t = s j

299

µˆ j =

300

σˆω2 , j =

1 Tj

ej

∑ (ε

t =s j

t

− µˆ j − ρˆ j (ε t −1 − µˆ j ) )

(10) 2

(11)

301 302

where ε t −1 is the sample mean of the one-lagged errors from t = s j − 1 to t = e j − 1 , εt is the

303

sample mean of the present errors from t = s j to t = e j ; and the hat operator denotes the estimator

304

of the corresponding variable.

305 306

307

308 309 310 311

Fig. 3. The estimation procedure of TSM parameter uncertainty. (a) The procedure when using the SM likelihood based on the DREAM sampling technique. (b) The procedure when using the RMSE-based informal likelihood based on the LHS sampling techniques.

312 313 314 315

2.2.3 Prediction uncertainty estimation In Bayesian inference, the posterior distribution of predicted concentrations can be derived after obtaining the posterior distribution of parameters by the following equation (Zellner, 1971):

316 317

% ) = f (P θ, C % ) f (θ C % )dθ f (P C ∫ θ

(12)

318 319

where P = {P1 ,..., PT } is the vector of the predicted concentrations in the main free-flowing

320

% ) is the posterior probability density function for the water zone at times t = {1,..., T } ; f (P C

321

% ) is the posterior probability density function for % ; and f (P θ, C predicted concentration given C

322

% and θ . When using the SM likelihood to estimate the the predicted concentration given C prediction uncertainty based on Eq. (12), the predicted concentrations were calculated by adding error series, ε = {ε 1 ,..., ε T } to the simulated concentrations, C = {C1 ,..., CT } , in order to represent

323 324 325 326 327 328

that the uncertainty of predicted concentrations contained both parameter and model uncertainties explicitly (Stedinger et al., 2008; Vrugt et al., 2009). Herein, the error at time t was generated from the assumed error structure of the SM likelihood, and added to the simulated concentration at time t , as follows:

329 330

 N ( µ1 , σ ω2 ,1 (1 − ρ12 ) ) , when 1 ≤ t ≤ e1  Pt = Ct + ε t , ε t ~  2 2  N ( µ2 , σ ω ,2 (1 − ρ 2 ) ) , when s2 ≤ t ≤ T

(13)

331 332 333 334 335

However, in the case of informal likelihood, the simulated concentration became the predicted concentration, i.e. Pt = C t , since it was assumed that the model uncertainty was implicitly included in the parameter uncertainty estimated by the informal likelihood, so that the model errors were not considered explicitly.

336 337 338 339 340 341 342 343 344

3. Case study 3.1 Synthetic tracer test data Three synthetic data were generated to see whether the suggested SM likelihood, Eq. (8) can spot the true parameters within the 95 % confidence interval of parameters, regardless of the level of parameter uncertainty. Thus three synthetic data were devised according to the Damkohler number (Dal) that represents the relationship between the rate of storage effect and the rate of advection, as follows:

345

α 1 + AF A  L S  Dal =  UF

(14)

346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363

when Dal was one, due to the dominance of advection. Meanwhile, when Dal was larger than 10, CoVs of AS and α increased by more than 10 times, due to the dominance of storage effect. Three types of synthetic data, referred to as Syn1, Syn2, and Syn3, were generated corresponding to Dal of 1, 0.01, and 10, respectively, as shown in Fig. 4. In this study, the procedure to generate synthetic BTC data could be largely divided into two steps. First, BTCs were generated by TSM using the assumed input variables and parameters given in Table 1. These values were determined within the reasonable parameter ranges according to the previous studies, to satisfy the specific Dal values. The numerical parameters, ∆x and ∆t , which represent the length of spatial intervals consisting of the target reach and the time step, respectively, were also selected to satisfy the numerical stability condition based on Courant number (Ward et al., 2017). In the second step, the errors were generated by sampling from the assumed error distributions that had the parameters given in Table 1, and added to generated BTCs to mimic the characteristics of the real data.

Conc. (ppm)

364

where L is the length of the targeted reach. Wagner and Harvey (1997) showed that parameters were estimated with the smallest uncertainty when Dal was close to one, while the parameter uncertainty increased when Dal was much higher or much lower than one, as either the advection or storage effect became to prevail. In detail, when Dal was less than 0.01, the coefficients of variance (CoVs) of AS and α became more than 100 times larger, compared to the situation

365 366 367

Fig. 4. The synthetic breakthrough curves generated for parameter uncertainty estimation.

368

Table 1. Input Variables and Parameters used for Generating Synthetic Data Q U L KF AF AS △x α △t ρ Case Dal 3 2 2 2 (m /s)

(m/s)

(m)

(m /s)

(m )

(m )

(1/s)

Syn1

0.3

0.075

250

1.5

4

2

1×10-4

Syn2

1

0.5

250

1

2

2

Syn3

0.3

0.075

250

2

4

2

σ ω2

(m)

(s)

1

2

20

0.7

0.05 Ct

1×10-5

0.01

5

10

0.7

0.03 Ct

1×10-3

10

2

20

0.7

0.01 Ct

369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390

3.2 Real tracer test data The real tracer test data used in this study was achieved in Uvas Creek, which is a small stream located in Santa Clara County, California, USA (Zand et al., 1976; Avanzino et al., 1984). The Uvas Creek tracer test data was selected for application of the proposed method, because it has been a benchmark data in the study of river mixing since Bencala and Walters (1983) used it for TSM application. Ward et al. (2017) also employed the Uvas Creek data for uncertainty estimation of TSM using the GLUE. As shown in Fig. 5, the tracer test was conducted in a 635 m experimental reach of Uvas Creek on September 26, 1972 (Avanzino et al., 1984). The study reach consists of pool and riffle sequences, and the riverbed is composed of bed materials greater than 4 mm in diameter. During the tracer test, the average daily streamflow was relatively constant at 0.0125 m3/s. As a tracer, chloride was injected at a constant rate for 3 hours at the upstream of the study reach. Before the injection, the background concentration of chloride was measured as 3.7 mg/L. Chloride concentrations were monitored at five sections of UV-S1, UVS2, UV-S3, UV-S4, and UV-S5, which were located at 38, 105, 281, 433, and 619 m downstream from the injection point, respectively, as shown in Fig. 5 (a). Water samples were collected to measure Chloride concentration from 30 minutes before the injection to 27 hours after the injection. Fig. 5 (b) shows the measured Chloride concentration curves at the five sections. For uncertainty estimation, the data achieved in the UV-S1 and UV-S5 sections were used as an upstream boundary condition and downstream observation data, respectively. The modeling scheme was chosen the same as the study of Ward et al. (2017).

391 392

Fig. 5. The description of the tracer test in Uvas Creek. (a) The study reach, Uvas Creek. (b) The measured breakthrough curves at five sections.

393 394 395 396 397 398 399 400 401 402

4. Results 4.1 Application to Synthetic data The SM likelihood based on the mixture distribution, Eq. (8), was applied for uncertainty estimation of the synthesized data. The prior distributions of parameters were assumed as uniform distributions, and the parameter ranges are given in Table 2. The DREAM was applied to sample the parameter combinations. The numbers of sampling required to satisfy the convergence condition were 3,000, 6,500, and 20,800 for Syn1, Syn2, and Syn3 cases, respectively. After the convergence, the respective 10,000 parameter combinations for the three cases were sampled to obtain the posterior distributions.

403 404

Table 2. Sampling Ranges of TSM Parameters for Uncertainty Estimation

Dispersion coefficient (m /s) Main free- flowing channel area (m2)

KF

Syn1 0.01 ~ 5

Sampling range Syn2 Syn3 0.01 ~ 5 0.01 ~ 5

AF

0.01 ~ 10

0.01 ~ 10

0.01 ~ 10

0.01 ~ 1

Storage zone area (m2)

AS

0.01 ~ 5

0.01 ~ 5

0.01 ~ 5

0.01 ~ 1

1×10-6

1×10-7 ~ 1×10-3

1×10-5 ~ 1×10-1

1×10-5 ~ 1×10-3

Parameter

Symbol 2

Mass exchange rate (/s)

α

~ 1×10-2

Uvas Creek 0.01 ~10

405 406 407

The statistics of parameter posterior distributions were calculated as listed in Table 3. The calculated CoVs in Table 3 demonstrate that the parameter uncertainty changed according to the

409

Dal, which was consistent with the study of Wagner and Harvey (1997). When Dal changed from 1 to 0.01, the CoVs of AS and α increased by almost seven and sixty times, respectively,

410

and the CoV of AF decreased by almost eight times. These results indicate that the uncertainty of

411

parameters related to the storage effect significantly increased while the advection related parameter became less uncertain when the advection became dominant. Meanwhile, when Dal increased from 1 to 10, the CoV of AF increased six times due to the dominance of the storage

408

412 413 414 415 416 417

effect. The uncertainty of the other parameters, especially α , also increased greatly since the mass exchange rate increased to the level of equilibrium condition in which the storage mechanism can be explained by the longitudinal dispersion, leading to difficulties in identifying the parameters (Harvey et al., 1996; Wagner and Harvey, 1997).

418 419 420

Table 3. The Statistical Summary of Posterior Distribution of Parameters and the Coverage Rate of the Prediction Uncertainty for the Synthetic Data Observation Case

Parameter

KF (m2/s) Syn1 (Dal = 1)

Syn2 (Dal = 0.01)

Syn3 (Dal = 10)

True value

Parameter uncertainty Simulation Percent Expected error* 95% CI value (%)

1.5

1.532

2.13

1.362~1.702

0.056

2

4

4.026

0.65

3.907~4.152

0.015

2

AS (m )

2

1.927

3.65

1.618~2.262

0.084

α (/s)

1×10-4

0.98×10-4

2.00

8.42×10-5 ~1.12×10-4

0.073

KF (m2/s)

1

0.956

3.70

0.807~1.116

0.083

AF (m )

2

1.995

0.20

1.989~2.002

0.002

AS (m2)

2

2.527

26.35

0.175~4.863

0.564

α (/s)

1×10-5

0.69×10-5

31.00

1.71×10-7 ~1.77×10-5

4.500

KF (m2/s)

2

2.144

7.20

1.639~2.892

0.151

AF (m2)

4

3.805

4.88

3.129~4.507

0.091

2

2.219

10.95

1.617~2.859

0.138

AF (m )

2

2

AS (m )

-3

0.72×10 ~3.72×10-3 * Percent error (%) = (True value-Expected value)/True value×100

α (/s)

1×10-3

1.48×10-3

** CoV = Standard deviation/Expected value 421

CoV**

48.00

0.530

Prediction uncertainty Coverage rate (%)

93.3

93.5

94.2

425

In order to compare the true and the expected values of the posterior distributions, the percent errors were calculated as shown in Table 3. In the case of Syn1, the percent errors of four parameters were less than 4 %, which showed that the expected values were estimated similarly to the true values of parameters. The case of Syn2 demonstrated that the percent errors of AS and

426

α significantly increased while those values of AF and KF were less than 4 %. Meanwhile, in

422 423 424

427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459

the case of Syn3, the percent errors of all parameters increased in general, and in particular, the expected values of α deviated substantially from the true values. The trend of the percent errors was in accord with that of the parameter uncertainty, in which the highly uncertain parameters showed large percent errors. However, despite the discrepancy between the true and estimated values of parameters in Syn2 and Syn3, it was shown that the 95 % confidence interval of parameters in all cases included the true values within the ranges. This indicated that the 95 % confidence interval estimated by using the SM likelihood could capture the true parameter value well, regardless of the level of parameter uncertainty. In addition, the uncertainty estimation was validated by calculating the coverage rate, which could be calculated as the percentage of the observations that are included within the 95 % prediction interval (Yadav et al., 2007; Smith and Marshall, 2010). Herein, the ideal value of the coverage rate corresponded to the value of the desired interval, 95 % (Smith et al., 2010). The coverage rates were calculated as 93.3, 93.5, and 94.2 % for Syn1, Syn2, and Syn3, respectively, as shown in Table 3. The calculated coverage rates could approximate the ideal values, which suggests that the SM likelihood can perform the uncertainty estimation of TSM well, when applying it to the various synthetic data having different levels of uncertainty. 4.2 Application to the real tracer test in Uvas Creek In application to the Uvas Creek tracer test data, besides the SM likelihood, four types of RMSE-based informal likelihood were also applied. Among them, two RMSE-based informal likelihoods were applied in particular, to compare the results of parameter uncertainty estimation of the SM likelihood with those of the RMSE-based likelihoods that were employed in the previous study of Ward et al. (2017): ‘RMSE-1-1’, which had the cutoff value of the top 1 % and the shaping factor of 1; and ‘RMSE-10-1', which had the cutoff value of the top 10 % and the shaping factor of 1. Besides, two additional RMSE-based informal likelihoods were applied to investigate whether the capacity of the informal likelihood could be improved by controlling the shaping factor of the RMSE-based informal likelihood. These informal likelihoods, ‘RMSE-1-6’ and ‘RMSE-10-6’ had the same shaping factors of 6, but had the cutoff values of the top 1 and 10 %, respectively. The prior distributions of parameters were assumed as uniform distributions spanning the parameter ranges, as shown in Table 2. The parameter ranges used in this study were determined based on the previous study of Ward et al. (2017). From the dotty plots of parameters and the corresponding RMSE that were estimated by Ward et al. (2017), we could narrow down the ranges of parameters, by discarding the non-behavioral parameter ranges where the distribution

460 461 462 463 464 465 466

of scatter points appears flat. When DREAM sampling technique was employed to use the SM likelihood, the convergence was fulfilled when the number of sampled parameters reached 48,000. Additional parameter combinations of 10,000 were sampled after the convergence was met. The number of sampling necessary for the convergence was larger in this case than the case of synthetic data, because there could be more uncertain sources in the data of real tracer test, than in the synthetic data sets. Meanwhile, LHS sampling technique was used to sample the 50,000 parameters, in order to apply the RMSE-based informal likelihoods.

467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491

4.2.1 Estimation of parameter uncertainty By applying the suggested SM likelihoods and the four RMSE-based informal likelihoods, we obtained the cumulative distributions of the four parameters of TSM. From the cumulative distributions, the statistics of the four parameters were calculated, and plotted in boxplots, as shown in Fig. 6. In the boxplots, the center vertical line in the box is the median of the data; the left and right of the box are the 25th and 75th percentiles, respectively; and the ends of the whiskers are the 2.5th and 97.5th percentiles, so that the range between whiskers shows the 95 % confidence interval. This figure shows that the uncertainties of the TSM parameters estimated by the SM likelihood were smaller than those estimated by any type of the informal likelihoods, which was also reported in the previous studies of hydrological modeling (Stedinger et al., 2008; Vrugt et al., 2009). In the cases of informal likelihoods, as the cutoff values changed from the top 1 % to the top 10 %, the parameter uncertainty increased as the 95 % confidence interval got wider. These results happened because the number of parameters that were treated as the behavioral parameters increased, which resulted in wider confidence intervals. On the other hand, when the shaping factor increased from 1 to 6, the parameters with smaller likelihood became neglected in constructing the confidence intervals by emphasizing the parameters with the larger likelihoods, which led to the decrease of the 95 % confidence intervals of the parameters. As with these results, when using the informal likelihood, the parameter uncertainty can be estimated differently, depending on the cutoff values and the shaping factors that are usually determined subjectively. Another notable result on the parameter uncertainty is shown in the statistical summary listed in Table 4. Among the four parameters, the least uncertain parameter was found to be AF , which showed the smallest CoV, regardless of the type of likelihood. In the meantime, two parameters of the storage effect, AS and α , were revealed to be the most uncertain parameters,

493

of which the CoV showed different values ranging from 0.267 to 2.138 depending upon the type of likelihood. Among the five likelihoods, the SM likelihood produced the smallest CoVs for AS

494

and α . The uncertainty of K F was found to show CoV ranging from 0.248 to 1.530.

492

495 496 497

Table 4 also lists the Dal calculated using the expected values of the parameters, in which Dal calculated by SM likelihood was less than unity while Dal of all RMSE-based informal likelihood ranged from 1.11 to 10.93. This result indicates that the SM likelihood characterized

498 499 500 501 502 503 504 505 506 507 508 509 510

the mixing processes in the Uvas Creek as the advection-dominated mechanism; whereas, the RMSE-based informal likelihood interpreted the mixing process as a storage effect-dominated process. These results demonstrate that the selection of likelihood type was significantly important in characterizing the mixing characteristics of natural rivers based on parameter uncertainty estimation. In particular, when using the informal likelihoods, subjective judgments, which were taken to determine the form of the informal likelihood such as cutoff value and shaping factor, could result in the incorrect representation of the mixing characteristics of natural rivers.

511 512 513 514

Fig. 6. The boxplots of TSM parameters estimated by the formal and the RMSE-based informal likelihoods. (a) Longitudinal dispersion coefficient. (b) Area of the main free-flowing water zone. (c) Area of the storage zone. (d) Mass exchange rate.

515 516

Table 4. The Statistical Summary of Posterior Distribution of Parameters and the Calculated Dal for Uvas Creek Case Statistics

Parameter 2

Expected value

CoV

SM

RMSE-10-6

KF (m /s)

0.241

1.001

2.581

0.352

0.614

AF (m2)

0.401

0.445

0.440

0.405

0.409

AS (m2)

0.692

0.470

0.356

0.680

0.637

α (/s)

2.67×10

6.94×10

2.39×10

3.70×10

6.35×10-5

KF (m2/s)

-5

-4

-5

0.248

0.575

0.680

0.930

1.530

0.018

0.124

0.321

0.069

0.136

2

0.320

0.622

0.706

0.267

0.349

0.269 0.79

2.138 2.79

1.170 10.93

1.193 1.11

2.031 1.98

AF (m )

α (/s) -

-5

2

AS (m ) Dal

Type of likelihood RMSE-1-1 RMSE-10-1 RMSE-1-6

517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541

4.2.2 Validation based on model prediction uncertainty In this study, to analyze the model prediction uncertainties by the SM and RMSE-based informal likelihoods, the 95 % prediction interval of each likelihood was estimated. Herein, in addition to coverage rates, the mean widths of the 95 % prediction intervals were also calculated. The mean width was obtained by averaging the width of the 95 % prediction interval over time, which can stand for the size of estimated uncertainty of the prediction results. Table 5 summarizes the estimated coverage rates and mean widths. The uncertainty of prediction results generally followed the trend of the parameter uncertainty according to the type of likelihood. Among the informal likelihoods, RMSE-1-1, RMSE-10-1, and RMSE-10-6, which showed larger parameter uncertainty than the SM likelihood, also presented wider 95 % prediction intervals than the SM likelihood. On the other hand, the mean width of RMSE-1-6 was similar to that of the SM likelihood despite the larger parameter uncertainty. This was because the predicted concentration of SM likelihood was estimated by explicitly adding the model error to the simulated concentration while that of RMSE-1-6 was not, as shown in Eq. (13). When comparing the four informal likelihoods with each other, the mean width increased as the cutoff value and the shaping factor varied from the top 1 % to the top 10 %, and from 6 to 1, respectively, which reflected the trend of the parameter uncertainty. In general, it can be thought that the wider prediction interval leads to the increase of coverage rate. However, the four informal likelihoods yielded lower coverage rates than the SM likelihood despite their wider or similar mean widths. The coverage rates of the four informal likelihoods were much less than the ideal coverage rate of 95 %. In particular, the RMSE-10-1 showed the largest mean width, but the lowest coverage rate. Given that the mean widths were large but the coverage rates were low in the informal likelihoods, it seems that the informal likelihoods estimated the posterior distributions of the parameters incorrectly in deviated

542 543 544 545 546 547

directions from the spaces where the true behavioral parameters existed. When applying the SM likelihood, the coverage rate increased to 92.5 %, which was close to 95 %, although it yielded the narrower 95 % prediction interval than the informal likelihoods. This result meant that the SM likelihood could result in more appropriate uncertainty estimation of the parameter and the prediction because the 95 % prediction interval of the SM likelihood covered the observation well, and did so tightly, without overestimation.

548 549

Table 5. Validation Results of the Prediction Uncertainty for Uvas Creek Case Case UV1 UV2 UV3 UV4 UV5

Type of likelihood

Formula

SM

Eq. (8)

RMSE -1-1 RMSE -10-1 RMSE -1-6 RMSE -10-6

Shaping Cutoff factor value 1 1

Eq. (4) 6 6

Top 1% Top 10% Top 1% Top 10%

Coverage rate (%) Total Body Tail

Mean width (ppm) Total Body Tail

92.5

100

85.2

1.2

1.7

0.7

86.8

100

74.1

1.5

2.1

1.0

75.5

80.8

70.4

1.9

2.7

1.2

81.1

100

63.0

1.1

1.5

0.8

83.0

100

66.7

1.5

2.1

0.9

550 551

5. Discussion

552

The validation in the previous section showed that the SM likelihood yielded the better uncertainty estimation of TSM, compared to the RMSE-based informal likelihood. The reasons for the successful application of the SM likelihood can be summarized as (1) appropriate assumptions of the error distribution; (2) better reflection of information in the entire segments of BTC; and (3) sound judgment on the behavioral parameters.

553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568

The SM likelihood could result in successful uncertainty estimation, thanks to the proper assumptions of the error distribution. The accurate assumptions of the error distribution are important for robust uncertainty estimation since false assumptions can lead to erroneous uncertainty estimation (Stedinger et al., 2008; Smith et al., 2010). Fig. 7 shows the error diagnostics for two formal likelihoods representing different characteristics of the error, in which the NS (not-segmented) likelihood of Eq. (7) assumed the homoscedastic residual error distribution along the whole BTC, while the SM likelihood of Eq. (8) segmented the BTC and assumed the heteroscedasticity of the residual error distribution between the body and tail segments. Figs. 7 (a) and (b) show the error diagnostics for the NS likelihood. Herein, quantilequantile (QQ) plots were used to check whether the residual errors follow the assumed distribution. In the case of NS likelihood, the variance of residual errors, as shown in Fig. 7 (a), fluctuated quite substantially for different values of concentration of the BTC, which implied that

569 570 571 572 573 574 575 576 577

the residual errors did not follow the assumptions of homoscedasticity. The QQ plot in Fig. 7 (b) also shows that the residual errors did not satisfy the assumed distribution either by not following the linear trend of the dotted line. However, the results of SM likelihood revealed that the heteroscedasticity was successfully rectified by dividing the BTC into the body and tail segments, as shown in Figs. 7 (c) and (e). Not only that, the assumption of the distribution of the residual errors became more satisfied by showing the more linear trend of the QQ plots, as shown in Figs. 7 (d) and (f). The error diagnostics showed that it was rational to divide the BTC into the body and tail segment, and to use the different distribution of the errors for each segment as the SM likelihood adopted.

601

Table 5 demonstrates that the SM likelihood generated higher coverage rate in the model prediction compared to the RMSE-based informal likelihood by reflecting the information of the entire segments of BTC. This table shows that the coverage rates of SM likelihood for both body and tail segments were higher than those of RMSE-based informal likelihoods, even though it was seen that the 95 % prediction intervals of both likelihoods missed a lot of observations in the tail segment. In particular, the prediction intervals estimated by the RMSE-based informal likelihoods missed the larger number of the observations in the tail segment. The failures of the RMSE-based informal likelihoods in explaining the observation of the tail segment can be attributed to the property of RMSE, in which it tended to focus on fitting the higher concentration in the peak of BTC while ignoring the information contained in the low concentration in the tail segment (Kelleher et al., 2013). Therefore, the RMSE-based informal likelihoods could not explain well the observation in the tail segment. On the other hand, the SM likelihood could capture more observations in the tail segment than the RMSE-based informal likelihood, which implied that the SM likelihood was better able to reflect the information in the tail segment. The better performance of the SM likelihood compared to the RMSE-based informal likelihood was manifested by estimating the parameter uncertainty using the tail-truncated BTC data which corresponded to the data of only the body segment as shown in Fig. 5 (b). Fig. 8 shows the boxplots of the posterior distributions of the parameters that were estimated by assuming that the tail of BTC was truncated during the tracer test. In this case, the parameter uncertainty estimation was performed using the two likelihoods, the SM and the RMSE-10-1 likelihood, and the results were compared with the results of two likelihoods using the original whole BTC which are shown in Fig. 6. As a result, when the parameter uncertainty was estimated by the SM likelihood, there were significant changes in the medians and the 95 %

602

confidence intervals of the storage zone effect related parameters, AS and α . However, when

603

using the RMSE-10-1 likelihood, the medians and the 95 % confidence intervals were not

604

varying, regardless of the existence of the data of the tail segment. This insensitivity of AS and

605

α to the absence of the data in the tail segment when using the RMSE-10-1 was unreasonable,

578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600

606 607

since the data in the tail segment contained significant information about the storage effect compared to the body segment of the BTC. Wagner and Harvey (1997), Wagener et al. (2002),

608 609 610 611

and Kelleher et al. (2013) have demonstrated that the parameters related to the storage effect, AS and α were the most sensitive to the data of the tail segment via the sensitivity analysis. Thus, this result demonstrates that the use of the SM likelihood estimates the parameter uncertainty accurately, by properly reflecting the information contained in the tail segment of BTC.

619

The robust uncertainty estimation of the SM likelihood could also be attributed to the ability to distinguish the behavioral parameters from the non-behavioral parameters. On the other hand, the RMSE-based informal likelihood misjudged the non-behavioral parameter as behavioral parameter due to the predetermined cutoff value, as indicated by the wider confidence intervals of parameters and larger mean widths compared to the SM likelihood. This argument can be corroborated by the interaction plots shown in Fig. 9. Fig. 9 (a) shows the interaction of the parameters that were included in the 95% confidence intervals estimated by the SM likelihood. In these plots, the clear correlation relationships appear in the relationships between

620

KF , AF , and α . In particular, KF and α had a negative correlation since these two parameters

612 613 614 615 616 617 618

621 622 623 624 625 626 627 628 629

are in a relationship of trade-off interaction. This result clearly demonstrates that the mixing in natural rivers happens by two causes, the shear flow and storage zone effect, which in turn results in a trade-off between the dispersion due to the shear flow, and mixing due to the storage zone effect. Meanwhile, Fig. 9 (b) presents the interaction relationships between the parameters that were included in the 95% confidence intervals when applying the RMSE-10-1. Unlike the case of SM likelihood, no clear relationships were shown. This is because the non-behavioral parameters were selected as the behavioral parameters, which obscured the obvious interaction relationships among the parameters.

630 631 632 633 634 635 636 637 638 639 640 641

Fig. 7. The diagnostic tests on the assumptions on the error distribution. (a) Scatter plot of residual errors of non-segmented BTC when using NS likelihood. (b) QQ plot when using NS likelihood. (c) Scatter plot of residual errors of body segment when using SM likelihood. (d) QQ plot of body segment when using SM likelihood. (e) Scatter plot of residual errors of tail segment when using SM likelihood. (f) QQ plot of tail segment when using SM likelihood.

642 643 644 645 646 647 648 649 650 651

Fig. 8. The boxplots of TSM parameters estimated using the whole BTC and the tail-truncated BTC. (a) Longitudinal dispersion coefficient. (b) Area of the main free-flowing water zone. (c) Area of the storage zone. (d) Mass exchange rate.

652 653 654 655 656 657

Fig. 9. Interaction plots between paired behavioral parameters. (a) The black scatter dots are behavioral parameters determined by the SM likelihood. The red line is the best-fitting linear regression line. (b) The black scatter dots are behavioral parameters determined by the RMSE10-1 likelihood. The red line is the best-fitting linear regression line.

658

5. Conclusion

659

In this study, the Bayesian inference with formal likelihood was applied to estimate the parameter uncertainty of TSM, and the result was compared with that of informal likelihood. This study proposed the segment mixture (SM) likelihood, based on the finite mixture distribution that consists of two error distributions that correspond to the body and tail segments of BTC, respectively. The proposed SM likelihood was tested by applying it to the synthetic tracer test data and the real stream data of the Uvas Creek, and by comparing the uncertainty estimation results with the results of the informal likelihood.

660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696

The results when applying the SM likelihood to the three synthetic data showed that the SM likelihood could capture the true parameters within the estimated 95 % confidence intervals regardless of the level of parameter uncertainty. The properness of the SM likelihood was also validated by calculating the coverage rates of 95 % prediction intervals. The coverage rates of all synthetic data were calculated as close to the ideal value of 95 %. Similarly, in the case of Uvas Creek, the SM likelihood yielded 95 % prediction interval which strictly and well covered the observation with the coverage rate of about 92.5 %. Therefore, the SM likelihood was found to perform appropriately for the parameter and prediction uncertainty of TSM. In contrast, the RMSE-based informal likelihoods estimated the parameter uncertainty incorrectly when compared to the SM likelihood. All of the RMSE-based informal likelihoods gave much less coverage rates than the 95 %, although these likelihoods yielded larger mean widths than the SM likelihood. The successful uncertainty estimation of SM likelihood compared to the RMSE-based likelihoods could be attributed to three properties of the SM likelihood. First, the SM likelihood represented the error characteristics of BTC for TSM appropriately as shown in the error diagnostic tests. In particular, the SM likelihood could represent the heteroscedasticity of the residual errors by the segmenting the BTC into the body and tail parts. Second, the SM likelihood explained the observation along the BTC well by combining comprehensive information from the whole segments of BTC such as the body and tail segments. On the other hand, the RMSE-based informal likelihoods could not employ the information of the tail segments which is useful in estimating the storage effect-related parameters. When comparing the parameter uncertainties estimated by the SM and RMSE-based likelihoods using the tailtruncated BTC, the 95 % confidence intervals of the parameters estimated by the RMSE-based likelihood were insensitive to the nonexistence of the data in the tail segment, while the SM likelihood was not. Also, the SM likelihood was successful in distinguishing the behavioral parameter and non-behavioral parameter, and showed the clear correlations between the parameters in the interaction plots. Meanwhile, the RMSE-based informal likelihood tended to misclassify the non-behavioral parameters as behavioral parameters, which could lead to the incorrect parameter and prediction uncertainty estimation. When it comes to the uncertainty estimation of TSM, it is important to choose appropriate likelihood since the selection of inappropriate likelihood can give rise to incorrect

697 698 699

interpretation of the river mixing characterization, and degradation of the accurate prediction of the pollutants mixing. In that sense, the SM likelihood that was suggested in this study can be a good option for the likelihood in uncertainty estimation of the TSM.

700 701

Software/data availability

702

The TSM developed in MATLAB is available from the link http://ehlab.snu.ac.kr/ and the MATLAB toolbox of DREAM (Vrugt, 2016) is available from the link https://faculty.sites.uci.edu/jasper/software/. The tracer data of Uvas Creek was taken from Avanzino et al. (1984) and the synthetic data can be downloaded from the link http://ehlab.snu.ac.kr/.

703 704 705 706 707 708

Acknowledgments

709

This research is supported by Korea Ministry of Environment (MOE) as “Chemical Accident Response R&D Program (2018001960001)”, and the BK21 PLUS research program of the National Research Foundation of Korea. This research work was conducted at the Institute of Engineering Research and Institute of Construction and Environmental Engineering in Seoul National University, Seoul, Korea.

710 711 712 713 714 715

References

716

Avanzino, R. J., Zellweger, G. W., Kennedy, V. C., Zand, S. M., Bencala, K. E., 1984. Results of a solute transport experiment at Uvas Creek, September 1972 (No. 84-236). Menlo Park, CA: US Geological Survey.

717 718 719 720 721 722 723 724 725 726 727 728

Bencala, K. E., 1984. Interactions of solutes and streambed sediment, 2, A dynamic analysis of coupled hydrologic and chemical processes that determine solute transport. Water Resour. Res. 20, 1804–1814. http://doi.org/10.1029/WR020i012p01804. Beven, K., Binley, A., 1992. The future of distributed models: Model calibration and uncertainty prediction. Hydrol. Process. 6, 279-298. http://doi.org/10.1002/hyp.3360060305. Beven, K., Binley, A., 2014. GLUE: 20 years on. Hydrol. Process. 28, 5897-5918. http://doi.org/10.1002/hyp.10082. Bencala, K. E., Walters, R. A., 1983. Simulation of solute transport in a mountain pool‐and‐ riffle stream: A transient storage model. Water Resour. Res. 19, 718-724. http://doi.org/10.1029/WR019i003p00718.

729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763

Camacho, L. A., González, R. A., 2008. Calibration and predictive ability analysis of longitudinal solute transport models in mountain streams. Environ. Fluid. Mech. 8, 597604. http://doi.org/10.1007/s10652-008-9109-0. Chaudhary, A., Hantush, M. M., 2017. Bayesian Monte Carlo and maximum likelihood approach for uncertainty estimation and risk management: Application to lake oxygen recovery model. Water Res. 108, 301-311. http://doi.org/10.1016/j.watres.2016.11.012. Cheong, T. S., Seo, I. W., 2003. Parameter estimation of the transient storage model by a routing method for river mixing processes. Water Resour. Res. 39, HWC11-HWC111. http://doi.org/10.1029/2001WR000676. Cheong, T. S., Younis, B. A., Seo, I. W., 2007. Estimation of key parameters in model for solute transport in rivers and streams. Water Resour. Manag. 21, 1165-1186. http://doi.org/10.1007/s11269-006-9074-7. Clark, M. P., Kavetski, D., Fenicia, F., 2011. Pursuing the method of multiple working hypotheses for hydrological modeling. Water Resour. Res. 47. http://doi.org/10.1029/2010WR009827. Dennis, J. E., Jr., Gay, D. M., Walsh, R. E., 1981. An adaptive nonlinear least-squares algorithm. ACM Trans. Math Softw. 7, 348-368. http://doi.org/10.1145/355958.355965. Donaldson, J. R., Tryon, P. V., 1990. User’s Guide To STARPAC-The standards, time series, and regression package: National Institute of Standards and Technology Internal Report NBSIR 86-3448. Technical report, National Institute of Standards and Technology. Fischer, H. B., List, J. E., Koh, C. R., Imberger, J., Brooks, N. H., 1979. Mixing in Inland and Coastal Waters, San Diego, CA: Academic. Freni, G., Mannina, G., 2012. Uncertainty estimation of a complex water quality model: The influence of box-cox transformation on Bayesian approaches and comparison with a nonBayesian method. Phys. Chem. Earth. 42, 31-41. http://doi.org/10.1016/j.pce.2011.08.024. Hantush, M. M., Chaudhary, A., 2014. Bayesian framework for water quality model uncertainty estimation and risk management. J. Hydrol. Eng. 19, 04014015. http://doi.org/10.1061/(ASCE)HE.1943-5584.0000900. Harvey, J. W., Bencala, K. E., 1993. The effect of streambed topography on surface‐subsurface water exchange in mountain catchments. Water Resour. Res. 29, 89-98. http://doi.org/10.1029/92WR01960. Harvey, J. W., Wagner, B. J., Bencala, K. E., 1996. Evaluating the reliability of the stream tracer approach to characterize stream-subsurface water exchange. Water Resour. Res. 32, 2441-2451. http://doi.org/10.1029/96WR01268.

764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799

Hays, J. R., 1966. Mass transport mechanisms in open channel flow, Ph.D. dissertation, Department of Civil Engineering, Vanderbilt University, Nashville, TN. Hutton, C. J., Kapelan, Z., Vamvakeridou-Lyroudia, L., Savić, D., 2014. Application of formal and informal bayesian methods for water distribution hydraulic model calibration. J. Water Res. Plan. Man. 140, 04014030. http://doi.org/10.1061/(ASCE)WR.19435452.0000412. Kelleher, C., Wagener, T., McGlynn, B., Ward, A. S., Gooseff, M. N., Payn, R. A., 2013. Identifiability of transient storage model parameters along a mountain stream. Water Resour. Res. 49, 5290-5306. http://doi.org/10.1002/wrcr.20413. Mantovan, P., Todini, E., 2006. Hydrological forecasting uncertainty assessment: Incoherence of the GLUE methodology. J. Hydrol. 330, 368-381. http://doi.org/10.1016/j.jhydrol.2006.04.046. McKay, M. D., Beckman, R. J., Conover, W. J., 1979. Comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics. 21, 239-245. http://doi.org/10.1080/00401706.1979.10489755. McMillan, H., Clark, M., 2009. Rainfall-runoff model calibration using informal likelihood measures within a markov chain monte carlo sampling scheme. Water Resour. Res. 45. http://doi.org/10.1029/2008WR007288. Nordin, C. F., Troutman, B. M., 1980. Longitudinal dispersion in rivers: The persistence of skewness in observed data. Water Resour. Res. 16, 123-128. Runkel, R. L., 1998. One-dimensional Transport with Inflow and Storage (OTIS): a solute transport model for streams and rivers, Water-Resources Investigation Report 98-4018, Reston, VA: US Geological Survey. Runkel, R. L., Broshears, R. E., 1991. One-dimensional transport with inflow and storage (OTIS): a solute transport model for small streams. CADSWES, Center for Advanced Decision Support for Water and Environmental Systems, Department of Civil Engineering, University of Colorado. Runkel, R. L., Chapra, S. C., 1993. An efficient numerical solution of the transient storage equations for solute transport in small streams. Water Resour. Res. 29, 211-215. http://doi.org/10.1029/92WR02217. Schaefli, B., Talamba, D. B., Musy, A., 2007. Quantifying hydrological modeling errors through a mixture of normal distributions. J. Hydrol. 332, 303-315. http://doi.org/10.1016/j.jhydrol.2006.07.005. Schoups, G., Vrugt, J. A., 2010. A formal likelihood function for parameter and predictive inference of hydrologic models with correlated, heteroscedastic, and non-Gaussian errors. Water Resour. Res. 46. http://doi.org/10.1029/2009WR008933.

800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834

Seo, I. W., Cheong, T. S., 2001. Moment-based calculation of parameters for the storage zone model for river dispersion. J. Hydraul. Eng. 127, 453-465. Seo, I. W., Choi, H. J., Kim, Y. D., Han, E. J., 2016. Analysis of two-dimensional mixing in natural streams based on transient tracer tests. J. Hydraul. Eng. 142. http://doi.org/10.1061/(ASCE)HY.1943-7900.0001118. Seo, I. W., Maxwell, W. H. C., 1992. Modeling low-flow mixing through pools and riffles. J. Hydraul. Eng. 118, 1406-1423. http://doi.org/10.1061/(ASCE)07339429(1992)118:10(1406). Smith, T. J., Marshall, L. A., 2010. Exploring uncertainty and model predictive performance concepts via a modular snowmelt-runoff modeling framework. Environ. Model. Softw 25, 691-701. http://doi.org/10.1016/j.envsoft.2009.11.010. Smith, T., Sharma, A., Marshall, L., Mehrotra, R., Sisson, S., 2010. Development of a formal likelihood function for improved Bayesian inference of ephemeral catchments. Water Resour. Res. 46. http://doi.org/10.1029/2010WR009514. Sorooshian, S., Dracup, J. A., 1980. Stochastic parameter estimation procedures for hydrologie rainfall‐runoff models: Correlated and heteroscedastic error cases. Water Resour. Res. 16, 430-442. Stedinger, J. R., Vogel, R. M., Lee, S. U., Batchelder, R., 2008. Appraisal of the generalized likelihood uncertainty estimation (GLUE) method. Water Resour. Res. 44. http://doi.org/10.1029/2008WR006822. Taylor, G. I., 1954. The dispersion of matter in turbulent flow through a pipe, Proc. R. Soc. London, Ser. A. 223, 446–468. Vrugt, J. A., 2016. Markov chain Monte Carlo simulation using the DREAM software package: Theory, concepts, and MATLAB implementation. Environ. Model. Softw. 75, 273-316. https://doi.org/10.1016/j.envsoft.2015.08.013. Vrugt, J. A., Ter Braak, C. J., Gupta, H. V., Robinson, B. A., 2009. Equifinality of formal (DREAM) and informal (GLUE) Bayesian approaches in hydrologic modeling?. Stochastic Environ. Res. Risk Assess. 23, 1011-1026. http://doi.org/10.1007/s00477-0080274-y. Wagener, T., Camacho, L. A., Wheater, H. S., 2002. Dynamic identifiability analysis of the transient storage model for solute transport in rivers. J. of Hydroinform. 4, 199-211. https://doi.org/10.2166/hydro.2002.0019. Wagener, T., Gupta, H. V., 2005. Model identification for hydrological forecasting under uncertainty. Stochastic Environ. Res. Risk Assess. 19, 378-387. http://doi.org/10.1007/s00477-005-0006-5.

835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854

Wagner, B. J., Gorelick, S. M., 1986. A Statistical Methodology for Estimating Transport Parameters: Theory and Applications to One‐Dimensional Advectivec‐Dispersive Systems. Water Resour. Res. 22, 1303-1315. Wagner, B. J., Harvey, J. W., 1997. Experimental design for estimating parameters of rate‐ limited mass transfer: Analysis of stream tracer studies. Water Resour. Res. 33, 17311741. Ward, A. S., Kelleher, C. A., Mason, S. J., Wagener, T., McIntyre, N., McGlynn, B., Payn, R. A., 2017. A software tool to assess uncertainty in transient-storage model parameters using Monte Carlo simulations. Freshw. Sci. 36, 195-217. http://doi.org/ 10.1086/690444. Yadav, M., Wagener, T., Gupta, H., 2007. Regionalization of constraints on expected watershed response behavior for improved predictions in ungauged basins. Adv. Water Resour. 30, 1756-1774. http://doi.org/10.1016/j.advwatres.2007.01.005. Yang, J., Reichert, P., Abbaspour, K. C., Xia, J., Yang, H., 2008. Comparing uncertainty analysis techniques for a SWAT application to the Chaohe Basin in China. J. Hydrol. 358, 1-23. https://doi.org/10.1016/j.jhydrol.2008.05.012. Zand, S. M., Kennedy, V. C., Zellweger, G. W., Avanzino, R. J., 1976. Solute transport and modeling of water quality in a small stream. J. Res. US Geol. Surv. 4, 233-240. Zellner, A., 1971. An introduction to Bayesian inference in econometrics. NY: Wiley.