- Email: [email protected]
Joint analysis of input and parametric uncertainties in watershed water quality modeling: A formal Bayesian approach
Accepted Manuscript
Joint analysis of input and parametric uncertainties in watershed water quality modeling: a formal Bayesian approach Feng Han , Yi Zheng PII: DOI: Reference:
S0309-1708(17)30982-X 10.1016/j.advwatres.2018.04.006 ADWR 3126
To appear in:
Advances in Water Resources
Received date: Revised date: Accepted date:
30 October 2017 11 April 2018 11 April 2018
Please cite this article as: Feng Han , Yi Zheng , Joint analysis of input and parametric uncertainties in watershed water quality modeling: a formal Bayesian approach, Advances in Water Resources (2018), doi: 10.1016/j.advwatres.2018.04.006
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.
ACCEPTED MANUSCRIPT
Highlights
AC
CE
PT
ED
M
AN US
CR IP T
Joint analysis of input and parametric uncertainties by a rigorous Bayesian method Markov Chain Monte Carlo analysis tightly coupled with Bayesian Model Averaging An inverse modeling approach to joint inference of uncertain parameters and inputs Impact of parameter-input interaction on modeling uncertainty addressed
1
ACCEPTED MANUSCRIPT
Joint analysis of input and parametric uncertainties in watershed water quality modeling: a formal Bayesian approach
a
CR IP T
Feng Han a, b, c and Yi Zheng a, c *
School of Environmental Science and Engineering, Southern University of Science and Technology, Shenzhen 518055, Guangdong Province, PR China.
b
School of Water Resources and Hydropower Engineering, Wuhan University, Wuhan
c
AN US
430072, Hubei Province, PR China.
Guangdong Provincial Key Laboratory of Soil and Groundwater Pollution Control, Southern University of Science and Technology, Shenzhen 518055, Guangdong
ED
M
Province, PR China.
* Corresponding author:
PT
Dr. Yi Zheng, Professor
Room 506, Block 2, Wisdom Valley, Southern University of Science and Technology,
CE
Shenzhen 518055, Guangdong Province, PR China
AC
Email: [email protected] Work phone/fax: 86-755-88018030
2
ACCEPTED MANUSCRIPT
Abstract Significant Input uncertainty is a major source of error in watershed water quality (WWQ) modeling. It remains challenging to address the input uncertainty in a rigorous Bayesian framework. This study develops the Bayesian Analysis of Input and
CR IP T
Parametric Uncertainties (BAIPU), an approach for the joint analysis of input and parametric uncertainties through a tight coupling of Markov Chain Monte Carlo
(MCMC) analysis and Bayesian Model Averaging (BMA). The formal likelihood
AN US
function for this approach is derived considering a lag-1 autocorrelated,
heteroscedastic, and Skew Exponential Power (SEP) distributed error model. A series of numerical experiments were performed based on a synthetic nitrate pollution case
M
and on a real study case in the Newport Bay Watershed, California. The Soil and
ED
Water Assessment Tool (SWAT) and Differential Evolution Adaptive Metropolis (DREAM(ZS)) were used as the representative WWQ model and MCMC algorithm,
PT
respectively. The major findings include the following: 1) the BAIPU can be
CE
implemented and used to appropriately identify the uncertain parameters and
AC
characterize the predictive uncertainty; 2) the compensation effect between the input and parametric uncertainties can seriously mislead the modeling based management decisions, if the input uncertainty is not explicitly accounted for; 3) the BAIPU accounts for the interaction between the input and parametric uncertainties and therefore provides more accurate calibration and uncertainty results than a sequential analysis of the uncertainties; and 4) the BAIPU quantifies the credibility of different 3
ACCEPTED MANUSCRIPT
input assumptions on a statistical basis and can be implemented as an effective inverse modeling approach to the joint inference of parameters and inputs. Keywords: Water quality modeling; watershed; Bayesian Model Averaging; Markov
AC
CE
PT
ED
M
AN US
CR IP T
Chain Monte Carlo; input uncertainty; Bayesian calibration
4
ACCEPTED MANUSCRIPT
1. Introduction Watershed water quality (WWQ) models typically integrate distributed hydrological modeling with water quality modeling, and depict the interrelated watershed processes of water, heat, sediments and chemicals [1]. Some representative
CR IP T
WWQ models include the Soil and Water Assessment Tool (SWAT) [2, 3], Watershed Analysis Risk Management Framework (WARMF) [4, 5], Hydrological Simulation Program - Fortran (HSPF) [6], and Annualized Agricultural Non-Point Source
AN US
Pollution model (AnnAGNPS) [7, 8]. These models have been widely used to address a variety of water and environmental issues, including agricultural non-point source (NPS) pollution [9-13], water quality response to climate change [14, 15],
M
food-energy-water nexus (FEW) [16, 17], low impact development (LID) [18], and
ED
many others. However, simulation results of the WWQ models often have notable uncertainty [19, 20]. First, owing to current knowledge gaps, the existing models
PT
often represent complex water quality processes (e.g., soil erosion, heat exchange,
CE
pollutant transport and chemical reactions) in a highly simplified way, which leads to
AC
significant model structural errors [20-22]. Second, while the quantity and quality of the hydrological forcing data (e.g., precipitation and temperatures) have been greatly improved due to the development of remote sensing [23], water quality drivers (e.g., point source (PS) and NPS loadings) still lack good data with the required temporal and spatial resolutions [24], which represents a significant input uncertainty in the WWQ models. Third, WWQ models involve a larger number of model parameters, 5
ACCEPTED MANUSCRIPT
and the parametric uncertainty would be substantial when observational data for model calibration are scarce and/or inaccurate. The multiple sources of uncertainty significantly hinder the application of WWQ models in scientific research and management practices, and effective approaches to uncertainty analysis are highly
CR IP T
desired for the models.
Bayesian calibration using Markov chain Monte Carlo (MCMC) sampling [25, 26] is a formal approach for model calibration and uncertainty assessment and has been
AN US
applied widely in recent years [27-32]. In the approach, model residuals (i.e., the
difference between observations and corresponding model outputs) are described by a statistical residual error model, based on which the likelihood function can be derived
M
[33]. When the prior distribution is provided, the posterior parameter distribution can
ED
be inferred by adopting the Bayes’ rule. The approach has explicit assumptions and a rigorous mathematical basis, but its effectiveness is dependent on many issues. One
PT
issue is the appropriateness of the error model. For example, in hydrological modeling,
CE
model residuals are often correlated, heteroscedastic and not normally distributed [34].
AC
There are a series of studies which discuss the residual model [28, 30, 35-39]. Another important issue is the consideration of input uncertainty and model structural errors. A classic Bayesian calibration approach directly addresses the parametric uncertainty and leaves out the other two uncertainty sources, which may bias the inference of the posterior distributions. Our previous study shows that such degradation can be alleviated by including additional model responses that are less impacted by the input 6
ACCEPTED MANUSCRIPT
and model structural errors in a multiple-response Bayesian calibration (MRBC) framework [24]. Another solution is to jointly assess the input, structural and parametric uncertainties, although this remains challenging [40]. Hierarchical Bayesian approaches, such as Bayesian Total Error Analysis
CR IP T
(BATEA) [41], Integrated Bayesian Uncertainty Estimator (IBUNE) [42], Integrated Parameter Estimation and Uncertainty Analysis Tool (IPEAT) [43], and some others [44, 45] have been developed to quantify the model input uncertainty along with
AN US
model calibration. In the existing approaches, input errors are usually represented by latent variables (e.g., rainfall multipliers) with prescribed hyperdistributions. The hierarchical Bayesian inference addresses all uncertain model parameters and latent
M
variables, as well as all hyperparameters [46]. Although these approaches can separate
ED
the input and parametric uncertainties, their hierarchical inference framework can lead to high dimensional posterior distributions, and sampling of the posterior distributions
PT
can be computationally prohibitive, especially for complex WWQ models. More
CE
importantly, in WWQ modeling, it is very difficult to represent the input uncertainty
AC
of pollutant loadings in terms of latent variables. Bayesian Model Averaging (BMA) has been increasingly used to account for
model structural uncertainty [42, 47-53]. BMA derives an aggregated simulation by averaging the simulations of multiple models based on their posterior probabilities (also referred to as the model weights) [54, 55]. BMA has also been used for model ranking and selection, where the model weights should balance goodness of fit with 7
ACCEPTED MANUSCRIPT
complexity following the law of parsimony [56-58]. In hydrological and water quality modeling, it is a common situation for multiple plausible datasets or assumptions to be used to prepare a specific model input. For example, precipitation data can be simulated by a regional climate model, interpolated from observations at
CR IP T
meteorological stations, or estimated based on radar data. Land cover data can be interpreted from different types of satellite data, and different agencies may have
different estimates of pollutant loadings. In all these cases, the input uncertainty can
AN US
hardly be represented as continuous probability distributions, and needs to be
addressed in a discrete form. BMA appears to be an appropriate way to address such input uncertainty. For example, Strauch et al. [59] established four SWAT models
M
with four different sets of precipitation data, and then performed BMA to account for
ED
the precipitation uncertainty in streamflow modeling. Some studies have attempted the coupling of BMA with Bayesian calibration in a sequential manner, that is,
PT
performing Bayesian calibration for individual models first and then achieving an
CE
aggregated simulation using BMA. For example, Zeng et al. [52] sequentially
AC
performed MCMC and BMA analyses to assess groundwater conceptualization uncertainty. However, this sequential approach ignores the interaction between the parametric uncertainty and conceptual uncertainty, which may lead to inaccurate results under certain circumstances. This study aims to develop an approach to investigating parametric and input uncertainties simultaneously through a tight coupling of MCMC and BMA. The 8
ACCEPTED MANUSCRIPT
approach is specifically designed for input uncertainty in a discrete form, as exemplified above. A previously developed SWAT model for the nitrate pollution in the Newport Bay Watershed (Southern California, U.S.) was used as the study case, and a state-of-art MCMC algorithm, the Differential Evolution Adaptive Metropolis
CR IP T
(DREAM(ZS)) [44, 60, 61], was employed. A series of numerical experiments were
designed and implemented to answer the following key questions: 1) whether the new approach can effectively address the discrete form input uncertainty in MCMC-based
AN US
Bayesian calibration; 2) how the discrete form input uncertainty in pollutant loadings would impact the Bayesian calibration results and potential water quality management decisions based on the results; and 3) what information the joint analysis of input and
M
parametric uncertainties can provide to measure the relative importance of the
ED
competing model input assumptions.
2. Coupling MCMC and BMA
PT
2.1. MCMC uncertainty analysis
CE
Markov Chain Monte Carlo (MCMC) uncertainty analysis represents a classic
AC
Bayesian calibration approach. In this approach, the residual errors, i.e., ε Z F , are described by a statistical error model, where F and Z are the model outputs and the observations, respectively. This lumped error accounts for all sources of uncertainty including model input, model structure, model parameters and observations for model calibration [36]. With the mathematical form of the error model specified, the likelihood function of model parameters can be derived, and the posterior 9
ACCEPTED MANUSCRIPT
distributions of the parameters can be obtained based on Bayes’ rule, as
p(θ, φ | X, Z) p(θ) p(φ) L(θ, φ, X, Z)
(1)
where θ represents model parameters, φ denotes the parameters of the error model, X denotes the model input, p(θ) and p(φ) are the prior distributions, and
CR IP T
L(θ, φ, X, Z) is the likelihood function. In an MCMC algorithm, samples are
obtained by constructing Markov chains with desired posterior distributions as their
stationary distributions [25, 26]. The samples obtained by the Markov chains are used
AN US
to derive the posterior distribution of uncertain parameters, p(θ, φ | X, Z) , as well as
the predictive distribution, p(Y | θ, φ, Z) , where Y denotes a vector of quantities to be predicted. Explanation of the approach to derive the predictive distribution can be
ED
hydrological modeling [27-32].
M
found elsewhere [62]. This paradigm has been widely used in uncertainty analysis for
This study considers an autocorrelated, heteroscedastic, and non-Gaussian error
PT
model [24, 62]. The autocorrelation of residuals is depicted by a first order
CE
autoregressive model (i.e., AR(1)), which involves a parameter, lag-1 autoregressive parameter ϕ. The heteroscedasticity is reflected by a linear equation, t b1 b2 ft ,
AC
where t is the estimated standard deviation, f t is the model output at time step t, and b1 and b2 are the two parameters. The Skew Exponential Power (SEP) distribution [36] is employed to represent two typical non-Gaussian effects, i.e., skewness and kurtosis, where the skewness parameter ξ and the kurtosis parameter β are introduced. Thus, the error model has five parameters in this case, that is, 10
ACCEPTED MANUSCRIPT
φ [ , b1 , b2 , , ] . The log-likelihood function can be formulated as logL(θ, φ, Z) J log
2
where J is the data length of Z; a , j
J
J
j 1
j 1
log j c a , j 1
sign ( a j )
2/(1 )
(2)
( a j ) , in which a j stands
CR IP T
for the independent and identically distributed random errors with a zero mean and unit standard deviation; and , , and c are functions of ξ and β. More details about the error model and likelihood function can be found elsewhere [24, 36, 38]. Note
AN US
that this highly complicated likelihood function may lead to an undesired
computational burden on the uncertainty analysis. To simplify the error model, as well as the corresponding likelihood function, a trail-and-error approach can be adopted to
M
examine the fulfillment of the statistical assumptions a posteriori. One can gradually increase the complexity of the error model, as was done in previous studies [24, 62].
ED
Differential Evolution Adaptive Metropolis (DREAM(ZS)) [44, 60, 61] is
PT
employed in this study as the MCMC algorithm. It is a population evolution MCMC algorithm which runs multiple Markov chains in parallel to facilitate the global
CE
exploration of the parameter space. It automatically tunes the scale and orientation of
AC
the proposal distribution. In DREAM(ZS), a candidate point can be generated by using a differential evolution with a randomized subspace sampling strategy, which can be described as follows: θ*,i θi (1d ed ) ( , d ) u r1 ( j ) u r2 ( n ) є d n 1 j 1
where θ*,i is the candidate point, θi is the current point of the i-th Markov chain, ed 11
(3)
ACCEPTED MANUSCRIPT
and ϵd are two d-dimensional random vectors, the coefficient γ controls the jumping size, and {ur1 ( j ) , ur2 ( n) | j, n 1,..., } represents δ pairs of previously sampled points ( r1 ( j ) r2 (n) ). In high dimension problems, it is often not optimal to sample all d dimensions simultaneously. In DREAM(ZS), only d’ dimensions of θ*,i, which are
CR IP T
randomly selected with a ―crossover probability‖, are updated in each sampling step, and the other dimensions remain the same as those of θi. To increase diversity, the
snooker update is included, which is an additional way to generate candidate points.
AN US
The way to propose a candidate point is described as follows:
θ*,i θi s θi u D(u R1 , u R 2 )
(4)
where γs is a random variable, u, uR1 and uR2 are three previously sampled points, and
M
D is a function of uR1 and uR2. Both D and γs are used to determine the distance of the
ED
jumping. More details of the snooker update can be found in ter Braak and Vrugt [63]. In DREAM(ZS), both differential evolution and snooker update are used to generate
PT
the candidate points. In this study, we updated the Markov chains in a mix of 10%
CE
snooker update and 90% differential evolution. The R statistics proposed by Gelman
AC
and Rubin [64] are used to assess the convergence of the Markov chains. The Markov chains are deemed to achieve the convergence when the R values of all dimensions (or parameters) are below 1.2, computed based on the last 50% points of the Markov chains.
2.2. Loose coupling: the sequential approach As mentioned in Section 1, MCMC-based Bayesian calibration can be coupled 12
ACCEPTED MANUSCRIPT
with BMA in a sequential manner to account for input uncertainty. Hereinafter, models with the same model structure but different inputs are considered as different models. When K plausible input datasets, X1, …, XK, are collected or assumed for a same model structure, K different models, M1, …, MK, can be established, respectively,
CR IP T
which represents the input uncertainty in a discrete form. MCMC uncertainty analysis can first be performed for each model to assess the parametric uncertainty. Next, the predictive distributions inferred based on the K models (i.e., p(Y | M k , Z) ) are
AN US
averaged into a total predictive distribution, representing both parametric and input
uncertainties. This two-step uncertainty analysis is referred to as the sequential BMA approach in this paper. Its mathematical details are provided below.
M
BMA provides a data-driven approach to statistically integrate the predictions of , yt ,
, yT denote the vector of quantities
ED
multiple models [54, 55]. Let Y y1 ,
to be predicted (i.e., the true values), and let Z denote the observations. In BMA, the
PT
predictive distribution of yt , based on the k-th model M k , can be deemed as the
CE
posterior distribution of yt conditional on Z, that is, p( yt | M k , Z) . Thus, the
AC
posterior distribution of yt conditional on the entire model ensemble can be written as
K
p( yt | Z) p( yt | M k , Z) p( M k | Z)
(5)
k 1
where p( M k | Z) stands for the posterior probability of M k and p( M k | Z) can be interpreted as the likelihood of M k being ―correct‖, given Z. Because the posterior probabilities of all the models have a sum of one, p( M k | Z) can also be 13
ACCEPTED MANUSCRIPT
interpreted as the weight of M k . Hereinafter, p( M k | Z) is denoted as wk and
w [w1 ,..., wK ] . The mean and variance of BMA predictions can be computed as K
E[ yt Z] wk E[ yt M k , Z]
(6)
k 1
K
K
Var[ yt Z] wk E[ yt M k , Z] E[ yt Z] wk Var[ yt M k , Z] 2
(7)
k 1
CR IP T
k 1
The predictive distribution p( yt | M k , Z) can be computed as
p( yt | M k , Z) p( yt | M k , θk , φk ) p(θk , φk | M k , Z)dθk dφk k
(8)
AN US
where θ k denotes the model parameters of M k ; φ k denotes the error parameters
associated with the residual error model, θk , φk k ; p( yt | M k , θk , φk ) represents the predictive distribution of yt given the specific values of θ k and φ k ; and
M
p(θk , φk | M k , Z) represents the posterior distribution of θ k and φ k . The
ED
calculation of the integral in Equation (8) is computationally intensive. Alternatively, the integral can be approximated using the posterior samples generated in the MCMC
CE
PT
analysis, as [26]
pˆ ( yt | M k , Z)
1 Nk
Nk
p( y n 1
t
| M k , θ(kn ) , φ(kn ) )
(9)
AC
where θ(kn ) and φ (kn ) denote the n-th sample of model parameters and error parameters, respectively, jointly generated by the Markov chain, and N k is the sample size. Assuming independence among the residual errors, the log-likelihood function of the model weights can be derived as
14
ACCEPTED MANUSCRIPT
K 1 logL(w, Z) log wk j 1 k 1 N k J
(n) (n) p ( z | M , θ , φ ) j k k k n 1 Nk
(10)
Note that the assumption of independence may not be valid. For example, when lag-1 autocorrelation of the residual errors exists, z j should be conditional on z j 1 . In
CR IP T
such cases, p( z j | M k , θ(kn ) , φ(kn) ) in Equation (10) should be replaced by p( z j | M k , θ(kn ) , φ(kn ) , z j 1 ) . Raftery et al. [47] suggested to estimate the weights using
the maximum likelihood estimation (MLE). For complex likelihood functions, the
AN US
expectation-maximization (EM) algorithm [65] can be employed to find a numerical solution for the MLE problem with a small computational cost. However, the EM method only leads to a set of deterministic weight estimates, and the results can be
M
unstable when the weights are sensitive to the observations. A more reliable way is to perform MCMC simulation [66], in which the means of the sampled weights can be
ED
used as the posterior probabilities of the respective models. There are also other
PT
options to determine the model weights such as the Bayesian model evidence method [58]. In this study, we adopted DREAM(ZS) to sample model weights in the sequential
CE
BMA approach. To guarantee a sum of one for the weights of a candidate point (i.e.,
AC
θ*,i), we modified the original operation of generating candidate points to the following steps: Step 1: Generate a jumping vector ∆ by Equation (3) or (4), with at least two
elements of ∆ being non-zero. Step 2: Normalize ∆ by subtracting the mean of the updated elements, and therefore the sum of the normalized ∆ is equal to zero. 15
ACCEPTED MANUSCRIPT
Step 3: Generate the candidate point θ*,i = θi + ∆. The sum of θ*,i is one, as the sum of θi is one. Step 4: If θ*,i is outside of the sampling domain, calculate a new point (θ**,i) using the ―returning‖ operation introduced below.
CR IP T
As θ*,i may be outside of sampling domain, a ―returning‖ operation was proposed, as illustrated in Figure 1. When the proposed jump from θi to θ*,i hits the boundary (i.e., θb is reached), the jump turns back along the original path and ends at θ**,i. In
AN US
this case, the total jumping distance, in its absolute value, remains the same. Because the sum of ∆ is zero, the sum of θ**,i remains one. Similarly, a jump from θ**,i with the same jumping vector ∆ will end at θi. It suggests that, with the ―returning‖ operation,
M
the conditional probability of a jump from θi to θ**,i is equal to that of the reverse
ED
jump (i.e., from θ**,i to θi). Thus, the detailed balance of DREAM(ZS) can be guaranteed by accepting the candidate point θ**,i (or θ*,i), based on the Metropolis
PT
acceptance probability [61].
CE
2.3. Tight coupling: BAIPU
AC
This study develops a new approach to jointly analyze the input and parametric
uncertainties by tightly coupling BMA with MCMC, which is named the Bayesian Analysis of Input and Parametric Uncertainties (BAIPU). Figure 2 illustrates the framework of BAIPU. It begins with the preparation of K sets of plausible input data, which result in K models, as defined in Section 2.2. Thus, given the same model structure, one realization of model parameters θ would lead to K sets of model outputs, 16
ACCEPTED MANUSCRIPT
F1, …, FK. In BAIPU, a shared residual error model with parameters φ is also assumed for the model ensemble. Hence, the multiple models share the same set of uncertain parameters (i.e., θ plus φ). The predictive distributions of yt derived based on different models can be further averaged by BMA, as K
,K
, θ, φ) wk (θ, φ) p( yt | M k , θ, φ) k 1
CR IP T
p( yt | M1,
(11)
where wk (θ, φ) represents the weight of the k-th model given specific values of θ and
AN US
φ, and p( yt | M k , θ, φ) is the predictive distribution of yt based on the k-th model
(i.e., the k-th input dataset). In Equation (11), the dependency of the model weights on θ and φ reflects the interaction between the input and parametric uncertainties, which is
M
ignored in the sequential BMA approach. Assuming independence among the residual errors, the log-likelihood function of the model parameters and error parameters can
ED
be defined as
(12)
PT
J K logL(θ, φ, w, Ζ) log wk (θ, φ) p( z j | M k , θ, φ) j 1 k 1
CE
If lag-1 autocorrelation among the residual errors exists, the probability of z j in
AC
Equation (12) should also be conditional on z j 1 , that is, p( z j | M k , θ, φ, z j 1 ) . The MCMC analysis is then performed to derive the posterior distributions of model parameters and predictive variables. As discussed in Section 2.1, there are different strategies to infer the model weights. The framework of BAIPU accommodates two (Figure 2). One strategy is to estimate the weights by the EM algorithm. When a realization of θ and φ is generated 17
ACCEPTED MANUSCRIPT
in one sampling step, the corresponding model weights wk (θ, φ) (k=1, ..., K) are then determined by maximizing this likelihood function using the EM optimization. The other strategy is to sample the model weights jointly with the model parameters using MCMC. To implement this strategy, the ―returning‖ operation introduced in
CR IP T
Section 2.2 needs to be incorporated in DREAM(ZS), such that the model weights in each candidate point add up to one. Because the MCMC sampling is more reliable
than the EM optimization, all the BAIPU results presented in this paper were based on
AN US
the MCMC sampling. For comparison, we also performed all the BAIPU analyses
using the EM approach. Section 4.1 briefly discusses the comparison, and readers are referred to the Auxiliary Materials for more details.
M
One difference between the BAIPU and sequential BMA approach is that the
ED
BAIPU derives a common posterior parameter distribution based on the entire set of plausible model inputs, while the sequential approach leads to multiple posterior
PT
parameter distributions with respect to different model inputs. Note that, the multiple
CE
posterior parameter distributions could be further averaged into a single posterior
AC
distribution by resampling the multiple distributions according to their posterior probabilities (i.e., the model weights). However, although such a posterior distribution is relevant to the parameter calibration, it cannot be directly assigned to each model for making predictions, which is an advantage of BAIPU. In addition, in the BAIPU, because a common posterior parameter distribution can be obtained, the model weights solely reflect the model input uncertainty. However, in the sequential BMA 18
ACCEPTED MANUSCRIPT
approach, the input and the parametric uncertainties can compensate for each other in the Bayesian calibration; therefore, the weights are impacted by the parametric uncertainty. The error model introduced in Section 2.1 is considered for the model ensemble,
CR IP T
which involves five error parameters, that is, φ [ , b1 , b2 , , ] . Let kt yt f kt
denote the residual error of the k-th model (i.e., input X k ) at time step t, where f kt is the corresponding model output. As mentioned before, the standard deviation of kt
AN US
is estimated by kt b1 b2 f kt , and hence the standardized residuals can be
calculated as kt kt kt . The AR(1) model is applied to the standardized residuals to account for the lag-1 autocorrelation, that is, kt k ,t 1 okt . It is a standard
1 2 should be a unit. Similarly, we assume that akt follows the SEP
ED
of akt okt
M
homoscedastic AR(1) process and variance of okt equals 1 2 . Thus, the variance
distribution. Hence, given the specific values of θ and φ, the predictive probability of
PT
yt based on the k-th model can be computed as
CE
p( yt | M k , θ, φ, yt 1 ) sign ( akt )
1 2 kt
exp c a ,kt 1
2 /(1 )
(13)
( akt ) , and , , and c are functions of
AC
where a ,kt
2
1
and [36]. According to Equation (12), the log-likelihood function for this
residual error model can be derived as, K J 2 1 logL(θ, φ, w, Ζ) log wk (θ, φ) exp c a ,kj 1 k 1 j 1 1 2 kj
2/ (1 )
(14)
This likelihood function is quite general and flexible. It relaxes the commonly used 19
ACCEPTED MANUSCRIPT
independent, homoscedastic and Gaussian assumption about the residual errors and hence is more realistic for environmental modeling. As was revealed by Han and Zheng [24], when model input errors are very significant, the inferred posterior distributions of the model parameters would be prone to notable bias. This likelihood
CR IP T
function, which is based on multiple input datasets rather than one single set of input data, is expected to be more accurate and reliable. It is worth re-emphasizing that
BAIPU was specifically designed to address discrete form input uncertainty. If the
AN US
input uncertainty could be appropriately characterized as continuous probability distributions, the BMA step in BAIPU may not be an efficient solution.
In BAIPU, the sampler, DREAM(ZS), would execute K times more model runs
M
than in a traditional MCMC uncertainty analysis, if the chains have the same length.
ED
In this study, the computational cost was largely reduced by parallelizing the sampling process, that is, making DREAM(ZS) execute Nchain*K model runs (Nchain denotes the
PT
number of Markov chains) simultaneously in parallel. Similar parallel calculation was
CE
adopted for the sequential BMA as well.
AC
3. Data and Methods 3.1. Study area The Newport Bay Watershed (NBW) is in Orange County, southern California. It
is a highly urbanized watershed with an area of approximately 400 km2 (Figure 3). As of 2001, approximately 70% of the watershed was residential, commercial and industrial areas, and agricultural and orchard areas accounted for approximately 8%. 20
ACCEPTED MANUSCRIPT
The NBW has a typical Mediterranean climate with warm, dry summers and mild, wet winters. The annual rainfall is approximately 330 mm, and occurs mostly between November and April. Nutrient pollution was a significant water quality issue in this area in the 1990s. Excessive loadings of nutrients were delivered into the Newport
CR IP T
Bay through the San Diego Creek (see Figure 3), which had caused the eutrophication of the bay and created the recreational and aesthetic nuisance [67]. Major sources of
nutrients included fertilizer application (non-point source) and commercial nurseries
AN US
(point sources). The NBW and San Diego Creek were included in the 1998 California 303(d) list of impaired waters for Total Maximum Daily Loads (TMDL) actions, with nutrients as a major pollutant stressor [68]. More details about this watershed can be
ED
3.2. Nitrate pollution modeling
M
found in our previous studies [62, 69].
SWAT is a representative WWQ model which is semi-distributed and runs at a
PT
daily time step [3]. A SWAT model of the nitrate pollution in the NBW has been
CE
previously developed [24, 62]. In brief, the watershed was first delineated into 11
AC
sub-basins (see Figure 3), with a total of 58 hydrologic response units (HRUs). The 90-m Digital Elevation Model (DEM) of the study area was acquired from the U.S. Geological Survey (USGS) to delineate the modeling domain. The 30-m land use raster from the National Land Cover Database 2001 (NLCD 2001) and the soil type map from the U.S. State Soil Geographic (STATSGO) database were referred to in defining the HRUs. Meteorological data from 1998 to 2004 were obtained from U. S. 21
ACCEPTED MANUSCRIPT
National Climatic Data Center (NCDC) and used to drive the model simulations. The model was run from July 1998 to June 2004, with the first two years being the warm-up period and the last four years being the calibration period. Observations of daily flow and weekly nitrate concentration at the watershed outlet (Figure 3a) from
Control Board reports (available at
CR IP T
July 2000 to June 2004 were collected from the Santa Ana Regional Water Quality
http://prg.ocpublicworks.com/DocmgmtInternet/Search.aspx, last accessed on
AN US
October 8, 2017) and were used for the model calibration.
In this study, nitrate was the water quality parameter of concern, and nitrogen source loadings were considered as the uncertain model inputs. In the NBW, fertilizer
M
application on urban lawns and effluents from commercial nurseries were the two
ED
major nitrogen sources [67]. In the SWAT model, Bermuda was selected as the representative lawn grass and the nurseries were conceptualized into two point
PT
sources, A and B (see Figure 3a). As detailed information of the timing of fertilizer
CE
application is not available, the heat unit scheduling in SWAT was chosen, which
AC
allows the model to schedule operations based on temperature. To simulate the stormwater runoff and nitrate pollution in urban areas, SWAT offers two options. One is to use a set of linear empirical equations developed by the USGS (referred to as the USGS option), and the other is to simulate the buildup and washoff processes (referred to as the buildup-washoff option). Although the buildup-washoff option is conceptually appealing, local data to parameterize the equations is often lacking. 22
ACCEPTED MANUSCRIPT
Hence, the USGS option was chosen in this study. In our previous study [24], the nitrate loadings at sources A and B were estimated to be 61 kg N/day and 31 kg N/day, respectively, and the fertilizer usage was assumed to be 240 kg N/ha/year. The PS loadings were assumed to be time-invariant because
CR IP T
no detailed temporal information was available. Inevitably, significant input
uncertainty exists in the source loadings, which would degrade the calibration of the
uncertain model parameters. Thus, this study considers multiple plausible PS and NPS
AN US
loadings instead, which will be further discussed in Section 3.4. 3.3. Multiple-response Bayesian calibration
To alleviate the impact of input and model structural errors, the streamflow and
M
nitrate concentration at the watershed outlet are calibrated simultaneously by applying
ED
the MRBC approach developed in our previous study [24]. In the context of the MRBC, different error models are required for different responses, based on which
PT
individual likelihood functions can be established. In this study, residual errors of
CE
streamflow and nitrate concentration are assumed to be independent. Hence, the
AC
combined likelihood function is simply the product of their likelihood functions [24, 70]:
Lmultiple LH (θ, φH , ZH ) LN (θ, φ N , Z N )
(15)
where LH and LN are the likelihood functions of streamflow and nitrate concentration, respectively. In this study, the likelihood functions of flow and nitrate responses are determined by Equation (2) and Equation (14) in the sequential BMA 23
ACCEPTED MANUSCRIPT
approach and BAIPU, respectively. 3.4. Numerical experiments In the real-situation case, five loading scenarios were assumed to account for the input uncertainty, which leads to five real-situation models (RMs) (see Table 1).
CR IP T
Among the five models, RM2 is the basic model that was built in our previous study
[24] and represents an average loading situation. RM1 has 20% lower PS loadings and 10% lower fertilizer usage, compared to RM2, while RM3 has 20% higher PS
AN US
loadings and 10% higher fertilizer usage. RM4 and RM5 are two models with
temporal variability in the PS loadings, while RM4 has a smaller seasonal fluctuation (see Figure S1 in the Auxiliary Materials). Note that the overall loadings in RM4 and
M
RM5 are the same as those in RM2. This ensemble of models effectively embodies a
ED
common situation encountered in WWQ modeling, that is, detailed and accurate loading data are often lacking, and the modeler may consider different assumptions.
PT
The selection of random parameters in the Bayesian calibration was the same as
CE
in our previous study [24]. In brief, 32 SWAT parameters that are physically relevant
AC
for flow and/or nitrate simulations were identified first. Next, 12 most important parameters (Table 2) were further determined after an initial sensitivity analysis on the 32 parameters using Morris screening [71]. Each of the 12 parameters ranks top 8 among the 32 initial parameters at least for one of the two responses (i.e., flow and nitrate). Full details of the selection procedure can be found in our previous study [24]. As most of the parameters have spatially distributed values, three strategies [30] were 24
ACCEPTED MANUSCRIPT
considered to vary them: varying the parameter value directly (Type I), adding a deviation value to the distributed parameters (Type II), and applying a multiplier to the distributed parameters (Type III). For parameters of Type II and III, the deviations and multipliers other than the original distributed parameters are treated as uncertain
CR IP T
model parameters in subsequent calibration approaches. All the uncertain parameters are assumed to be uniformly distributed, and their ranges are provided in Table 2.
Initialization of other non-random parameters is also the same as in the previous study
AN US
[24].
To investigate the applicability and performance of BAIPU, a synthetic modeling case was also constructed, in which a true model (denoted as TM) with true model
M
parameters, true model inputs and true model structure was hypothesized (Table 1).
ED
The true values of the 12 uncertain parameters (the last column in Table 2) were randomly picked, which can yield realistic simulation results. True nitrate loadings at
PT
sources A and B were assumed to be time-variant at the levels of 132 and 68 kg N/day,
CE
respectively (see Figure S2 in the Auxiliary Materials). The fertilizer usage in TM was
AC
assumed to be 240 kg N/ha/year and scheduled on specific dates. The buildup-washoff method was selected as the true model structure for the simulation of stormwater pollution from urban areas. Other model settings were the same as that in the real-situation case. True model responses of flow and nitrate concentration at the watershed outlet were achieved by running the true model with the above settings. In addition to the perfect model (i.e., TM), five imperfect models (IMs) with different 25
ACCEPTED MANUSCRIPT
source loadings were also created, denoted as IM1, IM2, IM3, IM4 and IM5 (see Table 1). The PS loadings were set to be time-invariant, which reflects the reality that the temporal information of PS loadings is often unavailable. Among the models, IM3 represents the average loading level (i.e., 150 kg/day and 75 kg/day at sources A and
CR IP T
B, respectively), and IM1, IM2, IM4 and IM5 represent the 66.7%, 83.3%, 116.7% and 133.3% of the average loading level, respectively. Similarly, different fertilizer
application rates were hypothesized for the five imperfect models (Table 1). The five
AN US
imperfect models collectively represent the hypothetical input uncertainty. To better mimic the reality, the model structural error was also introduced in the IMs. The USGS option was used instead of the buildup-washoff option. Daily flow
M
observations and weekly nitrate observations at the outlet were synthesized by
ED
corrupting the true model responses with hypothetic observational errors (see Figure S3 in the Auxiliary Materials). The hypothetic errors were generated from an error
PT
model with b1H 0.2 , b2H 0.05 , H 0.5 , H 3 , H 0.2 , b1N 0.2 , b2N 0.3 ,
CE
N 0.5 , N 3 , and N 0 . This procedure represents a common strategy to
AC
synthesize observational data [19, 36, 37]. As in the real-situation case, the artificial observations were generated for the period from July 2000 to June 2004, and all the models (i.e., the TM and IMs) were run from July 1998 to June 2004. In both the synthetic and real-situation cases, the BAIPU and the sequential BMA approach were carried out. In the BAIPU, the model weights were jointly sampled with the model and error parameters (i.e., the joint sampling method), and in the 26
ACCEPTED MANUSCRIPT
sequential BMA approach, the model weights were sampled by DREAM(ZS). Hereinafter, BAIPU_IMs denotes the BAIPU analysis on the five imperfect models, and MCMC_IM1, MCMC_IM2, MCMC_IM3, MCMC_IM4 and MCMC_IM5 denote the MCMC uncertainty analyses for the five individual models, respectively.
CR IP T
These individual MCMC analyses constitute the first step in the sequential BMA approach. The same terminology is applied to the real-situation case. The
computational efficiency, posterior distributions of the uncertain SWAT parameters
AN US
and the uncertainty results of the two approaches were evaluated and compared. It is
worth pointing out that each model (e.g., IM1) in the sequential BMA approach has its individual error model, while the multiple models in the BAIPU share the same error
M
model. In determining the proper error model, the complexity of the error models was
ED
gradually increased, from simple independent, homoscedastic and Gaussian type to the most complicated type, until the posterior checks confirmed that the residual
PT
errors are consistent with the error model assumptions. The error models for different
CE
calibration procedures are summarized in Table S1 in the Auxiliary Materials.
AC
Posterior justification of these error models is illustrated by Figures S4 to S8 in the Auxiliary Materials. Eventually, nine error parameters were included in the uncertainty analysis, i.e., φ [ H , b1H , b2H , H , H , b1N , b2N , N , N ] . In this study, the error parameters were jointly inferred with the uncertain SWAT parameters (i.e., θ) in the MCMC sampling. To further compare the sequential BMA approach and the BAIPU, three future 27
ACCEPTED MANUSCRIPT
loading scenarios with increased, steady and decreased loadings were proposed for the period from July 2004 to June 2008 (Table 3). The calibrated models were run to predict the nitrate concentration with the proposed future loading scenarios. To investigate how the calibration procedure would impact management decisions, two
CR IP T
management-relevant variables, average nitrate concentration (AVC) and frequency of exceeding (FOE) the water quality target (6 mg N/L), were evaluated based on the predictions. In the synthetic case, the true nitrate concentrations under different
AN US
loading scenarios were obtained by running the true model (i.e., TM), whereby the true values of AVC and FOE can be achieved. In the real-situation case, the true model was unknown, and hence the true AVCs and FOEs were not obtainable.
M
4. Results and Discussion
ED
4.1. Efficiency of DREAM(ZS) in BAIPU
Twelve Markov chains were run simultaneously in DREAM(ZS), and 12×K CPUs
PT
(2.67 GHz) were utilized in parallel to reduce the computing time, where K is the size
CE
of model ensemble. With the parallel computing, the cost of implementing
AC
DREAM(ZS) is affordable. For example, evaluating 200,000 samples in the BAIPU typically takes one day for the modeling cases. The sequential BMA approach has two stages of DREAM(ZS). The second stage is to sample the model weights, and the computational cost is only a few minutes because this stage involves no SWAT runs. To ensure that the posterior distributions and uncertainty results were reliable and representative, DREAM(ZS) explored 400,000 and 300,000 samples for the synthetic 28
ACCEPTED MANUSCRIPT
and real-situation cases, respectively. Using the fixed sample numbers can rule out the effect of computational cost when interpreting the analysis results of different models. The numbers are large enough to ensure the convergence of Markov chains in most of the experiments.
CR IP T
Table 4 shows the number of samples for Markov chains to achieve convergence (NC) and the acceptance rate (AR) of DREAM(ZS) in different numerical experiments. As indicated by the NC values, the convergence speed of Markov chains in BAIPU
AN US
was comparable to that in the MCMC uncertainty analyses in the real-situation case, but was moderately slower in the synthetic case. One obvious reason for the speed
decrease is that BAIPU leads to a higher dimensional posterior distribution. AR is the
M
probability that one candidate point is accepted in stationarity. It can be estimated as
ED
T'/T, where T and T' denote the numbers of all and the accepted candidate points in the second halves of the Markov chains, respectively. A higher AR simply means that it is
PT
easier to find an acceptable candidate point with relatively high posterior probability.
CE
As Table 4 shows, in the two-step sequential BMA, DREAM(ZS) obtained one order of
AC
magnitude higher AR values in estimating the model weights than in deriving the posterior parameter distributions. This is because the likelihood function of the model and error parameters (Equation (2)) is much more complex than that of the model weights (Equation (10)). Another interesting finding is that the two BAIPU analyses yielded several-fold higher AR than the respective MCMC uncertainty analyses of individual models. This is probably because the multi-model ensemble allows a 29
ACCEPTED MANUSCRIPT
broader range of predictions than a single model, making it easier to find an acceptable candidate point. Overall, the results in Table 4 indicate that the sampling efficiency of DREAM(ZS) in BAIPU was not significantly reduced by the integration with BMA.
CR IP T
As mentioned in Section 2.3, we also tested the EM strategy to infer the model weights in the BAIPU. It has been found that the EM strategy led to a faster
convergence of Markov chains (NC is 363,120 and 105,120 in the synthetic and
AN US
real-situation cases, respectively) but lower acceptance rates (AR is 2.75% and 7.14% in the synthetic and real-situation cases, respectively). The EM and joint sampling strategies yielded very similar posterior parameter distributions, uncertainty results,
M
and future predictions (refer to Figures S9 to S12 in the Auxiliary Materials). Yet the
ED
model weights are quite different with the two strategies (Figure S13 in the Auxiliary Materials). These results demonstrate that the joint sampling is more reliable for the
PT
weight estimation. However, if the convergence speed is an important concern, and
CE
accurate weight estimation is not a goal, the EM strategy can be desirable.
AC
4.2. The synthetic modeling case 4.2.1. Posterior parameter distributions Figure 4 illustrates the posterior parameter distributions derived by the BAIPU
(red lines) and the MCMC analyses (dotted lines) for the individual models in the synthetic case. The true values of the model parameters are indicated by the dashed vertical lines in the figure. All the parameter values were rescaled to [0, 1] by a linear 30
ACCEPTED MANUSCRIPT
min-max normalization. The posterior parameter distributions derived by the BAIPU (red lines) are relatively narrow and generally embrace the respective true parameter values, except for the initial depth of water in the shallow aquifer (SHALLST) and the denitrification threshold water content (SDNCO) parameter. The posterior distribution
CR IP T
of SHALLST is wide, but still embraces the true value; while the posterior
distribution of SDNCO is far off the true value. This is probably because the
uncertainty of this parameter compensated the model structure error not explicitly
AN US
accounted for in this study, and/or the interaction between the two parameters and the input assumptions in the MCMC sampling is significant. Figure 4 indicates that the model parameters have been appropriately identified by the BAIPU in general. For
M
most parameters, the posterior distributions based on the individual MCMC analyses
ED
are different from each other, as well as from the respective distribution derived by the BAIPU. The difference is most evident for IM5, the model with the most
PT
significant input error. This indicates that the calibration results can be highly
CE
dependent on the input assumption; therefore, an ensemble approach such as BAIPU
AC
is valuable for addressing the parametric uncertainty. Similar calibration results of the error parameters are presented in Figure S14 in the Auxiliary Materials. 4.2.2. Model weights Figure 5 compares the model weights based on the BAIPU and the sequential BMA approach. The weight results in the two approaches are similar. In general, the lower the input loading, the higher the model weights. Among the five models, the 31
ACCEPTED MANUSCRIPT
input scenario of IM2 is closest to the true one (see Table 1), but it only achieves the second highest weight in both approaches. IM1 receives the highest weight in both approaches, while its input scenario is the third closest to the true one. The hypothesized true PS loadings have skew distributions featured by a limited number
CR IP T
of loading pulses (see Figure S2 in the Auxiliary Materials); therefore, the mean values of the true nitrate loadings (e.g., 131.95 kg N/day at Source A) are much
higher than the respective medians (e.g., 88.58 kg N/day at Source A). This explains
AN US
why the inference favors the lowest input loading scenario IM1s, because the constant loadings in IM1 are much closer to the medians of the true PS loadings. Both BAIPU and sequential BMA enable the correlation analysis of the model
M
weights, which helps in understanding the tradeoff between different input scenarios.
ED
It was found that both approaches lead to very similar correlation results (see Table S2 in the Auxiliary Materials). BAIPU also allows a correlation analysis between the
PT
model weights and the parameters because they are jointly sampled by DREAM(ZS).
CE
As Table 5 shows, strong correlations exist in many parameter-weight pairs,
AC
especially for SHALLST and SDNCO, the two parameters that were not well identified. This implies that the significant interaction between the two parameters and the input assumptions in the MCMC sampling may be one reason for the weak identifiability of the two parameters. It also confirms the necessity of sampling model parameters and weights jointly in BAIPU.
32
ACCEPTED MANUSCRIPT
4.2.3. Uncertainty results Figure 6 displays the uncertainty results of the nitrate concentration achieved by the BAIPU and the sequential BMA approach. The medians and confidence intervals in Figure 6a and 6b are very similar. In both approaches, the seasonal variation is well
CR IP T
captured by the time series of the median concentrations. The Nash–Sutcliffe
coefficients of the median predictions by BAIPU and the sequential BMA approach
are 0.2932 and 0.2804, respectively. The observations are also adequately embraced
AN US
by the uncertainty bands. As demonstrated by Figures 6a and 6b, 48.33% and 50.72% of the observations are within the 50% confidence intervals, respectively, and 97.13% and 97.61% of the observations are within the 95% confidence intervals, respectively.
M
Figure 6c and 6d show that the means and standard deviations of the predictive
ED
distributions derived in the two approaches are very close. In the sequential BMA approach, the predictive means of individual models are first calculated, and then the
PT
ensemble predictive means are derived based on Equation (6), while in the BAIPU,
CE
the ensemble predictive means for each sampled realization of uncertainty parameters
AC
are first calculated according to Equation (6), and then the total predictive means are derived by averaging over all the sampled realizations of uncertainty parameters. The standard deviations illustrated in Figure 6d can be derived in a similar way based on Equation (7). The Quantile-Quantile (QQ) plots [46] of nitrate concentration in the synthetic case were also produced for the BAIPU, individual MCMC analyses and the 33
ACCEPTED MANUSCRIPT
sequential BMA approach (Figure 7a). The QQ plots can provide intuitive information on the consistency of predictive distribution with observations. In a QQ plot, if all the points are on the 45-degree line, the predictive distribution agrees perfectly with the observations. In Figure 7a, all the predictive distributions demonstrate adequate
CR IP T
consistency. In general, the BAIPU and the sequential approach perform better than individual MCMC analyses, because both approaches account for the input
uncertainty. Among the individual models, IM2 leads to the best consistency, because
AN US
its loading scenario is closest to the true one. It is worth pointing out that, although
the predictive distributions of different approaches are similar (Figure 6 and Figure 7a), the corresponding posterior parameter distributions significantly vary, which may
M
have a significant impact on future predictions and management decisions based on
ED
the stochastic modeling. This issue will be discussed further in Section 4.2.4. 4.2.4. Future predictions and management decisions
PT
The WWQ models are commonly used to assess hypothetical loading scenarios,
CE
based on which management decisions can be made with a strong scientific basis. To
AC
investigate the performance of the different approaches in predicting future water quality and therefore supporting decision-making, three future loading scenarios were assumed in the synthetic case, as discussed in Section 3.4 and summarized in Table 3. The future loading scenarios are free of input errors. Two management variables, average nitrate concentration (AVC) and frequency of exceeding (FOE) the water quality target (i.e., 6 mg N/L), were evaluated using the true parameter values, as well 34
ACCEPTED MANUSCRIPT
as the parameter values sampled in the Bayesian analyses. Figure 8 illustrates the cumulative frequency curves (CFC) of AVC and FOE derived by different approaches. Note that, in the BAIPU and MCMC analyses for individual models, the frequencies of AVC and FOE were calculated using the
CR IP T
parameter realizations in the second halves of the Markov chains. In the sequential BMA approach, the frequencies are weighted averages of the corresponding
frequencies derived by individual MCMC analyses. Thus, the range of corresponding
AN US
CFC depends on both the spread of the CFCs by the MCMC analyses and the model
weights. As Figure 8 shows, in general, the CFCs derived by BAIPU (solid red lines) and the sequential BMA approach (solid blue lines) are both close to the respective
M
true values. It indicates that the two approaches are relatively effective in reducing the
ED
bias caused by the input errors in the calibration stage. The CFCs by BAIPU are relatively narrower because BAIPU derived a common posterior parameter
PT
distribution. Moreover, the BAIPU demonstrates a better performance than the
CE
sequential approach in predicting the management variables, which reflects the fact
AC
that the interaction between input uncertainty and parametric uncertainty is considered. In contrast, the five individual models lead to distinct predictions, some of which are far away from the respective true values. This implies that the potential management decisions based on the Bayesian calibration of individual models would be seriously biased if the model inputs have large errors. Furthermore, the impact of the bias on potential decision-making could be dependent on the confidence level desired by the 35
ACCEPTED MANUSCRIPT
decision maker. Take FOE in Scenario 3 as an example (refer to Figure 8f). If a confidence level of 50% is considered, the value of FOE would be 0.28 according to IM5; therefore, the discrepancy with the true value is 0.14. However, given a confidence level of 95%, the discrepancy would instead be 0.07.
CR IP T
Figure 8 also shows that, IM1, which underestimates the observations in the calibration stage (i.e., MCMC_IM1), overpredicts AVC and FOE; IM5, which
overestimates the observations, now underpredicts AVC and FOE. This is because, in
AN US
the Bayesian calibration of individual models, the model parameters compensate for the effect of the input errors, such that the observations can be better matched. For example, underestimated PS loadings would lead to parameter values prone to
M
overestimate the water quality response (i.e., the case of IM1). Thus, when the model
ED
inputs are free of error, the compensation in the calibration stage causes an ―overreaction‖ of the model at the prediction stage. Similar overreaction behaviors,
PT
although not as significant, have also been observed for other individual models.
CE
Interestingly, IM3 has the best performance for predicting the two management
AC
variables, instead of IM2 (with the smallest input errors), which implies the compensation effect is complex and nonlinear. All the above results suggest that a Bayesian approach accounting for multiple
plausible input assumptions can provide a valuable mechanism for alleviating the impact of input uncertainty and therefore enhancing the reliability of decision-making in the context of real-world management. In addition, with the interaction between the 36
ACCEPTED MANUSCRIPT
input and parametric uncertainties systematically accounted for, the BAIPU has an improved performance in predicting management variables under the future loading scenarios compared to the straightforward sequential BMA approach. It is also worth emphasizing that one needs to carefully make input assumptions based on his/her best
CR IP T
knowledge when implementing the two approaches. This may require extra work such as data collection and compiling, field survey, literature review, etc. In theory, it is not a problem for BAIPU to accommodate many inappropriate input assumptions, but in
4.3. The real-situation modeling case
AN US
practice it may significantly increase the computational burden.
The performances of BAIPU and sequential BMA were also examined using the
M
real-situation modeling case (Table 1). Figure 9 shows the posterior parameter
ED
distributions derived by the BAIPU (red lines) and the MCMC analyses for individual models (dotted lines). Compared to the synthetic case (Figure 4), the difference
PT
between the BAIPU and the individual MCMC analyses is relatively small, which
CE
implies that the influence of input uncertainty on the parameter calibration is not as
AC
significant as that in the synthetic case. SHALLST and SDNCO are the two parameters with the most significant difference in the posterior distributions, which is also observed in Figure 4. In addition, the posterior distributions of the error parameters are presented in Figure S15 in the Auxiliary Materials. Figure 10 illustrates the model weights inferred by the BAIPU and the sequential BMA approach in this case. The overall weight patterns derived by the two 37
ACCEPTED MANUSCRIPT
approaches are similar. As introduced in Section 3.2, RM2 represents an average estimation of the nitrogen loadings based on the available data, while RM1 and RM3 represent a lower loading scenario and a higher loading scenario, respectively. Among the first three models (Figures 10a to 10c), RM1 receives the largest weights while
CR IP T
RM2 and RM3 receive the second and third largest weights in both approaches. It
suggests that the original loadings in RM2 has a large probability of overestimation.
On the other hand, in RM4 and RM5, the PS loadings are assumed to have small and
AN US
large seasonal fluctuations, respectively. The potential seasonal fluctuation could be
due to the temporal variability of tailwater discharges from the commercial nurseries. As shown in Figures 10d and 10e, both RM4 and RM5 were assigned significant
M
weights, which provides a unique evidence of the existence of seasonal variability in
ED
the PS loadings.
Table 6 shows the correlations between the model weights and parameters in
PT
BAIPU. Although statistically significant in general, the correlations in the
CE
real-situation case are not as strong as those in the synthetic case, indicating a weaker
AC
interaction between the input and parametric uncertainties. In addition, Table S3 in the Auxiliary Materials shows the correlations between different model weights. In this real-situation case, BAIPU and the sequential approaches led to slightly different results. The predictive uncertainty bands of the nitrate response are shown in Figure 11. The median predictions (black lines) in Figures 11a and 11b adequately match the 38
ACCEPTED MANUSCRIPT
observations. In the BAIPU and the sequential approach, the 50% confidence intervals encompasses 50.24% and 51.18% of the observations, respectively; the 95% confidence intervals encompasses 96.21% and 96.68% of the observations, respectively. Figures 11c and 11d further compare the predictive means and standard
CR IP T
deviations derived by the two approaches. In general, the BAIPU produces the same predictive means as the sequential approach but smaller standard deviations. This
indicates that the BAIPU reduces the predictive uncertainty in this case, compared to
AN US
the sequential approach.
The uncertainty results presented in Figure 11 and Figure 7b suggest that the two approaches, as well as those individual MCMC analyses, have a similar performance
M
during the Bayesian calibration. However, as revealed in the synthetic case, the
ED
difference could be more substantial in predicting the future water response. This real-situation modeling case also considers three future loading scenarios, as
PT
introduced in Table 3. Figure 12 shows that the BAIPU and the sequential BMA
CE
approach only have a slight difference in the frequency distributions, which indicates
AC
that the interaction between parameters and inputs in the Bayesian analysis is not as significant as that in the synthetic case. The difference between individual models are more significant, but not as significant as that in the synthetic case (Figure 8). This is because two time-variant loading scenarios are considered in this real-situation case, while the synthetic case considers five scenarios of constant loadings. Thus, in the WWQ modeling, adequately accounting for the temporal variability of source 39
ACCEPTED MANUSCRIPT
loadings is critical to the results of stochastic modeling. All the above results clearly demonstrate the applicability and usefulness of BAIPU in the context of real-world water quality management.
5. Conclusions
CR IP T
This study develops a new approach, BAIPU, for the joint analysis of input and parametric uncertainties through a tight coupling of MCMC and BMA. The formal
likelihood function for this approach is also derived considering a lag-1 autocorrelated,
AN US
heteroscedastic, and SEP distributed error model. Input uncertainty in a discrete form is specifically addressed, which is commonly encountered in environmental modeling for real-world problems. A series of numerical experiments were designed and
M
implemented based on a synthetic nitrate pollution case, as well as the real pollution
ED
case in the Newport Bay Watershed, California. The major conclusions include the following. First, with parallel computing, the BAIPU can be implemented with
PT
reasonable computational costs and appropriately identify the uncertain parameters
CE
and characterize the predictive uncertainty. Second, if significant input uncertainty is
AC
not accounted for, WWQ modeling may seriously mislead management decisions because the input uncertainty must be compensated for by the parametric uncertainty. Third, the BAIPU accounts for the interaction between the input and parametric uncertainties through joint analysis and is mathematically more rigorous than the sequential approach. Thus, in theory, it can provide more accurate calibration and uncertainty results than the sequential approach. Finally, the BAIPU quantifies the 40
ACCEPTED MANUSCRIPT
credibility (i.e., the model weights) of different input assumptions on a statistical basis, and it can be implemented as an inverse modeling approach for joint inference of both model parameters and inputs. The BAIPU was originally designed to address discrete form input uncertainty,
CR IP T
but future studies may expand the BAIPU to a joint analysis of model structural, input and parametric uncertainties. In addition, the BAIPU has a great potential of being applied in fields other than watershed water quality modeling. For example, in
AN US
basin-scale groundwater modeling, boundary conditions of subsurface inflow are
often highly uncertain [72], which could be addressed in the framework of BAIPU.
M
Acknowledgments
ED
This work was supported by the National Natural Science Foundation of China (No. 41622111; No. 91647201) and the China Postdoctoral Science Foundation
PT
(2017M612505). Additional support was provided by the Southern University of
CE
Science and Technology (No. G01296001). We thank Dr. Jasper A. Vrugt at the
AC
University of California, Irvine for kindly discussing with us on DREAM(ZS).
41
ACCEPTED MANUSCRIPT
References [1] Borah DK, M Bera. Watershed-scale hydrologic and nonpoint-source pollution models: Review of mathematical bases. American Society of Agricultural Engineers. 46 (2003) 1553-66, doi: 10.13031/2013.10398. [2] Arnold JG, R Srinivasan, RS Muttiah, JR Williams. Large area hydrologic modeling
Resources Association. 34 (1998) 73-89, doi: 10.1111/j.1752-1688.1998.tb05961.x.
CR IP T
and assessment part I: Model development. Journal of the American Water
[3] Neitsch SL, JG Arnold, JR Kiniry, JR Williams. Soil and water assessment tool
AN US
theoretical documentation version 2009. Texas Water Resources Institute & Texas A&M University 2011. http://swat.tamu.edu/media/99192/swat2009-theory.pdf. Accessed Oct. 14, 2017.
[4] Chen C, J Herr, L Weintraub. Decision support system for stakeholder involvement.
M
Journal of Environmental Engineering. 130 (2004) 714-21, doi: doi:10.1061/(ASCE)0733-9372(2004)130:6(714).
ED
[5] Zheng Y, W Wang, F Han, J Ping. Uncertainty assessment for watershed water quality modeling: A probabilistic collocation method based approach. Advances in
PT
Water Resources. 34 (2011) 887-98, doi: 10.1016/j.advwatres.2011.04.016. [6] Bicknell BR, JC Imhoff, JJL Kittle, JAS Donigian, RC Johanson. Hydrological
CE
simulation program -- Fortran, user's manual for version 11. U.S. Environmental Protection Agency 1997.
AC
[7] Young RA, CA Onstad, DD Bosch, WP Anderson. AGNPS: A nonpoint-source pollution model for evaluating agricultural watersheds. Journal of Soil and Water Conservation. 44 (1989) 168-73.
[8] Bingner RL, FD Theurer, Y Yuan. AnnAGNPS technical processes: version 5.4. 2015. https://www.wcc.nrcs.usda.gov/ftpref/wntsc/H&H/AGNPS/downloads/ AnnAGNPS_Technical_Documentation.pdf. Accessed Oct. 14, 2017. [9] Jeon J-H, CG Yoon, AS Donigian, K-W Jung. Development of the HSPF-Paddy 42
ACCEPTED MANUSCRIPT
model to estimate watershed pollutant loads in paddy farming regions. Agricultural Water Management. 90 (2007) 75-86, doi: 10.1016/j.agwat.2007.02.006. [10] Huang J, H Hong. Comparative study of two models to simulate diffuse nitrogen and phosphorus pollution in a medium-sized watershed, southeast China. Estuarine, Coastal and Shelf Science. 86 (2010) 387-94, doi: 10.1016/j.ecss.2009.04.003.
CR IP T
[11] Pease LM, P Oduor, G Padmanabhan. Estimating sediment, nitrogen, and phosphorous loads from the Pipestem Creek watershed, North Dakota, using AnnAGNPS. Computers & Geosciences. 36 (2010) 282-91, doi: 10.1016/j.cageo.2009.07.004.
[12] Volk M, D Bosch, V Nangia, B Narasimhan. SWAT: Agricultural water and
AN US
nonpoint source pollution management at a watershed scale. Agricultural Water Management. 175 (2016) 1-3, doi: 10.1016/j.agwat.2016.06.013. [13] Volk M, D Bosch, V Nangia, B Narasimhan. SWAT: Agricultural water and
nonpoint source pollution management at a watershed scale—Part II. Agricultural
M
Water Management. 180 (2017) 191-3, doi: 10.1016/j.agwat.2016.09.029.
ED
[14] Cho J, C Oh, J Choi, Y Cho. Climate change impacts on agricultural non-point source pollution with consideration of uncertainty in CMIP5. Irrigation and
PT
Drainage. 65 (2016) 209-20, doi: 10.1002/ird.2036. [15] Nguyen HH, F Recknagel, W Meyer, J Frizenschaf, MK Shrestha. Modelling the
CE
impacts of altered management practices, land use and climate changes on the water quality of the Millbrook catchment-reservoir system in South Australia.
AC
Journal of Environmental Management. 202 (2017) 1-11, doi: 10.1016/j.jenvman.2017.07.014.
[16] Karabulut A, BN Egoh, D Lanzanova, B Grizzetti, G Bidoglio, L Pagliero, et al. Mapping water provisioning services to support the ecosystem–water–food– energy nexus in the Danube river basin. Ecosystem Services. 17 (2016) 278-92, doi: 10.1016/j.ecoser.2015.08.002. [17] Psomas A, V Dagalaki, Y Panagopoulos, D Konsta, M Mimikou. Sustainable 43
ACCEPTED MANUSCRIPT
agricultural water management in Pinios River basin using remote sensing and hydrologic modeling. Procedia Engineering. 162 (2016) 277-83, doi: 10.1016/j.proeng.2016.11.059. [18] Her Y, J Jeong, J Arnold, L Gosselink, R Glick, F Jaber. A new framework for modeling decentralized low impact developments using Soil and Water
10.1016/j.envsoft.2017.06.005.
CR IP T
Assessment Tool. Environmental Modelling & Software. 96 (2017) 305-22, doi:
[19] Zheng Y, AA Keller. Uncertainty assessment in watershed-scale water quality modeling and management: 1. Framework and application of generalized
likelihood uncertainty estimation (GLUE) approach. Water Resources Research.
AN US
43 (2007), doi: 10.1029/2006wr005345.
[20] Rode M, G Arhonditsis, D Balin, T Kebede, V Krysanova, A van Griensven, et al. New challenges in integrated water quality modelling. Hydrological Processes. 24 (2010) 3447-61, doi: 10.1002/hyp.7766.
M
[21] Beck MB. Water quality modeling: A review of the analysis of uncertainty. Water
ED
Resources Research. 23 (1987) 1393-442, doi: 10.1029/WR023i008p01393. [22] van Griensven A, T Meixner. Methods to quantify and identify the sources of
PT
uncertainty for river basin water quality models. Water Science & Technology. 53 (2006) 51-9, doi: 10.2166/wst.2006.007.
CE
[23] Knoche M, C Fischer, E Pohl, P Krause, R Merz. Combined uncertainty of hydrological model complexity and satellite-based forcing data evaluated in two
AC
data-scarce semi-arid catchments in Ethiopia. Journal of Hydrology. 519 (2014) 2049-66, doi: 10.1016/j.jhydrol.2014.10.003.
[24] Han F, Y Zheng. Multiple-response Bayesian calibration of watershed water quality models with significant input and model structure errors. Advances in Water Resources. 88 (2015) 109-23, doi: 10.1016/j.advwatres.2015.12.007. [25] Gelman A, JB Carlin, HS Stern, DB Rubin. Bayesian data analysis. Chapman & Hall/CRC, United States of America, 2003. 44
ACCEPTED MANUSCRIPT
[26] Robert CP, G Casella. Monte Carlo statistical methods. Springer, United States of America, 2004. [27] Kuczera G, E Parent. Monte Carlo assessment of parameter uncertainty in conceptual catchment models: the Metropolis algorithm. Journal of Hydrology. 211 (1998) 69-85, doi: 10.1016/s0022-1694(98)00198-x.
CR IP T
[28] Bates BC, EP Campbell. A Markov Chain Monte Carlo scheme for parameter estimation and inference in conceptual rainfall-runoff modeling. Water Resources Research. 37 (2001) 937-47, doi: 10.1029/2000wr900363.
[29] Marshall L, D Nott, A Sharma. A comparative study of Markov chain Monte Carlo methods for conceptual rainfall-runoff modeling. Water Resources Research. 40
AN US
(2004), doi: 10.1029/2003wr002378.
[30] Yang J, P Reichert, KC Abbaspour, H Yang. Hydrological modelling of the Chaohe basin in China: Statistical model formulation and Bayesian inference. Journal of Hydrology. 340 (2007) 167-82, doi: 10.1016/j.jhydrol.2007.04.006.
M
[31] Smith TJ, LA Marshall. Bayesian methods in hydrologic modeling: A study of
ED
recent advancements in Markov chain Monte Carlo techniques. Water Resources Research. 44 (2008), doi: 10.1029/2007wr006705.
PT
[32] Vrugt JA, CJF ter Braak, HV Gupta, BA Robinson. Equifinality of formal (DREAM) and informal (GLUE) Bayesian approaches in hydrologic modeling?
CE
Stochastic Environmental Research and Risk Assessment. 23 (2008) 1011-26, doi: 10.1007/s00477-008-0274-y.
AC
[33] Box GEP, GC Tiao. Bayesian inference in statistical analysis. Wiley, United States of America & Canada, 1992.
[34] Kuczera G. Improved parameter inference in catchment models: 1. Evaluating parameter uncertainty. Water Resources Research. 19 (1983) 1151-62, doi: 10.1029/WR019i005p01151. [35] Schaefli B, DB Talamba, A Musy. Quantifying hydrological modeling errors through a mixture of normal distributions. Journal of Hydrology. 332 (2007) 45
ACCEPTED MANUSCRIPT
303-15, doi: 10.1016/j.jhydrol.2006.07.005. [36] Schoups G, JA Vrugt. A formal likelihood function for parameter and predictive inference of hydrologic models with correlated, heteroscedastic, and non-Gaussian errors. Water Resources Research. 46 (2010), doi: 10.1029/2009wr008933. [37] Smith T, A Sharma, L Marshall, R Mehrotra, S Sisson. Development of a formal
CR IP T
likelihood function for improved Bayesian inference of ephemeral catchments. Water Resources Research. 46 (2010), doi: 10.1029/2010wr009514.
[38] Evin G, D Kavetski, M Thyer, G Kuczera. Pitfalls and improvements in the joint inference of heteroscedasticity and autocorrelation in hydrological model calibration. Water Resources Research. 49 (2013) 4518-24, doi:
AN US
10.1002/wrcr.20284.
[39] Honti M, C Stamm, P Reichert. Integrated uncertainty assessment of discharge predictions with a statistical error model. Water Resources Research. 49 (2013) 4866-84, doi: 10.1002/wrcr.20374.
M
[40] Renard B, D Kavetski, G Kuczera, M Thyer, SW Franks. Understanding predictive
ED
uncertainty in hydrologic modeling: The challenge of identifying input and structural errors. Water Resources Research. 46 (2010), doi:
PT
10.1029/2009wr008328.
[41] Kavetski D, G Kuczera, SW Franks. Bayesian analysis of input uncertainty in
CE
hydrological modeling: 1. Theory. Water Resources Research. 42 (2006), doi: 10.1029/2005wr004368.
AC
[42] Ajami NK, Q Duan, S Sorooshian. An integrated hydrologic Bayesian multimodel combination framework: Confronting input, parameter, and model structural uncertainty in hydrologic prediction. Water Resources Research. 43 (2007), doi: 10.1029/2005wr004745. [43] Yen H, X Wang, DG Fontane, RD Harmel, M Arabi. A framework for propagation of uncertainty contributed by parameterization, input data, model structure, and calibration/validation data in watershed modeling. Environmental Modelling & 46
ACCEPTED MANUSCRIPT
Software. 54 (2014) 211-21, doi: 10.1016/j.envsoft.2014.01.004. [44] Vrugt JA, CJF ter Braak, MP Clark, JM Hyman, BA Robinson. Treatment of input uncertainty in hydrologic modeling: Doing hydrology backward with Markov chain Monte Carlo simulation. Water Resources Research. 44 (2008), doi: 10.1029/2007wr006720.
CR IP T
[45] Li M, D Yang, J Chen, SS Hubbard. Calibration of a distributed flood forecasting model with input uncertainty using a Bayesian framework. Water Resources Research. 48 (2012), doi: 10.1029/2010wr010062.
[46] Thyer M, B Renard, D Kavetski, G Kuczera, SW Franks, S Srikanthan. Critical evaluation of parameter consistency and predictive uncertainty in hydrological
AN US
modeling: A case study using Bayesian total error analysis. Water Resources Research. 45 (2009), doi: 10.1029/2008wr006825.
[47] Raftery AE, T Gneiting, F Balabdaoui, M Polakowski. Using Bayesian Model Averaging to calibrate forecast ensembles. Monthly Weather Review. 133 (2005)
M
1155-74, doi: 10.1175/mwr2906.1.
ED
[48] Duan Q, NK Ajami, X Gao, S Sorooshian. Multi-model ensemble hydrologic prediction using Bayesian model averaging. Advances in Water Resources. 30
PT
(2007) 1371-86, doi: 10.1016/j.advwatres.2006.11.014. [49] Rojas R, L Feyen, A Dassargues. Conceptual model uncertainty in groundwater
CE
modeling: Combining generalized likelihood uncertainty estimation and Bayesian model averaging. Water Resources Research. 44 (2008), doi:
AC
10.1029/2008wr006908.
[50] Zhang X, R Srinivasan, D Bosch. Calibration and uncertainty analysis of the SWAT model using Genetic Algorithms and Bayesian Model Averaging. Journal of Hydrology. 374 (2009) 307-17, doi: 10.1016/j.jhydrol.2009.06.023. [51] Bastola S, C Murphy, J Sweeney. The role of hydrological modelling uncertainties in climate change impact assessments of Irish river catchments. Advances in Water Resources. 34 (2011) 562-76, doi: 10.1016/j.advwatres.2011.01.008. 47
ACCEPTED MANUSCRIPT
[52] Zeng X, J Wu, D Wang, X Zhu, Y Long. Assessing Bayesian model averaging uncertainty of groundwater modeling based on information entropy method. Journal of Hydrology. 538 (2016) 689-704, doi: 10.1016/j.jhydrol.2016.04.038. [53] Wellen C, GB Arhonditsis, T Long, D Boyd. Quantifying the uncertainty of nonpoint source attribution in distributed water quality models: A Bayesian
(2014) 3353-68, doi: 10.1016/j.jhydrol.2014.10.007.
CR IP T
assessment of SWAT’s sediment export predictions. Journal of Hydrology. 519
[54] Leamer EE. Specification searches: Ad hoc inference with nonexperimental data. Wiley, United States of America, 1978.
[55] Hoeting JA, D Madigan, AE Raftery, CT Volinsky. Bayesian model averaging: A
AN US
tutorial. Statistical Science. 14 (1999) 382-417, doi: 10.1214/ss/1009212519.
[56] Schoniger A, T Wohling, W Nowak. A statistical concept to assess the uncertainty in Bayesian model weights and its impact on model ranking. Water Resources Research. 51 (2015) 7524-46, doi: 10.1002/2015WR016918.
M
[57] Tsai FTC, X Li. Inverse groundwater modeling for hydraulic conductivity
ED
estimation using Bayesian model averaging and variance window. Water Resources Research. 44 (2008) W09434, doi: 10.1029/2007WR006576.
PT
[58] Volpi E, G Schoups, G Firmani, JA Vrugt. Sworn testimony of the model evidence: Gaussian Mixture Importance (GAME) sampling. Water Resources Research. 53
CE
(2017) 6133-58, doi: 10.1002/2016WR020167. [59] Strauch M, C Bernhofer, S Koide, M Volk, C Lorz, F Makeschin. Using
AC
precipitation data ensemble for uncertainty analysis in SWAT streamflow simulation. Journal of Hydrology. 414-415 (2012) 413-24, doi: 10.1016/j.jhydrol.2011.11.014.
[60] Vrugt JA, CJF ter Braak, CGH Diks, BA Robinson, JM Hyman, D Higdon. Accelerating Markov Chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling. Int J Nonlin Sci Num. 10 (2009) 273-90, doi: 10.1515/ijnsns.2009.10.3.273. 48
ACCEPTED MANUSCRIPT
[61] Laloy E, JA Vrugt. High-dimensional posterior exploration of hydrologic models using multiple-try DREAM(ZS) and high-performance computing. Water Resources Research. 48 (2012), doi: 10.1029/2011wr010608. [62] Zheng Y, F Han. Markov Chain Monte Carlo (MCMC) uncertainty analysis for watershed water quality modeling and management. Stochastic Environmental
CR IP T
Research and Risk Assessment. 30 (2015) 293-308, doi: 10.1007/s00477-015-1091-8.
[63] ter Braak CJF, JA Vrugt. Differential Evolution Markov Chain with snooker updater and fewer chains. Statistics and Computing. 18 (2008) 435-46, doi: 10.1007/s11222-008-9104-9.
AN US
[64] Gelman A, DB Rubin. Inference from iterative simulation using multiple
sequences. Statistical Science. 7 (1992) 457-72, doi: 10.1214/ss/1177011136. [65] Dempster AP, NM Laird, DB Rubin. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B
M
(Methodological)1977. pp. 1-38.
ED
[66] Vrugt JA, CGH Diks, MP Clark. Ensemble Bayesian model averaging using Markov Chain Monte Carlo sampling. Environmental Fluid Mechanics. 8 (2008)
PT
579-595, doi: 10.1007/s10652-008-9106-3. [67] U.S. Environmental Protection Agency. Total Maximum Daily Loads for Nutrients
CE
San Diego Creek and Newport Bay, California. 1998. http://www. epa.gov/owow/tmdl/examples/nutrients/ca_sdnbay.pdf. Accessed Sep. 10, 2014.
AC
[68] U.S. Environmental Protection Agency. 1998 California 303(d) list and TMDL priority schedule. 1999. http://www.krisweb.com/biblio/ ncc_swrcb_ncrwqcb_2001_wqcontrolplanapdx.pdf. Accessed Oct. 14, 2017.
[69] Zheng Y, AA Keller. Stochastic watershed water quality simulation for TMDL development - a case study in the Newport Bay watershed. Journal of the American Water Resources Association. 44 (2008) 1397-410, doi: 10.1111/j.1752-1688.2008.00232.x. 49
ACCEPTED MANUSCRIPT
[70] Balin Talamba D, E Parent, A Musy. Bayesian multiresponse calibration of TOPMODEL: Application to the Haute-Mentue catchment, Switzerland. Water Resources Research. 46 (2010), doi: 10.1029/2007wr006449. [71] Campolongo F, J Cariboni, A Saltelli. An effective screening design for sensitivity analysis of large models. Environmental Modelling & Software. 22 (2007)
CR IP T
1509-18, doi: 10.1016/j.envsoft.2006.10.004. [72] Wu B, Y Zheng, Y Tian, X Wu, Y Yao, F Han, et al. Systematic assessment of the uncertainty in integrated surface water-groundwater modeling based on the
probabilistic collocation method. Water Resources Research. 50 (2014) 5848-65,
AC
CE
PT
ED
M
AN US
doi: 10.1002/2014wr015366.
50
ACCEPTED MANUSCRIPT
Tables Table 1. Different SWAT models in the synthetic and real-situation modeling cases True model
Model 1
Model 2
Model 3
Model 4
Model 5
/
RM1
RM2
RM3
RM4
RM5
Point source A a
/
48.8 (C)
61 (C)
73.2 (C)
61 (SF)
61 (BF)
a
/
24.8 (C)
31 (C)
37.2 (C)
31 (SF)
31 (BF)
Fertilizer use b
/
216 (HU)
240 (HU)
264 (HU)
240 (HU)
240 (HU)
Stormwater model
/
USGS
USGS
USGS
USGS
USGS
Model name
Point source B
Synthetic case TM
IM1
IM2
IM3
IM4
IM5
Point source A
a
132 (D)
100 (C)
125 (C)
150 (C)
175 (C)
200 (C)
Point source B
a
68 (D)
50 (C)
62.5 (C)
75 (C)
87.5 (C)
100 (C)
240 (Date)
165 (HU)
206.3 (HU)
247.5 (HU)
288.8 (HU)
330 (HU)
buildup-washoff
USGS
USGS
USGS
USGS
USGS
Fertilizer use
b
Stormwater model a
AN US
Model ID
CR IP T
Real-situation case
The average PS loading of nitrate pollution (kg N/day). "D" means the true PS loading in the synthetic case is Dynamic or time-variant. "C" represents a Constant loading input. "SF" and Fluctuations, respectively.
b
M
"BF" in the real-situation case indicate artificial time-variant loading inputs with Small and Big Fertilization rate in the NBW (kg N/ha/year). "Date" means the fertilizer applications are
AC
CE
PT
ED
scheduled on specific dates, and "HU" indicates the heat unit scheduling approach.
51
ACCEPTED MANUSCRIPT
Table 2. Uncertain SWAT parameters considered in the Bayesian inference Definition
Type
Range
ESCO
Soil evaporation compensation factor (/)
I
[0.01, 1]
Synthetic True value 0.02312
CN2
SCS runoff curve number for moisture condition II (/)
II
[-20, 15]
-9.770
SURLAG
Surface runoff lag coefficient (/)
I
[0.1, 12]
2.235
SOL_AWC
Available water capacity of the soil layer (mm H2O/mm soil)
II
[-0.2, 0.2]
0.1103
CH_K1
Hydraulic conductivity in tributary channel alluvium (mm/h)
I
[1, 25]
1.458
CH_N1
Manning's "n" value for the tributary channels (/)
I
[0.01, 0.2]
0.0511
CH_K2
Hydraulic conductivity in main channel alluvium (mm/h)
I
[1, 25]
7.094
ALPHA_BNK
Baseflow alpha factor for bank storage (day)
I
[0.001, 1]
0.0047
SHALLST
Initial depth of water in the shallow aquifer (mm H2O)
I
[500, 2000]
644.60
GW_DELAY
Groundwater delay time (day)
I
[10, 800]
723.69
NPERCO
Nitrogen percolation coefficient (/)
I
[0.01, 1]
0.0843
SDNCO
Denitrification threshold water content (/)
I
[0.2, 1.2]
0.2556
AC
CE
PT
ED
M
AN US
CR IP T
Parameters
52
ACCEPTED MANUSCRIPT
Table 3. Hypothetical future loading scenarios for predicting potential water quality responses Loading scenario
Synthetic case
Real-situation case
Point source A a
210 kg N/day
90 kg N/day
a
105 kg N/day
45 kg N/day
372 kg N/ha/year
356 kg N/ha/year
Point source B Fertilizer use
b
Scenario 2 (steady loading) Point source A a Point source B Fertilizer use
a
b
140 kg N/day
60 kg N/day
70 kg N/day
30 kg N/day
248 kg N/ha/year
238 kg N/ha/year
Scenario 3 (decreased loading) 70 kg N/day
30 kg N/day
a
35 kg N/day
15 kg N/day
Point source B Fertilizer use
b
124 kg N/ha/year
AN US
Point source A a
CR IP T
Scenario 1 (increased loading)
119 kg N/ha/year
Constant loading.
b
Fertilizer application is scheduled on specific dates.
AC
CE
PT
ED
M
a
53
ACCEPTED MANUSCRIPT
Table 4. Efficiency of DREAM(ZS) in both the synthetic and real-situation cases AR a
NC b
[θ, φ, w]
5.17%
489,840
w
62.83%
2640
MCMC_IM1
[θ, φ]
2.99%
350,160
MCMC_IM2
[θ, φ]
2.96%
328,320
MCMC_IM3
[θ, φ]
4.02%
288,000
MCMC_IM4
[θ, φ]
3.97%
357,600
MCMC_IM5
[θ, φ]
1.88%
325,680
BAIPU_RMs
[θ, φ, w]
23.19%
143,040
w
66.25%
2,400
MCMC_RM1
[θ, φ]
4.69%
141,840
MCMC_RM2
[θ, φ]
6.27%
143,760
MCMC_RM3
[θ, φ]
5.47%
107,760
MCMC_RM4
[θ, φ]
5.39%
173,040
MCMC_RM5
[θ, φ]
6.84%
100,080
BMA_IMs Synthetic case
BMA_RMs Real-situation case
AN US
BAIPU_IMs
CR IP T
Parameters
Approach
AR is the acceptance rate.
b
NC is the number of samples for the Markov chains to achieve convergence.
AC
CE
PT
ED
M
a
54
ACCEPTED MANUSCRIPT Table 5. Correlation coefficients a between model parameters and the weights derived by BAIPU in the synthetic case w2
w3
w4
w5
ESCO
0.08
0.06
0.07
-0.12
0.02
CN2
0.08
0.06
0.00
0.01
-0.08
SURLAG
0.02
0.04
-0.08
0.09
-0.07
SOL_AWC
0.06
0.03
0.11
0.03
-0.17
CH_K1
-0.05
-0.13
0.00
0.00
0.10
CH_N1
0.02
0.05
-0.08
0.09
-0.08
CH_K2
-0.02
0.04
-0.04
0.05
ALPHA_BNK
0.03
-0.03
0.00
0.04
SHALLST
-0.27
-0.12
-0.14
0.02
GW_DELAY
-0.12
-0.07
-0.07
0.11
NPERCO
-0.08
-0.07
-0.09
-0.04
SDNCO
0.17
0.07
0.12
0.07
-0.03 -0.05 0.27 0.01 0.18
-0.29
AN US
a
CR IP T
w1
5
Because the sample size is huge (above 10 ), all the correlations are statistically significant at the
confidence level of 99% (i.e., p<0.01), except those with 0.00 coefficients. The coefficients larger
AC
CE
PT
ED
M
than 0.1 or smaller than -0.1 are shaded in the table to indicate relatively strong correlations.
55
ACCEPTED MANUSCRIPT Table 6. Correlation coefficients a between model parameters and the weights derived by the BAIPU in the real-situation case w2
w3
w4
w5
ESCO
0.04
-0.01
-0.05
-0.03
0.03
CN2
0.03
-0.02
0.01
0.00
-0.02
SURLAG
-0.05
0.00
0.01
0.00
0.04
SOL_AWC
0.09
0.00
-0.07
0.04
-0.09
CH_K1
-0.04
0.05
-0.02
-0.01
0.02
CH_N1
-0.03
0.02
0.00
0.00
0.02
CH_K2
-0.02
0.02
0.02
0.01
ALPHA_BNK
0.00
-0.01
-0.02
0.03
SHALLST
-0.09
0.02
-0.01
0.06
GW_DELAY
-0.11
-0.01
0.01
0.09
NPERCO
-0.23
-0.04
0.07
0.03
SDNCO
0.13
0.00
-0.10
-0.03
-0.02 -0.02 -0.01 -0.02 0.14
-0.03
AN US
a
CR IP T
w1
5
Because the sample size is huge (above 10 ), all the correlations are statistically significant at the
confidence level of 99% (i.e., p<0.01), except those with 0.00 coefficients. The coefficients larger
AC
CE
PT
ED
M
than 0.1 or smaller than -0.1 are shaded in the table to indicate relatively strong correlations.
56
ACCEPTED MANUSCRIPT
Figure captions Figure 1. The ―returning‖ operation proposed for sampling model weights using DREAM(ZS). Figure 2. The framework of the Bayesian Analysis of Input and Parametric
CR IP T
Uncertainties (BAIPU).
Figure 3. The study area and modeling domain. (a) The Newport Bay Watershed (NBW) and (b) location of the NBW.
AN US
Figure 4. Cumulative distribution functions (CDFs) of model parameters inferred by the BAIPU and MCMC analyses for the individual models in the synthetic case. (a)-(l) are for the twelve different SWAT parameters.
M
Figure 5. Marginal probability distribution functions of model weights derived based
ED
on BAIPU (red bars) and the sequential BMA approach (grey bars) in the synthetic case. Red points and black circles indicate the respective mean values. (a)-(e)
PT
correspond to the five models.
CE
Figure 6. Results of the uncertainty analysis in the synthetic case. (a) and (b) illustrate
AC
the uncertainty bands (50% and 95% confidence intervals) of nitrate concentration based on BAIPU and the sequential BMA approach, respectively. (c) and (d) compare the means and standard deviations of the predictive distributions based on the two approaches. Each circle in the subplots (c) and (d) represents a simulation day. Figure 7. Quantile-Quantile (QQ) plots of nitrate concentration based on the BAIPU, sequential BMA and individual MCMC analyses. (a) and (b) are for the synthetic and 57
ACCEPTED MANUSCRIPT
real-situation cases, respectively. The solid black lines are the 45-degree lines. Figure 8. True and simulated values of the two management variables under the three future loading scenarios in the synthetic modeling case. (a), (c) and (e) are the results
exceeding (FOE) 6 mg N/L.
CR IP T
of average nitrate concentration (AVC). (b), (d) and (f) are the results of frequency of
Figure 9. Cumulative distribution functions (CDFs) of uncertain model parameters
inferred by the BAIPU and individual MCMC analyses in the real-situation modeling
AN US
case. (a)-(i) are for the twelve model parameters.
Figure 10. Marginal probability distribution functions of model weights derived based on BAIPU (red bars) and the sequential BMA approach (grey bars) in the
M
real-situation modeling case. Red points and black circles indicate the respective
ED
mean values. (a)-(e) correspond to the five models. Figure 11. Results of the uncertainty analysis in the real-situation modeling case. (a)
PT
and (b) illustrate the uncertainty bands (50% and 95% confidence intervals) of nitrate
CE
concentration based on the BAIPU and the sequential BMA approach, respectively. (c)
AC
and (d) compare the means and standard deviations of the predictive distributions based on the two approaches. Each circle in the subplots (c) and (d) represents a simulation day. Figure 12. Simulated values of the two management-relevant variables under the three future loading scenarios in the real-situation modeling case. (a), (c) and (e) are the results of average nitrate concentration (AVC). (b), (d) and (f) are the results of 58
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
frequency of exceeding (FOE) 6 mg N/L.
59
ACCEPTED MANUSCRIPT
-∆o
∆o
θ*,i
θb
θ**,i ∆
boundary
boundary
CR IP T
θi
Figure 1. The ―returning‖ operation proposed for sampling model weights using
AC
CE
PT
ED
M
AN US
DREAM(ZS).
60
ACCEPTED MANUSCRIPT
Input X2
. . .
Input XK
Model M1
Model M2
. . .
Model MK
Hyperparameters φ
MLE: EM algorithm
The residual error model
p(Y|X1,θ,φ)
CR IP T
Model parameters θ
Input X1
p(Y|X2,θ,φ)
w1(θ,φ)
w2(θ,φ)
MCMC sampling
. . .
p(Y|XK,θ,φ)
. . . wK(θ,φ)
AN US
BMA prediction: p(Y|X1,...,K,θ,φ)
No
ED
MCMC sampling
Likelihood function: l(θ,φ,w1,...K,X1,...K,Z)
M
Observations Z
Has the Markov chain converged ? Yes
Stop the sampling
PT
Predictive distribution
CE
Figure 2. The framework of the Bayesian Analysis of Input and Parametric
AC
Uncertainties (BAIPU).
61
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 3. The study area and modeling domain. (a) The Newport Bay Watershed
AC
CE
PT
ED
M
(NBW) and (b) location of the NBW.
62
ACCEPTED MANUSCRIPT
1
1
1
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
(a) ESCO
0 0
0.1
0.2
0.3
0.4
0.5
CDF
0.8
(d) SOL_AWC
0.2
0.3
0.4
0.8
0.8
0.4
0.4
0.2
0.2 0.5
0.7
0.9 1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
(e) CH_K1
0.3
0.4
0.5
1 (j) GW_DELAY
0.6
0.1
0
0.4
PT
0
0.6
0.8
0.2
0.8
11
0
1
1
0.8
0.8
CE
0
0.4
0.6
0.8
1
0.6
0.8
1
(l) SDNCO 0.2 0
0
1 0
0.2
0.4 (k) NPERCO
0.05
0.1
0.15
0.2
0
MCMC_IM1
MCMC_IM2
MCMC_IM5
True value
0.04
AC
MCMC_IM4
(i) SHALLST
0.6
0
BAIPU_IMs
0.4
0
0.0070.6
0.2
1
0.3
0.4
0.4
0.2
0.2
0.6
0.6
0.4
0.1
0.8
(h) ALPHA_BNK
ED
0.8
0.025 0.05 0.075
M
0.2
(f) CH_N1
1
0 0.003
0
0.25
0
0.2
(g) CH_K2
0.2
0.4
0
1
0.15
0.6
0 0.3
0.1 1
0.6
0.1
CDF
0.1
1
0.6
0
CDF
0 0
1
(c) SURLAG
CR IP T
0
0.2
(b) CN2
AN US
CDF
0.8
0.2
0.4
MCMC_IM3
0.08
Figure 4. Cumulative distribution functions (CDFs) of model parameters inferred by the BAIPU and MCMC analyses for the individual models in the synthetic case. (a)-(l) are for the twelve different SWAT parameters.
63
ACCEPTED MANUSCRIPT
0.2
0.1
0.0 0.0
0.5
0.5
1.0
(c) w3 0.2
0.0 0.0
0.5
1.0
0.6
(d) w4 0.3
0.2
0.4
0.6
(e) w5
Sequential BMA BAIPU
0.3
CR IP T
Posterior density
Posterior density
0.1
0.0 0.0
1.0
0.6
0.0 0.0
0.4
(b) w2
Posterior density
(a) w1
Posterior density
Posterior density
0.2
Synthetic case
0.0 0.0
0.2
0.4
AN US
Figure 5. Marginal probability distribution functions of model weights derived based on BAIPU (red bars) and the sequential BMA approach (grey bars) in the synthetic case. Red points and black circles indicate the respective mean values. (a)-(e)
AC
CE
PT
ED
M
correspond to the five models.
64
ACCEPTED MANUSCRIPT
60
95% conf. int.
-2 -12 Median
NO3-N (mg N/L)
50
50% conf. int.
(a) BAIPU
Observations
40 30 20 10
CR IP T
0
60
(b) Sequential BMA 40
30
AN US
20 10 0
16
(c) Predictive Mean
M
20
ED
15 10 5 0 0
PT
By sequential BMA (mg N/L)
25
5
10
15
By sequential BMA (mg N/L)
NO3-N (mg N/L)
50
(d) Predictive std. dev.
12
8
4
0 20
0
25
4
8
12
16
By BAIPU (mg N/L)
CE
By BAIPU (mg N/L)
AC
Figure 6. Results of the uncertainty analysis in the synthetic case. (a) and (b) illustrate the uncertainty bands (50% and 95% confidence intervals) of nitrate concentration based on BAIPU and the sequential BMA approach, respectively. (c) and (d) compare the means and standard deviations of the predictive distributions based on the two approaches. Each circle in the subplots (c) and (d) represents a simulation day.
65
ACCEPTED MANUSCRIPT
1
0.8
0.6 0.4 0.2
0
(b) NO3; Real-situation case
0.8
0.6 0.4 0.2
0 0
0.2 0.4 0.6 0.8 Theoretical Quantile of U[0,1]
1
0
0.2 0.4 0.6 0.8 Theoretical Quantile of U[0,1]
1 BAIPU_IMs(RMs)
BMA_IMs(RMs)
MCMC_IM2(RM2) 0
MCMC_IM3(RM3)
0
CR IP T
(a) NO3; Synthetic case
Quantile of observed p-values
Quantile of observed p-values
1
1
MCMC_IM1(RM1) MCMC_IM4(RM4) 0.1
AN US
MCMC_IM5(RM5)
Figure 7. Quantile-Quantile (QQ) plots of nitrate concentration based on the BAIPU, sequential BMA and individual MCMC analyses. (a) and (b) are for the synthetic and
AC
CE
PT
ED
M
real-situation cases, respectively. The solid black lines are the 45-degree lines.
66
ACCEPTED MANUSCRIPT
1
(a) Scenario 1; AVC
0.8
Cumulative frequency
0.6 0.4 0.2 0 12
14
16
18
0.6 0.4 0.2
1
0.8
Cumulative frequency
(c) Scenario 2; AVC
0.6 0.4
0 8
10
0.8
0.94
0.95
12
0.6 0.4 0.2 0 0.75
14
1
0.96
(d) Scenario 2; FOE
AN US
0.2
0.8
0.85
0.9
0.95
0.5
0.6
1
0.8
M
0.6 0.4 0.2
ED
(e) Scenario 3; AVC
0 4
5
6
7
PT
0
(f) Scenario 3; FOE
0.6 0.4 0.2
0.2
0.3
0.4
Frequency of exceeding
BMA_IMs
MCMC_IM1
MCMC_IM2
MCMC_IM4
MCMC_IM50.1
True Model
CE
MCMC_IM3
0
0.8
0 8
Average NO3-N concentration BAIPU_IMs 1
Cumulative frequency
Cumulative frequency
(b) Scenario 1: FOE
0 0.93
20
1
Cumulative frequency
0.8
CR IP T
Cumulative frequency
1
AC
Figure 8. True and simulated values of the two management variables under the three future loading scenarios in the synthetic modeling case. (a), (c) and (e) are the results of average nitrate concentration (AVC). (b), (d) and (f) are the results of frequency of exceeding (FOE) 6 mg N/L.
67
1
1
1
0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
(a) ESCO
0
CDF
0
0.1
0.2
0.3
0
1
1
0.8
0.8
0.6
0.6
0.4
0.4 (d) SOL_AWC
0.2 0.3
0.4
0.5
0.6
0.4 0.2
(e) CH_K1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.25
0.3
0
0.2
(h) ALPHA_BNK
0
0.005
0.01
0.4
0.8
0.8
0.8
0.6
0.6
0.4
0.4
ED
1
PT
0.2
0.2
0
0.6
0.8
1 0
0 1
0.8
1
(l) SDNCO
0 0
0.1
0.2
BAIPU_RMs
MCMC_RM1
MCMC_RM3 0
MCMC_RM4
AC
CE
0.4
0.6
0.2
(k) NPERCO
(j) GW_DELAY
(i) SHALLST
0
1
0.4
0.1
0.4
1
0.6
0.05
0.6
0 0.2
0.01
0.8
M
0
0.005
1
0.2
(g) CH_K2
(f) CH_N1
0
0
1
0.2
0.8
0 0.2
CDF
0.1
1
0.2
0
CDF
0.05
(c) SURLAG
0 0.15
CR IP T
0
0.2
(b) CN2
AN US
CDF
ACCEPTED MANUSCRIPT
0.3
0
0.5
1
MCMC_RM2 0.1
MCMC_RM5
Figure 9. Cumulative distribution functions (CDFs) of uncertain model parameters inferred by the BAIPU and individual MCMC analyses in the real-situation modeling case. (a)-(i) are for the twelve model parameters.
68
ACCEPTED MANUSCRIPT
(a) w1 0.2
0.4
0.8
(d) w4 0.15
0.00 0.0
0.25
0.00 0.0
0.3
0.6
(c) w3 0.3
0.0 0.00
0.25
0.50
0.30
Posterior density
Posterior density
0.30
(b) w2
0.5
(e) w5 Sequential BMA BAIPU
0.15
CR IP T
0.0 0.0
0.6
Posterior density
0.50
Posterior density
Posterior density
0.4
Real-situation case
0.00 1.0 0.0
0.5
1.0
AN US
Figure 10. Marginal probability distribution functions of model weights derived based on BAIPU (red bars) and the sequential BMA approach (grey bars) in the
real-situation modeling case. Red points and black circles indicate the respective
AC
CE
PT
ED
M
mean values. (a)-(e) correspond to the five models.
69
ACCEPTED MANUSCRIPT
30
95% conf. int.
-2 -12 Median
NO3-N (mg N/L)
25
50% conf. int.
(a) BAIPU
Observations
20 15 10 5
CR IP T
0
30
(b) Sequential BMA 20
15
AN US
10 5 0
7
(c) Predictive Mean
M
20
ED
15 10 5 0 0
PT
By sequential BMA (mg N/L)
25
5
10
15
By sequential BMA (mg N/L)
NO3-N (mg N/L)
25
(d) Predictive std. dev.
6 5 4 3 2
20
2
25
3
4
5
6
7
By BAIPU (mg N/L)
CE
By BAIPU (mg N/L)
AC
Figure 11. Results of the uncertainty analysis in the real-situation modeling case. (a) and (b) illustrate the uncertainty bands (50% and 95% confidence intervals) of nitrate concentration based on the BAIPU and the sequential BMA approach, respectively. (c) and (d) compare the means and standard deviations of the predictive distributions based on the two approaches. Each circle in the subplots (c) and (d) represents a simulation day. 70
ACCEPTED MANUSCRIPT
1
(a) Scenario 1; AVC
0.8
Cumulative frequency
0.6 0.4 0.2 0 10
11
12
13
0.6 0.4 0.2
Cumulative frequency
1
0.8 0.6 0.4
0
0.9
0.8 0.6 0.4
0.95
(d) Scenario 2; FOE
0.2
AN US
(c) Scenario 2; AVC
0.2
0
7
7.5
8
8.5
1
9
0.6
0.65
0.7
0.75
0.8
1
0.8
M
0.6 0.4
(e) Scenario 3; AVC
0 4
ED
0.2
4.5
Cumulative frequency
Cumulative frequency
0.8
0 0.85
14
1
Cumulative frequency
(b) Scenario 1; FOE
CR IP T
Cumulative frequency
1
(f) Scenario 3; FOE
0.8 0.6 0.4 0.2
0 5
0
PT
Average NO3-N concentration BAIPU_RMs 1
0
0.2
0.3
BMA_RMs
MCMC_RM1
MCMC_RM4
MCMC_RM5
MCMC_RM2
0.08
CE
MCMC_RM3 0
0.1
Frequency of exceeding
AC
Figure 12. Simulated values of the two management-relevant variables under the three future loading scenarios in the real-situation modeling case. (a), (c) and (e) are the results of average nitrate concentration (AVC). (b), (d) and (f) are the results of frequency of exceeding (FOE) 6 mg N/L.
71
Joint analysis of input and parametric uncertainties in watershed water quality modeling: a formal Bayesian approach Feng Han , Yi Zheng PII: DOI: Reference:
S0309-1708(17)30982-X 10.1016/j.advwatres.2018.04.006 ADWR 3126
To appear in:
Advances in Water Resources
Received date: Revised date: Accepted date:
30 October 2017 11 April 2018 11 April 2018
Please cite this article as: Feng Han , Yi Zheng , Joint analysis of input and parametric uncertainties in watershed water quality modeling: a formal Bayesian approach, Advances in Water Resources (2018), doi: 10.1016/j.advwatres.2018.04.006
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.
ACCEPTED MANUSCRIPT
Highlights
AC
CE
PT
ED
M
AN US
CR IP T
Joint analysis of input and parametric uncertainties by a rigorous Bayesian method Markov Chain Monte Carlo analysis tightly coupled with Bayesian Model Averaging An inverse modeling approach to joint inference of uncertain parameters and inputs Impact of parameter-input interaction on modeling uncertainty addressed
1
ACCEPTED MANUSCRIPT
Joint analysis of input and parametric uncertainties in watershed water quality modeling: a formal Bayesian approach
a
CR IP T
Feng Han a, b, c and Yi Zheng a, c *
School of Environmental Science and Engineering, Southern University of Science and Technology, Shenzhen 518055, Guangdong Province, PR China.
b
School of Water Resources and Hydropower Engineering, Wuhan University, Wuhan
c
AN US
430072, Hubei Province, PR China.
Guangdong Provincial Key Laboratory of Soil and Groundwater Pollution Control, Southern University of Science and Technology, Shenzhen 518055, Guangdong
ED
M
Province, PR China.
* Corresponding author:
PT
Dr. Yi Zheng, Professor
Room 506, Block 2, Wisdom Valley, Southern University of Science and Technology,
CE
Shenzhen 518055, Guangdong Province, PR China
AC
Email: [email protected] Work phone/fax: 86-755-88018030
2
ACCEPTED MANUSCRIPT
Abstract Significant Input uncertainty is a major source of error in watershed water quality (WWQ) modeling. It remains challenging to address the input uncertainty in a rigorous Bayesian framework. This study develops the Bayesian Analysis of Input and
CR IP T
Parametric Uncertainties (BAIPU), an approach for the joint analysis of input and parametric uncertainties through a tight coupling of Markov Chain Monte Carlo
(MCMC) analysis and Bayesian Model Averaging (BMA). The formal likelihood
AN US
function for this approach is derived considering a lag-1 autocorrelated,
heteroscedastic, and Skew Exponential Power (SEP) distributed error model. A series of numerical experiments were performed based on a synthetic nitrate pollution case
M
and on a real study case in the Newport Bay Watershed, California. The Soil and
ED
Water Assessment Tool (SWAT) and Differential Evolution Adaptive Metropolis (DREAM(ZS)) were used as the representative WWQ model and MCMC algorithm,
PT
respectively. The major findings include the following: 1) the BAIPU can be
CE
implemented and used to appropriately identify the uncertain parameters and
AC
characterize the predictive uncertainty; 2) the compensation effect between the input and parametric uncertainties can seriously mislead the modeling based management decisions, if the input uncertainty is not explicitly accounted for; 3) the BAIPU accounts for the interaction between the input and parametric uncertainties and therefore provides more accurate calibration and uncertainty results than a sequential analysis of the uncertainties; and 4) the BAIPU quantifies the credibility of different 3
ACCEPTED MANUSCRIPT
input assumptions on a statistical basis and can be implemented as an effective inverse modeling approach to the joint inference of parameters and inputs. Keywords: Water quality modeling; watershed; Bayesian Model Averaging; Markov
AC
CE
PT
ED
M
AN US
CR IP T
Chain Monte Carlo; input uncertainty; Bayesian calibration
4
ACCEPTED MANUSCRIPT
1. Introduction Watershed water quality (WWQ) models typically integrate distributed hydrological modeling with water quality modeling, and depict the interrelated watershed processes of water, heat, sediments and chemicals [1]. Some representative
CR IP T
WWQ models include the Soil and Water Assessment Tool (SWAT) [2, 3], Watershed Analysis Risk Management Framework (WARMF) [4, 5], Hydrological Simulation Program - Fortran (HSPF) [6], and Annualized Agricultural Non-Point Source
AN US
Pollution model (AnnAGNPS) [7, 8]. These models have been widely used to address a variety of water and environmental issues, including agricultural non-point source (NPS) pollution [9-13], water quality response to climate change [14, 15],
M
food-energy-water nexus (FEW) [16, 17], low impact development (LID) [18], and
ED
many others. However, simulation results of the WWQ models often have notable uncertainty [19, 20]. First, owing to current knowledge gaps, the existing models
PT
often represent complex water quality processes (e.g., soil erosion, heat exchange,
CE
pollutant transport and chemical reactions) in a highly simplified way, which leads to
AC
significant model structural errors [20-22]. Second, while the quantity and quality of the hydrological forcing data (e.g., precipitation and temperatures) have been greatly improved due to the development of remote sensing [23], water quality drivers (e.g., point source (PS) and NPS loadings) still lack good data with the required temporal and spatial resolutions [24], which represents a significant input uncertainty in the WWQ models. Third, WWQ models involve a larger number of model parameters, 5
ACCEPTED MANUSCRIPT
and the parametric uncertainty would be substantial when observational data for model calibration are scarce and/or inaccurate. The multiple sources of uncertainty significantly hinder the application of WWQ models in scientific research and management practices, and effective approaches to uncertainty analysis are highly
CR IP T
desired for the models.
Bayesian calibration using Markov chain Monte Carlo (MCMC) sampling [25, 26] is a formal approach for model calibration and uncertainty assessment and has been
AN US
applied widely in recent years [27-32]. In the approach, model residuals (i.e., the
difference between observations and corresponding model outputs) are described by a statistical residual error model, based on which the likelihood function can be derived
M
[33]. When the prior distribution is provided, the posterior parameter distribution can
ED
be inferred by adopting the Bayes’ rule. The approach has explicit assumptions and a rigorous mathematical basis, but its effectiveness is dependent on many issues. One
PT
issue is the appropriateness of the error model. For example, in hydrological modeling,
CE
model residuals are often correlated, heteroscedastic and not normally distributed [34].
AC
There are a series of studies which discuss the residual model [28, 30, 35-39]. Another important issue is the consideration of input uncertainty and model structural errors. A classic Bayesian calibration approach directly addresses the parametric uncertainty and leaves out the other two uncertainty sources, which may bias the inference of the posterior distributions. Our previous study shows that such degradation can be alleviated by including additional model responses that are less impacted by the input 6
ACCEPTED MANUSCRIPT
and model structural errors in a multiple-response Bayesian calibration (MRBC) framework [24]. Another solution is to jointly assess the input, structural and parametric uncertainties, although this remains challenging [40]. Hierarchical Bayesian approaches, such as Bayesian Total Error Analysis
CR IP T
(BATEA) [41], Integrated Bayesian Uncertainty Estimator (IBUNE) [42], Integrated Parameter Estimation and Uncertainty Analysis Tool (IPEAT) [43], and some others [44, 45] have been developed to quantify the model input uncertainty along with
AN US
model calibration. In the existing approaches, input errors are usually represented by latent variables (e.g., rainfall multipliers) with prescribed hyperdistributions. The hierarchical Bayesian inference addresses all uncertain model parameters and latent
M
variables, as well as all hyperparameters [46]. Although these approaches can separate
ED
the input and parametric uncertainties, their hierarchical inference framework can lead to high dimensional posterior distributions, and sampling of the posterior distributions
PT
can be computationally prohibitive, especially for complex WWQ models. More
CE
importantly, in WWQ modeling, it is very difficult to represent the input uncertainty
AC
of pollutant loadings in terms of latent variables. Bayesian Model Averaging (BMA) has been increasingly used to account for
model structural uncertainty [42, 47-53]. BMA derives an aggregated simulation by averaging the simulations of multiple models based on their posterior probabilities (also referred to as the model weights) [54, 55]. BMA has also been used for model ranking and selection, where the model weights should balance goodness of fit with 7
ACCEPTED MANUSCRIPT
complexity following the law of parsimony [56-58]. In hydrological and water quality modeling, it is a common situation for multiple plausible datasets or assumptions to be used to prepare a specific model input. For example, precipitation data can be simulated by a regional climate model, interpolated from observations at
CR IP T
meteorological stations, or estimated based on radar data. Land cover data can be interpreted from different types of satellite data, and different agencies may have
different estimates of pollutant loadings. In all these cases, the input uncertainty can
AN US
hardly be represented as continuous probability distributions, and needs to be
addressed in a discrete form. BMA appears to be an appropriate way to address such input uncertainty. For example, Strauch et al. [59] established four SWAT models
M
with four different sets of precipitation data, and then performed BMA to account for
ED
the precipitation uncertainty in streamflow modeling. Some studies have attempted the coupling of BMA with Bayesian calibration in a sequential manner, that is,
PT
performing Bayesian calibration for individual models first and then achieving an
CE
aggregated simulation using BMA. For example, Zeng et al. [52] sequentially
AC
performed MCMC and BMA analyses to assess groundwater conceptualization uncertainty. However, this sequential approach ignores the interaction between the parametric uncertainty and conceptual uncertainty, which may lead to inaccurate results under certain circumstances. This study aims to develop an approach to investigating parametric and input uncertainties simultaneously through a tight coupling of MCMC and BMA. The 8
ACCEPTED MANUSCRIPT
approach is specifically designed for input uncertainty in a discrete form, as exemplified above. A previously developed SWAT model for the nitrate pollution in the Newport Bay Watershed (Southern California, U.S.) was used as the study case, and a state-of-art MCMC algorithm, the Differential Evolution Adaptive Metropolis
CR IP T
(DREAM(ZS)) [44, 60, 61], was employed. A series of numerical experiments were
designed and implemented to answer the following key questions: 1) whether the new approach can effectively address the discrete form input uncertainty in MCMC-based
AN US
Bayesian calibration; 2) how the discrete form input uncertainty in pollutant loadings would impact the Bayesian calibration results and potential water quality management decisions based on the results; and 3) what information the joint analysis of input and
M
parametric uncertainties can provide to measure the relative importance of the
ED
competing model input assumptions.
2. Coupling MCMC and BMA
PT
2.1. MCMC uncertainty analysis
CE
Markov Chain Monte Carlo (MCMC) uncertainty analysis represents a classic
AC
Bayesian calibration approach. In this approach, the residual errors, i.e., ε Z F , are described by a statistical error model, where F and Z are the model outputs and the observations, respectively. This lumped error accounts for all sources of uncertainty including model input, model structure, model parameters and observations for model calibration [36]. With the mathematical form of the error model specified, the likelihood function of model parameters can be derived, and the posterior 9
ACCEPTED MANUSCRIPT
distributions of the parameters can be obtained based on Bayes’ rule, as
p(θ, φ | X, Z) p(θ) p(φ) L(θ, φ, X, Z)
(1)
where θ represents model parameters, φ denotes the parameters of the error model, X denotes the model input, p(θ) and p(φ) are the prior distributions, and
CR IP T
L(θ, φ, X, Z) is the likelihood function. In an MCMC algorithm, samples are
obtained by constructing Markov chains with desired posterior distributions as their
stationary distributions [25, 26]. The samples obtained by the Markov chains are used
AN US
to derive the posterior distribution of uncertain parameters, p(θ, φ | X, Z) , as well as
the predictive distribution, p(Y | θ, φ, Z) , where Y denotes a vector of quantities to be predicted. Explanation of the approach to derive the predictive distribution can be
ED
hydrological modeling [27-32].
M
found elsewhere [62]. This paradigm has been widely used in uncertainty analysis for
This study considers an autocorrelated, heteroscedastic, and non-Gaussian error
PT
model [24, 62]. The autocorrelation of residuals is depicted by a first order
CE
autoregressive model (i.e., AR(1)), which involves a parameter, lag-1 autoregressive parameter ϕ. The heteroscedasticity is reflected by a linear equation, t b1 b2 ft ,
AC
where t is the estimated standard deviation, f t is the model output at time step t, and b1 and b2 are the two parameters. The Skew Exponential Power (SEP) distribution [36] is employed to represent two typical non-Gaussian effects, i.e., skewness and kurtosis, where the skewness parameter ξ and the kurtosis parameter β are introduced. Thus, the error model has five parameters in this case, that is, 10
ACCEPTED MANUSCRIPT
φ [ , b1 , b2 , , ] . The log-likelihood function can be formulated as logL(θ, φ, Z) J log
2
where J is the data length of Z; a , j
J
J
j 1
j 1
log j c a , j 1
sign ( a j )
2/(1 )
(2)
( a j ) , in which a j stands
CR IP T
for the independent and identically distributed random errors with a zero mean and unit standard deviation; and , , and c are functions of ξ and β. More details about the error model and likelihood function can be found elsewhere [24, 36, 38]. Note
AN US
that this highly complicated likelihood function may lead to an undesired
computational burden on the uncertainty analysis. To simplify the error model, as well as the corresponding likelihood function, a trail-and-error approach can be adopted to
M
examine the fulfillment of the statistical assumptions a posteriori. One can gradually increase the complexity of the error model, as was done in previous studies [24, 62].
ED
Differential Evolution Adaptive Metropolis (DREAM(ZS)) [44, 60, 61] is
PT
employed in this study as the MCMC algorithm. It is a population evolution MCMC algorithm which runs multiple Markov chains in parallel to facilitate the global
CE
exploration of the parameter space. It automatically tunes the scale and orientation of
AC
the proposal distribution. In DREAM(ZS), a candidate point can be generated by using a differential evolution with a randomized subspace sampling strategy, which can be described as follows: θ*,i θi (1d ed ) ( , d ) u r1 ( j ) u r2 ( n ) є d n 1 j 1
where θ*,i is the candidate point, θi is the current point of the i-th Markov chain, ed 11
(3)
ACCEPTED MANUSCRIPT
and ϵd are two d-dimensional random vectors, the coefficient γ controls the jumping size, and {ur1 ( j ) , ur2 ( n) | j, n 1,..., } represents δ pairs of previously sampled points ( r1 ( j ) r2 (n) ). In high dimension problems, it is often not optimal to sample all d dimensions simultaneously. In DREAM(ZS), only d’ dimensions of θ*,i, which are
CR IP T
randomly selected with a ―crossover probability‖, are updated in each sampling step, and the other dimensions remain the same as those of θi. To increase diversity, the
snooker update is included, which is an additional way to generate candidate points.
AN US
The way to propose a candidate point is described as follows:
θ*,i θi s θi u D(u R1 , u R 2 )
(4)
where γs is a random variable, u, uR1 and uR2 are three previously sampled points, and
M
D is a function of uR1 and uR2. Both D and γs are used to determine the distance of the
ED
jumping. More details of the snooker update can be found in ter Braak and Vrugt [63]. In DREAM(ZS), both differential evolution and snooker update are used to generate
PT
the candidate points. In this study, we updated the Markov chains in a mix of 10%
CE
snooker update and 90% differential evolution. The R statistics proposed by Gelman
AC
and Rubin [64] are used to assess the convergence of the Markov chains. The Markov chains are deemed to achieve the convergence when the R values of all dimensions (or parameters) are below 1.2, computed based on the last 50% points of the Markov chains.
2.2. Loose coupling: the sequential approach As mentioned in Section 1, MCMC-based Bayesian calibration can be coupled 12
ACCEPTED MANUSCRIPT
with BMA in a sequential manner to account for input uncertainty. Hereinafter, models with the same model structure but different inputs are considered as different models. When K plausible input datasets, X1, …, XK, are collected or assumed for a same model structure, K different models, M1, …, MK, can be established, respectively,
CR IP T
which represents the input uncertainty in a discrete form. MCMC uncertainty analysis can first be performed for each model to assess the parametric uncertainty. Next, the predictive distributions inferred based on the K models (i.e., p(Y | M k , Z) ) are
AN US
averaged into a total predictive distribution, representing both parametric and input
uncertainties. This two-step uncertainty analysis is referred to as the sequential BMA approach in this paper. Its mathematical details are provided below.
M
BMA provides a data-driven approach to statistically integrate the predictions of , yt ,
, yT denote the vector of quantities
ED
multiple models [54, 55]. Let Y y1 ,
to be predicted (i.e., the true values), and let Z denote the observations. In BMA, the
PT
predictive distribution of yt , based on the k-th model M k , can be deemed as the
CE
posterior distribution of yt conditional on Z, that is, p( yt | M k , Z) . Thus, the
AC
posterior distribution of yt conditional on the entire model ensemble can be written as
K
p( yt | Z) p( yt | M k , Z) p( M k | Z)
(5)
k 1
where p( M k | Z) stands for the posterior probability of M k and p( M k | Z) can be interpreted as the likelihood of M k being ―correct‖, given Z. Because the posterior probabilities of all the models have a sum of one, p( M k | Z) can also be 13
ACCEPTED MANUSCRIPT
interpreted as the weight of M k . Hereinafter, p( M k | Z) is denoted as wk and
w [w1 ,..., wK ] . The mean and variance of BMA predictions can be computed as K
E[ yt Z] wk E[ yt M k , Z]
(6)
k 1
K
K
Var[ yt Z] wk E[ yt M k , Z] E[ yt Z] wk Var[ yt M k , Z] 2
(7)
k 1
CR IP T
k 1
The predictive distribution p( yt | M k , Z) can be computed as
p( yt | M k , Z) p( yt | M k , θk , φk ) p(θk , φk | M k , Z)dθk dφk k
(8)
AN US
where θ k denotes the model parameters of M k ; φ k denotes the error parameters
associated with the residual error model, θk , φk k ; p( yt | M k , θk , φk ) represents the predictive distribution of yt given the specific values of θ k and φ k ; and
M
p(θk , φk | M k , Z) represents the posterior distribution of θ k and φ k . The
ED
calculation of the integral in Equation (8) is computationally intensive. Alternatively, the integral can be approximated using the posterior samples generated in the MCMC
CE
PT
analysis, as [26]
pˆ ( yt | M k , Z)
1 Nk
Nk
p( y n 1
t
| M k , θ(kn ) , φ(kn ) )
(9)
AC
where θ(kn ) and φ (kn ) denote the n-th sample of model parameters and error parameters, respectively, jointly generated by the Markov chain, and N k is the sample size. Assuming independence among the residual errors, the log-likelihood function of the model weights can be derived as
14
ACCEPTED MANUSCRIPT
K 1 logL(w, Z) log wk j 1 k 1 N k J
(n) (n) p ( z | M , θ , φ ) j k k k n 1 Nk
(10)
Note that the assumption of independence may not be valid. For example, when lag-1 autocorrelation of the residual errors exists, z j should be conditional on z j 1 . In
CR IP T
such cases, p( z j | M k , θ(kn ) , φ(kn) ) in Equation (10) should be replaced by p( z j | M k , θ(kn ) , φ(kn ) , z j 1 ) . Raftery et al. [47] suggested to estimate the weights using
the maximum likelihood estimation (MLE). For complex likelihood functions, the
AN US
expectation-maximization (EM) algorithm [65] can be employed to find a numerical solution for the MLE problem with a small computational cost. However, the EM method only leads to a set of deterministic weight estimates, and the results can be
M
unstable when the weights are sensitive to the observations. A more reliable way is to perform MCMC simulation [66], in which the means of the sampled weights can be
ED
used as the posterior probabilities of the respective models. There are also other
PT
options to determine the model weights such as the Bayesian model evidence method [58]. In this study, we adopted DREAM(ZS) to sample model weights in the sequential
CE
BMA approach. To guarantee a sum of one for the weights of a candidate point (i.e.,
AC
θ*,i), we modified the original operation of generating candidate points to the following steps: Step 1: Generate a jumping vector ∆ by Equation (3) or (4), with at least two
elements of ∆ being non-zero. Step 2: Normalize ∆ by subtracting the mean of the updated elements, and therefore the sum of the normalized ∆ is equal to zero. 15
ACCEPTED MANUSCRIPT
Step 3: Generate the candidate point θ*,i = θi + ∆. The sum of θ*,i is one, as the sum of θi is one. Step 4: If θ*,i is outside of the sampling domain, calculate a new point (θ**,i) using the ―returning‖ operation introduced below.
CR IP T
As θ*,i may be outside of sampling domain, a ―returning‖ operation was proposed, as illustrated in Figure 1. When the proposed jump from θi to θ*,i hits the boundary (i.e., θb is reached), the jump turns back along the original path and ends at θ**,i. In
AN US
this case, the total jumping distance, in its absolute value, remains the same. Because the sum of ∆ is zero, the sum of θ**,i remains one. Similarly, a jump from θ**,i with the same jumping vector ∆ will end at θi. It suggests that, with the ―returning‖ operation,
M
the conditional probability of a jump from θi to θ**,i is equal to that of the reverse
ED
jump (i.e., from θ**,i to θi). Thus, the detailed balance of DREAM(ZS) can be guaranteed by accepting the candidate point θ**,i (or θ*,i), based on the Metropolis
PT
acceptance probability [61].
CE
2.3. Tight coupling: BAIPU
AC
This study develops a new approach to jointly analyze the input and parametric
uncertainties by tightly coupling BMA with MCMC, which is named the Bayesian Analysis of Input and Parametric Uncertainties (BAIPU). Figure 2 illustrates the framework of BAIPU. It begins with the preparation of K sets of plausible input data, which result in K models, as defined in Section 2.2. Thus, given the same model structure, one realization of model parameters θ would lead to K sets of model outputs, 16
ACCEPTED MANUSCRIPT
F1, …, FK. In BAIPU, a shared residual error model with parameters φ is also assumed for the model ensemble. Hence, the multiple models share the same set of uncertain parameters (i.e., θ plus φ). The predictive distributions of yt derived based on different models can be further averaged by BMA, as K
,K
, θ, φ) wk (θ, φ) p( yt | M k , θ, φ) k 1
CR IP T
p( yt | M1,
(11)
where wk (θ, φ) represents the weight of the k-th model given specific values of θ and
AN US
φ, and p( yt | M k , θ, φ) is the predictive distribution of yt based on the k-th model
(i.e., the k-th input dataset). In Equation (11), the dependency of the model weights on θ and φ reflects the interaction between the input and parametric uncertainties, which is
M
ignored in the sequential BMA approach. Assuming independence among the residual errors, the log-likelihood function of the model parameters and error parameters can
ED
be defined as
(12)
PT
J K logL(θ, φ, w, Ζ) log wk (θ, φ) p( z j | M k , θ, φ) j 1 k 1
CE
If lag-1 autocorrelation among the residual errors exists, the probability of z j in
AC
Equation (12) should also be conditional on z j 1 , that is, p( z j | M k , θ, φ, z j 1 ) . The MCMC analysis is then performed to derive the posterior distributions of model parameters and predictive variables. As discussed in Section 2.1, there are different strategies to infer the model weights. The framework of BAIPU accommodates two (Figure 2). One strategy is to estimate the weights by the EM algorithm. When a realization of θ and φ is generated 17
ACCEPTED MANUSCRIPT
in one sampling step, the corresponding model weights wk (θ, φ) (k=1, ..., K) are then determined by maximizing this likelihood function using the EM optimization. The other strategy is to sample the model weights jointly with the model parameters using MCMC. To implement this strategy, the ―returning‖ operation introduced in
CR IP T
Section 2.2 needs to be incorporated in DREAM(ZS), such that the model weights in each candidate point add up to one. Because the MCMC sampling is more reliable
than the EM optimization, all the BAIPU results presented in this paper were based on
AN US
the MCMC sampling. For comparison, we also performed all the BAIPU analyses
using the EM approach. Section 4.1 briefly discusses the comparison, and readers are referred to the Auxiliary Materials for more details.
M
One difference between the BAIPU and sequential BMA approach is that the
ED
BAIPU derives a common posterior parameter distribution based on the entire set of plausible model inputs, while the sequential approach leads to multiple posterior
PT
parameter distributions with respect to different model inputs. Note that, the multiple
CE
posterior parameter distributions could be further averaged into a single posterior
AC
distribution by resampling the multiple distributions according to their posterior probabilities (i.e., the model weights). However, although such a posterior distribution is relevant to the parameter calibration, it cannot be directly assigned to each model for making predictions, which is an advantage of BAIPU. In addition, in the BAIPU, because a common posterior parameter distribution can be obtained, the model weights solely reflect the model input uncertainty. However, in the sequential BMA 18
ACCEPTED MANUSCRIPT
approach, the input and the parametric uncertainties can compensate for each other in the Bayesian calibration; therefore, the weights are impacted by the parametric uncertainty. The error model introduced in Section 2.1 is considered for the model ensemble,
CR IP T
which involves five error parameters, that is, φ [ , b1 , b2 , , ] . Let kt yt f kt
denote the residual error of the k-th model (i.e., input X k ) at time step t, where f kt is the corresponding model output. As mentioned before, the standard deviation of kt
AN US
is estimated by kt b1 b2 f kt , and hence the standardized residuals can be
calculated as kt kt kt . The AR(1) model is applied to the standardized residuals to account for the lag-1 autocorrelation, that is, kt k ,t 1 okt . It is a standard
1 2 should be a unit. Similarly, we assume that akt follows the SEP
ED
of akt okt
M
homoscedastic AR(1) process and variance of okt equals 1 2 . Thus, the variance
distribution. Hence, given the specific values of θ and φ, the predictive probability of
PT
yt based on the k-th model can be computed as
CE
p( yt | M k , θ, φ, yt 1 ) sign ( akt )
1 2 kt
exp c a ,kt 1
2 /(1 )
(13)
( akt ) , and , , and c are functions of
AC
where a ,kt
2
1
and [36]. According to Equation (12), the log-likelihood function for this
residual error model can be derived as, K J 2 1 logL(θ, φ, w, Ζ) log wk (θ, φ) exp c a ,kj 1 k 1 j 1 1 2 kj
2/ (1 )
(14)
This likelihood function is quite general and flexible. It relaxes the commonly used 19
ACCEPTED MANUSCRIPT
independent, homoscedastic and Gaussian assumption about the residual errors and hence is more realistic for environmental modeling. As was revealed by Han and Zheng [24], when model input errors are very significant, the inferred posterior distributions of the model parameters would be prone to notable bias. This likelihood
CR IP T
function, which is based on multiple input datasets rather than one single set of input data, is expected to be more accurate and reliable. It is worth re-emphasizing that
BAIPU was specifically designed to address discrete form input uncertainty. If the
AN US
input uncertainty could be appropriately characterized as continuous probability distributions, the BMA step in BAIPU may not be an efficient solution.
In BAIPU, the sampler, DREAM(ZS), would execute K times more model runs
M
than in a traditional MCMC uncertainty analysis, if the chains have the same length.
ED
In this study, the computational cost was largely reduced by parallelizing the sampling process, that is, making DREAM(ZS) execute Nchain*K model runs (Nchain denotes the
PT
number of Markov chains) simultaneously in parallel. Similar parallel calculation was
CE
adopted for the sequential BMA as well.
AC
3. Data and Methods 3.1. Study area The Newport Bay Watershed (NBW) is in Orange County, southern California. It
is a highly urbanized watershed with an area of approximately 400 km2 (Figure 3). As of 2001, approximately 70% of the watershed was residential, commercial and industrial areas, and agricultural and orchard areas accounted for approximately 8%. 20
ACCEPTED MANUSCRIPT
The NBW has a typical Mediterranean climate with warm, dry summers and mild, wet winters. The annual rainfall is approximately 330 mm, and occurs mostly between November and April. Nutrient pollution was a significant water quality issue in this area in the 1990s. Excessive loadings of nutrients were delivered into the Newport
CR IP T
Bay through the San Diego Creek (see Figure 3), which had caused the eutrophication of the bay and created the recreational and aesthetic nuisance [67]. Major sources of
nutrients included fertilizer application (non-point source) and commercial nurseries
AN US
(point sources). The NBW and San Diego Creek were included in the 1998 California 303(d) list of impaired waters for Total Maximum Daily Loads (TMDL) actions, with nutrients as a major pollutant stressor [68]. More details about this watershed can be
ED
3.2. Nitrate pollution modeling
M
found in our previous studies [62, 69].
SWAT is a representative WWQ model which is semi-distributed and runs at a
PT
daily time step [3]. A SWAT model of the nitrate pollution in the NBW has been
CE
previously developed [24, 62]. In brief, the watershed was first delineated into 11
AC
sub-basins (see Figure 3), with a total of 58 hydrologic response units (HRUs). The 90-m Digital Elevation Model (DEM) of the study area was acquired from the U.S. Geological Survey (USGS) to delineate the modeling domain. The 30-m land use raster from the National Land Cover Database 2001 (NLCD 2001) and the soil type map from the U.S. State Soil Geographic (STATSGO) database were referred to in defining the HRUs. Meteorological data from 1998 to 2004 were obtained from U. S. 21
ACCEPTED MANUSCRIPT
National Climatic Data Center (NCDC) and used to drive the model simulations. The model was run from July 1998 to June 2004, with the first two years being the warm-up period and the last four years being the calibration period. Observations of daily flow and weekly nitrate concentration at the watershed outlet (Figure 3a) from
Control Board reports (available at
CR IP T
July 2000 to June 2004 were collected from the Santa Ana Regional Water Quality
http://prg.ocpublicworks.com/DocmgmtInternet/Search.aspx, last accessed on
AN US
October 8, 2017) and were used for the model calibration.
In this study, nitrate was the water quality parameter of concern, and nitrogen source loadings were considered as the uncertain model inputs. In the NBW, fertilizer
M
application on urban lawns and effluents from commercial nurseries were the two
ED
major nitrogen sources [67]. In the SWAT model, Bermuda was selected as the representative lawn grass and the nurseries were conceptualized into two point
PT
sources, A and B (see Figure 3a). As detailed information of the timing of fertilizer
CE
application is not available, the heat unit scheduling in SWAT was chosen, which
AC
allows the model to schedule operations based on temperature. To simulate the stormwater runoff and nitrate pollution in urban areas, SWAT offers two options. One is to use a set of linear empirical equations developed by the USGS (referred to as the USGS option), and the other is to simulate the buildup and washoff processes (referred to as the buildup-washoff option). Although the buildup-washoff option is conceptually appealing, local data to parameterize the equations is often lacking. 22
ACCEPTED MANUSCRIPT
Hence, the USGS option was chosen in this study. In our previous study [24], the nitrate loadings at sources A and B were estimated to be 61 kg N/day and 31 kg N/day, respectively, and the fertilizer usage was assumed to be 240 kg N/ha/year. The PS loadings were assumed to be time-invariant because
CR IP T
no detailed temporal information was available. Inevitably, significant input
uncertainty exists in the source loadings, which would degrade the calibration of the
uncertain model parameters. Thus, this study considers multiple plausible PS and NPS
AN US
loadings instead, which will be further discussed in Section 3.4. 3.3. Multiple-response Bayesian calibration
To alleviate the impact of input and model structural errors, the streamflow and
M
nitrate concentration at the watershed outlet are calibrated simultaneously by applying
ED
the MRBC approach developed in our previous study [24]. In the context of the MRBC, different error models are required for different responses, based on which
PT
individual likelihood functions can be established. In this study, residual errors of
CE
streamflow and nitrate concentration are assumed to be independent. Hence, the
AC
combined likelihood function is simply the product of their likelihood functions [24, 70]:
Lmultiple LH (θ, φH , ZH ) LN (θ, φ N , Z N )
(15)
where LH and LN are the likelihood functions of streamflow and nitrate concentration, respectively. In this study, the likelihood functions of flow and nitrate responses are determined by Equation (2) and Equation (14) in the sequential BMA 23
ACCEPTED MANUSCRIPT
approach and BAIPU, respectively. 3.4. Numerical experiments In the real-situation case, five loading scenarios were assumed to account for the input uncertainty, which leads to five real-situation models (RMs) (see Table 1).
CR IP T
Among the five models, RM2 is the basic model that was built in our previous study
[24] and represents an average loading situation. RM1 has 20% lower PS loadings and 10% lower fertilizer usage, compared to RM2, while RM3 has 20% higher PS
AN US
loadings and 10% higher fertilizer usage. RM4 and RM5 are two models with
temporal variability in the PS loadings, while RM4 has a smaller seasonal fluctuation (see Figure S1 in the Auxiliary Materials). Note that the overall loadings in RM4 and
M
RM5 are the same as those in RM2. This ensemble of models effectively embodies a
ED
common situation encountered in WWQ modeling, that is, detailed and accurate loading data are often lacking, and the modeler may consider different assumptions.
PT
The selection of random parameters in the Bayesian calibration was the same as
CE
in our previous study [24]. In brief, 32 SWAT parameters that are physically relevant
AC
for flow and/or nitrate simulations were identified first. Next, 12 most important parameters (Table 2) were further determined after an initial sensitivity analysis on the 32 parameters using Morris screening [71]. Each of the 12 parameters ranks top 8 among the 32 initial parameters at least for one of the two responses (i.e., flow and nitrate). Full details of the selection procedure can be found in our previous study [24]. As most of the parameters have spatially distributed values, three strategies [30] were 24
ACCEPTED MANUSCRIPT
considered to vary them: varying the parameter value directly (Type I), adding a deviation value to the distributed parameters (Type II), and applying a multiplier to the distributed parameters (Type III). For parameters of Type II and III, the deviations and multipliers other than the original distributed parameters are treated as uncertain
CR IP T
model parameters in subsequent calibration approaches. All the uncertain parameters are assumed to be uniformly distributed, and their ranges are provided in Table 2.
Initialization of other non-random parameters is also the same as in the previous study
AN US
[24].
To investigate the applicability and performance of BAIPU, a synthetic modeling case was also constructed, in which a true model (denoted as TM) with true model
M
parameters, true model inputs and true model structure was hypothesized (Table 1).
ED
The true values of the 12 uncertain parameters (the last column in Table 2) were randomly picked, which can yield realistic simulation results. True nitrate loadings at
PT
sources A and B were assumed to be time-variant at the levels of 132 and 68 kg N/day,
CE
respectively (see Figure S2 in the Auxiliary Materials). The fertilizer usage in TM was
AC
assumed to be 240 kg N/ha/year and scheduled on specific dates. The buildup-washoff method was selected as the true model structure for the simulation of stormwater pollution from urban areas. Other model settings were the same as that in the real-situation case. True model responses of flow and nitrate concentration at the watershed outlet were achieved by running the true model with the above settings. In addition to the perfect model (i.e., TM), five imperfect models (IMs) with different 25
ACCEPTED MANUSCRIPT
source loadings were also created, denoted as IM1, IM2, IM3, IM4 and IM5 (see Table 1). The PS loadings were set to be time-invariant, which reflects the reality that the temporal information of PS loadings is often unavailable. Among the models, IM3 represents the average loading level (i.e., 150 kg/day and 75 kg/day at sources A and
CR IP T
B, respectively), and IM1, IM2, IM4 and IM5 represent the 66.7%, 83.3%, 116.7% and 133.3% of the average loading level, respectively. Similarly, different fertilizer
application rates were hypothesized for the five imperfect models (Table 1). The five
AN US
imperfect models collectively represent the hypothetical input uncertainty. To better mimic the reality, the model structural error was also introduced in the IMs. The USGS option was used instead of the buildup-washoff option. Daily flow
M
observations and weekly nitrate observations at the outlet were synthesized by
ED
corrupting the true model responses with hypothetic observational errors (see Figure S3 in the Auxiliary Materials). The hypothetic errors were generated from an error
PT
model with b1H 0.2 , b2H 0.05 , H 0.5 , H 3 , H 0.2 , b1N 0.2 , b2N 0.3 ,
CE
N 0.5 , N 3 , and N 0 . This procedure represents a common strategy to
AC
synthesize observational data [19, 36, 37]. As in the real-situation case, the artificial observations were generated for the period from July 2000 to June 2004, and all the models (i.e., the TM and IMs) were run from July 1998 to June 2004. In both the synthetic and real-situation cases, the BAIPU and the sequential BMA approach were carried out. In the BAIPU, the model weights were jointly sampled with the model and error parameters (i.e., the joint sampling method), and in the 26
ACCEPTED MANUSCRIPT
sequential BMA approach, the model weights were sampled by DREAM(ZS). Hereinafter, BAIPU_IMs denotes the BAIPU analysis on the five imperfect models, and MCMC_IM1, MCMC_IM2, MCMC_IM3, MCMC_IM4 and MCMC_IM5 denote the MCMC uncertainty analyses for the five individual models, respectively.
CR IP T
These individual MCMC analyses constitute the first step in the sequential BMA approach. The same terminology is applied to the real-situation case. The
computational efficiency, posterior distributions of the uncertain SWAT parameters
AN US
and the uncertainty results of the two approaches were evaluated and compared. It is
worth pointing out that each model (e.g., IM1) in the sequential BMA approach has its individual error model, while the multiple models in the BAIPU share the same error
M
model. In determining the proper error model, the complexity of the error models was
ED
gradually increased, from simple independent, homoscedastic and Gaussian type to the most complicated type, until the posterior checks confirmed that the residual
PT
errors are consistent with the error model assumptions. The error models for different
CE
calibration procedures are summarized in Table S1 in the Auxiliary Materials.
AC
Posterior justification of these error models is illustrated by Figures S4 to S8 in the Auxiliary Materials. Eventually, nine error parameters were included in the uncertainty analysis, i.e., φ [ H , b1H , b2H , H , H , b1N , b2N , N , N ] . In this study, the error parameters were jointly inferred with the uncertain SWAT parameters (i.e., θ) in the MCMC sampling. To further compare the sequential BMA approach and the BAIPU, three future 27
ACCEPTED MANUSCRIPT
loading scenarios with increased, steady and decreased loadings were proposed for the period from July 2004 to June 2008 (Table 3). The calibrated models were run to predict the nitrate concentration with the proposed future loading scenarios. To investigate how the calibration procedure would impact management decisions, two
CR IP T
management-relevant variables, average nitrate concentration (AVC) and frequency of exceeding (FOE) the water quality target (6 mg N/L), were evaluated based on the predictions. In the synthetic case, the true nitrate concentrations under different
AN US
loading scenarios were obtained by running the true model (i.e., TM), whereby the true values of AVC and FOE can be achieved. In the real-situation case, the true model was unknown, and hence the true AVCs and FOEs were not obtainable.
M
4. Results and Discussion
ED
4.1. Efficiency of DREAM(ZS) in BAIPU
Twelve Markov chains were run simultaneously in DREAM(ZS), and 12×K CPUs
PT
(2.67 GHz) were utilized in parallel to reduce the computing time, where K is the size
CE
of model ensemble. With the parallel computing, the cost of implementing
AC
DREAM(ZS) is affordable. For example, evaluating 200,000 samples in the BAIPU typically takes one day for the modeling cases. The sequential BMA approach has two stages of DREAM(ZS). The second stage is to sample the model weights, and the computational cost is only a few minutes because this stage involves no SWAT runs. To ensure that the posterior distributions and uncertainty results were reliable and representative, DREAM(ZS) explored 400,000 and 300,000 samples for the synthetic 28
ACCEPTED MANUSCRIPT
and real-situation cases, respectively. Using the fixed sample numbers can rule out the effect of computational cost when interpreting the analysis results of different models. The numbers are large enough to ensure the convergence of Markov chains in most of the experiments.
CR IP T
Table 4 shows the number of samples for Markov chains to achieve convergence (NC) and the acceptance rate (AR) of DREAM(ZS) in different numerical experiments. As indicated by the NC values, the convergence speed of Markov chains in BAIPU
AN US
was comparable to that in the MCMC uncertainty analyses in the real-situation case, but was moderately slower in the synthetic case. One obvious reason for the speed
decrease is that BAIPU leads to a higher dimensional posterior distribution. AR is the
M
probability that one candidate point is accepted in stationarity. It can be estimated as
ED
T'/T, where T and T' denote the numbers of all and the accepted candidate points in the second halves of the Markov chains, respectively. A higher AR simply means that it is
PT
easier to find an acceptable candidate point with relatively high posterior probability.
CE
As Table 4 shows, in the two-step sequential BMA, DREAM(ZS) obtained one order of
AC
magnitude higher AR values in estimating the model weights than in deriving the posterior parameter distributions. This is because the likelihood function of the model and error parameters (Equation (2)) is much more complex than that of the model weights (Equation (10)). Another interesting finding is that the two BAIPU analyses yielded several-fold higher AR than the respective MCMC uncertainty analyses of individual models. This is probably because the multi-model ensemble allows a 29
ACCEPTED MANUSCRIPT
broader range of predictions than a single model, making it easier to find an acceptable candidate point. Overall, the results in Table 4 indicate that the sampling efficiency of DREAM(ZS) in BAIPU was not significantly reduced by the integration with BMA.
CR IP T
As mentioned in Section 2.3, we also tested the EM strategy to infer the model weights in the BAIPU. It has been found that the EM strategy led to a faster
convergence of Markov chains (NC is 363,120 and 105,120 in the synthetic and
AN US
real-situation cases, respectively) but lower acceptance rates (AR is 2.75% and 7.14% in the synthetic and real-situation cases, respectively). The EM and joint sampling strategies yielded very similar posterior parameter distributions, uncertainty results,
M
and future predictions (refer to Figures S9 to S12 in the Auxiliary Materials). Yet the
ED
model weights are quite different with the two strategies (Figure S13 in the Auxiliary Materials). These results demonstrate that the joint sampling is more reliable for the
PT
weight estimation. However, if the convergence speed is an important concern, and
CE
accurate weight estimation is not a goal, the EM strategy can be desirable.
AC
4.2. The synthetic modeling case 4.2.1. Posterior parameter distributions Figure 4 illustrates the posterior parameter distributions derived by the BAIPU
(red lines) and the MCMC analyses (dotted lines) for the individual models in the synthetic case. The true values of the model parameters are indicated by the dashed vertical lines in the figure. All the parameter values were rescaled to [0, 1] by a linear 30
ACCEPTED MANUSCRIPT
min-max normalization. The posterior parameter distributions derived by the BAIPU (red lines) are relatively narrow and generally embrace the respective true parameter values, except for the initial depth of water in the shallow aquifer (SHALLST) and the denitrification threshold water content (SDNCO) parameter. The posterior distribution
CR IP T
of SHALLST is wide, but still embraces the true value; while the posterior
distribution of SDNCO is far off the true value. This is probably because the
uncertainty of this parameter compensated the model structure error not explicitly
AN US
accounted for in this study, and/or the interaction between the two parameters and the input assumptions in the MCMC sampling is significant. Figure 4 indicates that the model parameters have been appropriately identified by the BAIPU in general. For
M
most parameters, the posterior distributions based on the individual MCMC analyses
ED
are different from each other, as well as from the respective distribution derived by the BAIPU. The difference is most evident for IM5, the model with the most
PT
significant input error. This indicates that the calibration results can be highly
CE
dependent on the input assumption; therefore, an ensemble approach such as BAIPU
AC
is valuable for addressing the parametric uncertainty. Similar calibration results of the error parameters are presented in Figure S14 in the Auxiliary Materials. 4.2.2. Model weights Figure 5 compares the model weights based on the BAIPU and the sequential BMA approach. The weight results in the two approaches are similar. In general, the lower the input loading, the higher the model weights. Among the five models, the 31
ACCEPTED MANUSCRIPT
input scenario of IM2 is closest to the true one (see Table 1), but it only achieves the second highest weight in both approaches. IM1 receives the highest weight in both approaches, while its input scenario is the third closest to the true one. The hypothesized true PS loadings have skew distributions featured by a limited number
CR IP T
of loading pulses (see Figure S2 in the Auxiliary Materials); therefore, the mean values of the true nitrate loadings (e.g., 131.95 kg N/day at Source A) are much
higher than the respective medians (e.g., 88.58 kg N/day at Source A). This explains
AN US
why the inference favors the lowest input loading scenario IM1s, because the constant loadings in IM1 are much closer to the medians of the true PS loadings. Both BAIPU and sequential BMA enable the correlation analysis of the model
M
weights, which helps in understanding the tradeoff between different input scenarios.
ED
It was found that both approaches lead to very similar correlation results (see Table S2 in the Auxiliary Materials). BAIPU also allows a correlation analysis between the
PT
model weights and the parameters because they are jointly sampled by DREAM(ZS).
CE
As Table 5 shows, strong correlations exist in many parameter-weight pairs,
AC
especially for SHALLST and SDNCO, the two parameters that were not well identified. This implies that the significant interaction between the two parameters and the input assumptions in the MCMC sampling may be one reason for the weak identifiability of the two parameters. It also confirms the necessity of sampling model parameters and weights jointly in BAIPU.
32
ACCEPTED MANUSCRIPT
4.2.3. Uncertainty results Figure 6 displays the uncertainty results of the nitrate concentration achieved by the BAIPU and the sequential BMA approach. The medians and confidence intervals in Figure 6a and 6b are very similar. In both approaches, the seasonal variation is well
CR IP T
captured by the time series of the median concentrations. The Nash–Sutcliffe
coefficients of the median predictions by BAIPU and the sequential BMA approach
are 0.2932 and 0.2804, respectively. The observations are also adequately embraced
AN US
by the uncertainty bands. As demonstrated by Figures 6a and 6b, 48.33% and 50.72% of the observations are within the 50% confidence intervals, respectively, and 97.13% and 97.61% of the observations are within the 95% confidence intervals, respectively.
M
Figure 6c and 6d show that the means and standard deviations of the predictive
ED
distributions derived in the two approaches are very close. In the sequential BMA approach, the predictive means of individual models are first calculated, and then the
PT
ensemble predictive means are derived based on Equation (6), while in the BAIPU,
CE
the ensemble predictive means for each sampled realization of uncertainty parameters
AC
are first calculated according to Equation (6), and then the total predictive means are derived by averaging over all the sampled realizations of uncertainty parameters. The standard deviations illustrated in Figure 6d can be derived in a similar way based on Equation (7). The Quantile-Quantile (QQ) plots [46] of nitrate concentration in the synthetic case were also produced for the BAIPU, individual MCMC analyses and the 33
ACCEPTED MANUSCRIPT
sequential BMA approach (Figure 7a). The QQ plots can provide intuitive information on the consistency of predictive distribution with observations. In a QQ plot, if all the points are on the 45-degree line, the predictive distribution agrees perfectly with the observations. In Figure 7a, all the predictive distributions demonstrate adequate
CR IP T
consistency. In general, the BAIPU and the sequential approach perform better than individual MCMC analyses, because both approaches account for the input
uncertainty. Among the individual models, IM2 leads to the best consistency, because
AN US
its loading scenario is closest to the true one. It is worth pointing out that, although
the predictive distributions of different approaches are similar (Figure 6 and Figure 7a), the corresponding posterior parameter distributions significantly vary, which may
M
have a significant impact on future predictions and management decisions based on
ED
the stochastic modeling. This issue will be discussed further in Section 4.2.4. 4.2.4. Future predictions and management decisions
PT
The WWQ models are commonly used to assess hypothetical loading scenarios,
CE
based on which management decisions can be made with a strong scientific basis. To
AC
investigate the performance of the different approaches in predicting future water quality and therefore supporting decision-making, three future loading scenarios were assumed in the synthetic case, as discussed in Section 3.4 and summarized in Table 3. The future loading scenarios are free of input errors. Two management variables, average nitrate concentration (AVC) and frequency of exceeding (FOE) the water quality target (i.e., 6 mg N/L), were evaluated using the true parameter values, as well 34
ACCEPTED MANUSCRIPT
as the parameter values sampled in the Bayesian analyses. Figure 8 illustrates the cumulative frequency curves (CFC) of AVC and FOE derived by different approaches. Note that, in the BAIPU and MCMC analyses for individual models, the frequencies of AVC and FOE were calculated using the
CR IP T
parameter realizations in the second halves of the Markov chains. In the sequential BMA approach, the frequencies are weighted averages of the corresponding
frequencies derived by individual MCMC analyses. Thus, the range of corresponding
AN US
CFC depends on both the spread of the CFCs by the MCMC analyses and the model
weights. As Figure 8 shows, in general, the CFCs derived by BAIPU (solid red lines) and the sequential BMA approach (solid blue lines) are both close to the respective
M
true values. It indicates that the two approaches are relatively effective in reducing the
ED
bias caused by the input errors in the calibration stage. The CFCs by BAIPU are relatively narrower because BAIPU derived a common posterior parameter
PT
distribution. Moreover, the BAIPU demonstrates a better performance than the
CE
sequential approach in predicting the management variables, which reflects the fact
AC
that the interaction between input uncertainty and parametric uncertainty is considered. In contrast, the five individual models lead to distinct predictions, some of which are far away from the respective true values. This implies that the potential management decisions based on the Bayesian calibration of individual models would be seriously biased if the model inputs have large errors. Furthermore, the impact of the bias on potential decision-making could be dependent on the confidence level desired by the 35
ACCEPTED MANUSCRIPT
decision maker. Take FOE in Scenario 3 as an example (refer to Figure 8f). If a confidence level of 50% is considered, the value of FOE would be 0.28 according to IM5; therefore, the discrepancy with the true value is 0.14. However, given a confidence level of 95%, the discrepancy would instead be 0.07.
CR IP T
Figure 8 also shows that, IM1, which underestimates the observations in the calibration stage (i.e., MCMC_IM1), overpredicts AVC and FOE; IM5, which
overestimates the observations, now underpredicts AVC and FOE. This is because, in
AN US
the Bayesian calibration of individual models, the model parameters compensate for the effect of the input errors, such that the observations can be better matched. For example, underestimated PS loadings would lead to parameter values prone to
M
overestimate the water quality response (i.e., the case of IM1). Thus, when the model
ED
inputs are free of error, the compensation in the calibration stage causes an ―overreaction‖ of the model at the prediction stage. Similar overreaction behaviors,
PT
although not as significant, have also been observed for other individual models.
CE
Interestingly, IM3 has the best performance for predicting the two management
AC
variables, instead of IM2 (with the smallest input errors), which implies the compensation effect is complex and nonlinear. All the above results suggest that a Bayesian approach accounting for multiple
plausible input assumptions can provide a valuable mechanism for alleviating the impact of input uncertainty and therefore enhancing the reliability of decision-making in the context of real-world management. In addition, with the interaction between the 36
ACCEPTED MANUSCRIPT
input and parametric uncertainties systematically accounted for, the BAIPU has an improved performance in predicting management variables under the future loading scenarios compared to the straightforward sequential BMA approach. It is also worth emphasizing that one needs to carefully make input assumptions based on his/her best
CR IP T
knowledge when implementing the two approaches. This may require extra work such as data collection and compiling, field survey, literature review, etc. In theory, it is not a problem for BAIPU to accommodate many inappropriate input assumptions, but in
4.3. The real-situation modeling case
AN US
practice it may significantly increase the computational burden.
The performances of BAIPU and sequential BMA were also examined using the
M
real-situation modeling case (Table 1). Figure 9 shows the posterior parameter
ED
distributions derived by the BAIPU (red lines) and the MCMC analyses for individual models (dotted lines). Compared to the synthetic case (Figure 4), the difference
PT
between the BAIPU and the individual MCMC analyses is relatively small, which
CE
implies that the influence of input uncertainty on the parameter calibration is not as
AC
significant as that in the synthetic case. SHALLST and SDNCO are the two parameters with the most significant difference in the posterior distributions, which is also observed in Figure 4. In addition, the posterior distributions of the error parameters are presented in Figure S15 in the Auxiliary Materials. Figure 10 illustrates the model weights inferred by the BAIPU and the sequential BMA approach in this case. The overall weight patterns derived by the two 37
ACCEPTED MANUSCRIPT
approaches are similar. As introduced in Section 3.2, RM2 represents an average estimation of the nitrogen loadings based on the available data, while RM1 and RM3 represent a lower loading scenario and a higher loading scenario, respectively. Among the first three models (Figures 10a to 10c), RM1 receives the largest weights while
CR IP T
RM2 and RM3 receive the second and third largest weights in both approaches. It
suggests that the original loadings in RM2 has a large probability of overestimation.
On the other hand, in RM4 and RM5, the PS loadings are assumed to have small and
AN US
large seasonal fluctuations, respectively. The potential seasonal fluctuation could be
due to the temporal variability of tailwater discharges from the commercial nurseries. As shown in Figures 10d and 10e, both RM4 and RM5 were assigned significant
M
weights, which provides a unique evidence of the existence of seasonal variability in
ED
the PS loadings.
Table 6 shows the correlations between the model weights and parameters in
PT
BAIPU. Although statistically significant in general, the correlations in the
CE
real-situation case are not as strong as those in the synthetic case, indicating a weaker
AC
interaction between the input and parametric uncertainties. In addition, Table S3 in the Auxiliary Materials shows the correlations between different model weights. In this real-situation case, BAIPU and the sequential approaches led to slightly different results. The predictive uncertainty bands of the nitrate response are shown in Figure 11. The median predictions (black lines) in Figures 11a and 11b adequately match the 38
ACCEPTED MANUSCRIPT
observations. In the BAIPU and the sequential approach, the 50% confidence intervals encompasses 50.24% and 51.18% of the observations, respectively; the 95% confidence intervals encompasses 96.21% and 96.68% of the observations, respectively. Figures 11c and 11d further compare the predictive means and standard
CR IP T
deviations derived by the two approaches. In general, the BAIPU produces the same predictive means as the sequential approach but smaller standard deviations. This
indicates that the BAIPU reduces the predictive uncertainty in this case, compared to
AN US
the sequential approach.
The uncertainty results presented in Figure 11 and Figure 7b suggest that the two approaches, as well as those individual MCMC analyses, have a similar performance
M
during the Bayesian calibration. However, as revealed in the synthetic case, the
ED
difference could be more substantial in predicting the future water response. This real-situation modeling case also considers three future loading scenarios, as
PT
introduced in Table 3. Figure 12 shows that the BAIPU and the sequential BMA
CE
approach only have a slight difference in the frequency distributions, which indicates
AC
that the interaction between parameters and inputs in the Bayesian analysis is not as significant as that in the synthetic case. The difference between individual models are more significant, but not as significant as that in the synthetic case (Figure 8). This is because two time-variant loading scenarios are considered in this real-situation case, while the synthetic case considers five scenarios of constant loadings. Thus, in the WWQ modeling, adequately accounting for the temporal variability of source 39
ACCEPTED MANUSCRIPT
loadings is critical to the results of stochastic modeling. All the above results clearly demonstrate the applicability and usefulness of BAIPU in the context of real-world water quality management.
5. Conclusions
CR IP T
This study develops a new approach, BAIPU, for the joint analysis of input and parametric uncertainties through a tight coupling of MCMC and BMA. The formal
likelihood function for this approach is also derived considering a lag-1 autocorrelated,
AN US
heteroscedastic, and SEP distributed error model. Input uncertainty in a discrete form is specifically addressed, which is commonly encountered in environmental modeling for real-world problems. A series of numerical experiments were designed and
M
implemented based on a synthetic nitrate pollution case, as well as the real pollution
ED
case in the Newport Bay Watershed, California. The major conclusions include the following. First, with parallel computing, the BAIPU can be implemented with
PT
reasonable computational costs and appropriately identify the uncertain parameters
CE
and characterize the predictive uncertainty. Second, if significant input uncertainty is
AC
not accounted for, WWQ modeling may seriously mislead management decisions because the input uncertainty must be compensated for by the parametric uncertainty. Third, the BAIPU accounts for the interaction between the input and parametric uncertainties through joint analysis and is mathematically more rigorous than the sequential approach. Thus, in theory, it can provide more accurate calibration and uncertainty results than the sequential approach. Finally, the BAIPU quantifies the 40
ACCEPTED MANUSCRIPT
credibility (i.e., the model weights) of different input assumptions on a statistical basis, and it can be implemented as an inverse modeling approach for joint inference of both model parameters and inputs. The BAIPU was originally designed to address discrete form input uncertainty,
CR IP T
but future studies may expand the BAIPU to a joint analysis of model structural, input and parametric uncertainties. In addition, the BAIPU has a great potential of being applied in fields other than watershed water quality modeling. For example, in
AN US
basin-scale groundwater modeling, boundary conditions of subsurface inflow are
often highly uncertain [72], which could be addressed in the framework of BAIPU.
M
Acknowledgments
ED
This work was supported by the National Natural Science Foundation of China (No. 41622111; No. 91647201) and the China Postdoctoral Science Foundation
PT
(2017M612505). Additional support was provided by the Southern University of
CE
Science and Technology (No. G01296001). We thank Dr. Jasper A. Vrugt at the
AC
University of California, Irvine for kindly discussing with us on DREAM(ZS).
41
ACCEPTED MANUSCRIPT
References [1] Borah DK, M Bera. Watershed-scale hydrologic and nonpoint-source pollution models: Review of mathematical bases. American Society of Agricultural Engineers. 46 (2003) 1553-66, doi: 10.13031/2013.10398. [2] Arnold JG, R Srinivasan, RS Muttiah, JR Williams. Large area hydrologic modeling
Resources Association. 34 (1998) 73-89, doi: 10.1111/j.1752-1688.1998.tb05961.x.
CR IP T
and assessment part I: Model development. Journal of the American Water
[3] Neitsch SL, JG Arnold, JR Kiniry, JR Williams. Soil and water assessment tool
AN US
theoretical documentation version 2009. Texas Water Resources Institute & Texas A&M University 2011. http://swat.tamu.edu/media/99192/swat2009-theory.pdf. Accessed Oct. 14, 2017.
[4] Chen C, J Herr, L Weintraub. Decision support system for stakeholder involvement.
M
Journal of Environmental Engineering. 130 (2004) 714-21, doi: doi:10.1061/(ASCE)0733-9372(2004)130:6(714).
ED
[5] Zheng Y, W Wang, F Han, J Ping. Uncertainty assessment for watershed water quality modeling: A probabilistic collocation method based approach. Advances in
PT
Water Resources. 34 (2011) 887-98, doi: 10.1016/j.advwatres.2011.04.016. [6] Bicknell BR, JC Imhoff, JJL Kittle, JAS Donigian, RC Johanson. Hydrological
CE
simulation program -- Fortran, user's manual for version 11. U.S. Environmental Protection Agency 1997.
AC
[7] Young RA, CA Onstad, DD Bosch, WP Anderson. AGNPS: A nonpoint-source pollution model for evaluating agricultural watersheds. Journal of Soil and Water Conservation. 44 (1989) 168-73.
[8] Bingner RL, FD Theurer, Y Yuan. AnnAGNPS technical processes: version 5.4. 2015. https://www.wcc.nrcs.usda.gov/ftpref/wntsc/H&H/AGNPS/downloads/ AnnAGNPS_Technical_Documentation.pdf. Accessed Oct. 14, 2017. [9] Jeon J-H, CG Yoon, AS Donigian, K-W Jung. Development of the HSPF-Paddy 42
ACCEPTED MANUSCRIPT
model to estimate watershed pollutant loads in paddy farming regions. Agricultural Water Management. 90 (2007) 75-86, doi: 10.1016/j.agwat.2007.02.006. [10] Huang J, H Hong. Comparative study of two models to simulate diffuse nitrogen and phosphorus pollution in a medium-sized watershed, southeast China. Estuarine, Coastal and Shelf Science. 86 (2010) 387-94, doi: 10.1016/j.ecss.2009.04.003.
CR IP T
[11] Pease LM, P Oduor, G Padmanabhan. Estimating sediment, nitrogen, and phosphorous loads from the Pipestem Creek watershed, North Dakota, using AnnAGNPS. Computers & Geosciences. 36 (2010) 282-91, doi: 10.1016/j.cageo.2009.07.004.
[12] Volk M, D Bosch, V Nangia, B Narasimhan. SWAT: Agricultural water and
AN US
nonpoint source pollution management at a watershed scale. Agricultural Water Management. 175 (2016) 1-3, doi: 10.1016/j.agwat.2016.06.013. [13] Volk M, D Bosch, V Nangia, B Narasimhan. SWAT: Agricultural water and
nonpoint source pollution management at a watershed scale—Part II. Agricultural
M
Water Management. 180 (2017) 191-3, doi: 10.1016/j.agwat.2016.09.029.
ED
[14] Cho J, C Oh, J Choi, Y Cho. Climate change impacts on agricultural non-point source pollution with consideration of uncertainty in CMIP5. Irrigation and
PT
Drainage. 65 (2016) 209-20, doi: 10.1002/ird.2036. [15] Nguyen HH, F Recknagel, W Meyer, J Frizenschaf, MK Shrestha. Modelling the
CE
impacts of altered management practices, land use and climate changes on the water quality of the Millbrook catchment-reservoir system in South Australia.
AC
Journal of Environmental Management. 202 (2017) 1-11, doi: 10.1016/j.jenvman.2017.07.014.
[16] Karabulut A, BN Egoh, D Lanzanova, B Grizzetti, G Bidoglio, L Pagliero, et al. Mapping water provisioning services to support the ecosystem–water–food– energy nexus in the Danube river basin. Ecosystem Services. 17 (2016) 278-92, doi: 10.1016/j.ecoser.2015.08.002. [17] Psomas A, V Dagalaki, Y Panagopoulos, D Konsta, M Mimikou. Sustainable 43
ACCEPTED MANUSCRIPT
agricultural water management in Pinios River basin using remote sensing and hydrologic modeling. Procedia Engineering. 162 (2016) 277-83, doi: 10.1016/j.proeng.2016.11.059. [18] Her Y, J Jeong, J Arnold, L Gosselink, R Glick, F Jaber. A new framework for modeling decentralized low impact developments using Soil and Water
10.1016/j.envsoft.2017.06.005.
CR IP T
Assessment Tool. Environmental Modelling & Software. 96 (2017) 305-22, doi:
[19] Zheng Y, AA Keller. Uncertainty assessment in watershed-scale water quality modeling and management: 1. Framework and application of generalized
likelihood uncertainty estimation (GLUE) approach. Water Resources Research.
AN US
43 (2007), doi: 10.1029/2006wr005345.
[20] Rode M, G Arhonditsis, D Balin, T Kebede, V Krysanova, A van Griensven, et al. New challenges in integrated water quality modelling. Hydrological Processes. 24 (2010) 3447-61, doi: 10.1002/hyp.7766.
M
[21] Beck MB. Water quality modeling: A review of the analysis of uncertainty. Water
ED
Resources Research. 23 (1987) 1393-442, doi: 10.1029/WR023i008p01393. [22] van Griensven A, T Meixner. Methods to quantify and identify the sources of
PT
uncertainty for river basin water quality models. Water Science & Technology. 53 (2006) 51-9, doi: 10.2166/wst.2006.007.
CE
[23] Knoche M, C Fischer, E Pohl, P Krause, R Merz. Combined uncertainty of hydrological model complexity and satellite-based forcing data evaluated in two
AC
data-scarce semi-arid catchments in Ethiopia. Journal of Hydrology. 519 (2014) 2049-66, doi: 10.1016/j.jhydrol.2014.10.003.
[24] Han F, Y Zheng. Multiple-response Bayesian calibration of watershed water quality models with significant input and model structure errors. Advances in Water Resources. 88 (2015) 109-23, doi: 10.1016/j.advwatres.2015.12.007. [25] Gelman A, JB Carlin, HS Stern, DB Rubin. Bayesian data analysis. Chapman & Hall/CRC, United States of America, 2003. 44
ACCEPTED MANUSCRIPT
[26] Robert CP, G Casella. Monte Carlo statistical methods. Springer, United States of America, 2004. [27] Kuczera G, E Parent. Monte Carlo assessment of parameter uncertainty in conceptual catchment models: the Metropolis algorithm. Journal of Hydrology. 211 (1998) 69-85, doi: 10.1016/s0022-1694(98)00198-x.
CR IP T
[28] Bates BC, EP Campbell. A Markov Chain Monte Carlo scheme for parameter estimation and inference in conceptual rainfall-runoff modeling. Water Resources Research. 37 (2001) 937-47, doi: 10.1029/2000wr900363.
[29] Marshall L, D Nott, A Sharma. A comparative study of Markov chain Monte Carlo methods for conceptual rainfall-runoff modeling. Water Resources Research. 40
AN US
(2004), doi: 10.1029/2003wr002378.
[30] Yang J, P Reichert, KC Abbaspour, H Yang. Hydrological modelling of the Chaohe basin in China: Statistical model formulation and Bayesian inference. Journal of Hydrology. 340 (2007) 167-82, doi: 10.1016/j.jhydrol.2007.04.006.
M
[31] Smith TJ, LA Marshall. Bayesian methods in hydrologic modeling: A study of
ED
recent advancements in Markov chain Monte Carlo techniques. Water Resources Research. 44 (2008), doi: 10.1029/2007wr006705.
PT
[32] Vrugt JA, CJF ter Braak, HV Gupta, BA Robinson. Equifinality of formal (DREAM) and informal (GLUE) Bayesian approaches in hydrologic modeling?
CE
Stochastic Environmental Research and Risk Assessment. 23 (2008) 1011-26, doi: 10.1007/s00477-008-0274-y.
AC
[33] Box GEP, GC Tiao. Bayesian inference in statistical analysis. Wiley, United States of America & Canada, 1992.
[34] Kuczera G. Improved parameter inference in catchment models: 1. Evaluating parameter uncertainty. Water Resources Research. 19 (1983) 1151-62, doi: 10.1029/WR019i005p01151. [35] Schaefli B, DB Talamba, A Musy. Quantifying hydrological modeling errors through a mixture of normal distributions. Journal of Hydrology. 332 (2007) 45
ACCEPTED MANUSCRIPT
303-15, doi: 10.1016/j.jhydrol.2006.07.005. [36] Schoups G, JA Vrugt. A formal likelihood function for parameter and predictive inference of hydrologic models with correlated, heteroscedastic, and non-Gaussian errors. Water Resources Research. 46 (2010), doi: 10.1029/2009wr008933. [37] Smith T, A Sharma, L Marshall, R Mehrotra, S Sisson. Development of a formal
CR IP T
likelihood function for improved Bayesian inference of ephemeral catchments. Water Resources Research. 46 (2010), doi: 10.1029/2010wr009514.
[38] Evin G, D Kavetski, M Thyer, G Kuczera. Pitfalls and improvements in the joint inference of heteroscedasticity and autocorrelation in hydrological model calibration. Water Resources Research. 49 (2013) 4518-24, doi:
AN US
10.1002/wrcr.20284.
[39] Honti M, C Stamm, P Reichert. Integrated uncertainty assessment of discharge predictions with a statistical error model. Water Resources Research. 49 (2013) 4866-84, doi: 10.1002/wrcr.20374.
M
[40] Renard B, D Kavetski, G Kuczera, M Thyer, SW Franks. Understanding predictive
ED
uncertainty in hydrologic modeling: The challenge of identifying input and structural errors. Water Resources Research. 46 (2010), doi:
PT
10.1029/2009wr008328.
[41] Kavetski D, G Kuczera, SW Franks. Bayesian analysis of input uncertainty in
CE
hydrological modeling: 1. Theory. Water Resources Research. 42 (2006), doi: 10.1029/2005wr004368.
AC
[42] Ajami NK, Q Duan, S Sorooshian. An integrated hydrologic Bayesian multimodel combination framework: Confronting input, parameter, and model structural uncertainty in hydrologic prediction. Water Resources Research. 43 (2007), doi: 10.1029/2005wr004745. [43] Yen H, X Wang, DG Fontane, RD Harmel, M Arabi. A framework for propagation of uncertainty contributed by parameterization, input data, model structure, and calibration/validation data in watershed modeling. Environmental Modelling & 46
ACCEPTED MANUSCRIPT
Software. 54 (2014) 211-21, doi: 10.1016/j.envsoft.2014.01.004. [44] Vrugt JA, CJF ter Braak, MP Clark, JM Hyman, BA Robinson. Treatment of input uncertainty in hydrologic modeling: Doing hydrology backward with Markov chain Monte Carlo simulation. Water Resources Research. 44 (2008), doi: 10.1029/2007wr006720.
CR IP T
[45] Li M, D Yang, J Chen, SS Hubbard. Calibration of a distributed flood forecasting model with input uncertainty using a Bayesian framework. Water Resources Research. 48 (2012), doi: 10.1029/2010wr010062.
[46] Thyer M, B Renard, D Kavetski, G Kuczera, SW Franks, S Srikanthan. Critical evaluation of parameter consistency and predictive uncertainty in hydrological
AN US
modeling: A case study using Bayesian total error analysis. Water Resources Research. 45 (2009), doi: 10.1029/2008wr006825.
[47] Raftery AE, T Gneiting, F Balabdaoui, M Polakowski. Using Bayesian Model Averaging to calibrate forecast ensembles. Monthly Weather Review. 133 (2005)
M
1155-74, doi: 10.1175/mwr2906.1.
ED
[48] Duan Q, NK Ajami, X Gao, S Sorooshian. Multi-model ensemble hydrologic prediction using Bayesian model averaging. Advances in Water Resources. 30
PT
(2007) 1371-86, doi: 10.1016/j.advwatres.2006.11.014. [49] Rojas R, L Feyen, A Dassargues. Conceptual model uncertainty in groundwater
CE
modeling: Combining generalized likelihood uncertainty estimation and Bayesian model averaging. Water Resources Research. 44 (2008), doi:
AC
10.1029/2008wr006908.
[50] Zhang X, R Srinivasan, D Bosch. Calibration and uncertainty analysis of the SWAT model using Genetic Algorithms and Bayesian Model Averaging. Journal of Hydrology. 374 (2009) 307-17, doi: 10.1016/j.jhydrol.2009.06.023. [51] Bastola S, C Murphy, J Sweeney. The role of hydrological modelling uncertainties in climate change impact assessments of Irish river catchments. Advances in Water Resources. 34 (2011) 562-76, doi: 10.1016/j.advwatres.2011.01.008. 47
ACCEPTED MANUSCRIPT
[52] Zeng X, J Wu, D Wang, X Zhu, Y Long. Assessing Bayesian model averaging uncertainty of groundwater modeling based on information entropy method. Journal of Hydrology. 538 (2016) 689-704, doi: 10.1016/j.jhydrol.2016.04.038. [53] Wellen C, GB Arhonditsis, T Long, D Boyd. Quantifying the uncertainty of nonpoint source attribution in distributed water quality models: A Bayesian
(2014) 3353-68, doi: 10.1016/j.jhydrol.2014.10.007.
CR IP T
assessment of SWAT’s sediment export predictions. Journal of Hydrology. 519
[54] Leamer EE. Specification searches: Ad hoc inference with nonexperimental data. Wiley, United States of America, 1978.
[55] Hoeting JA, D Madigan, AE Raftery, CT Volinsky. Bayesian model averaging: A
AN US
tutorial. Statistical Science. 14 (1999) 382-417, doi: 10.1214/ss/1009212519.
[56] Schoniger A, T Wohling, W Nowak. A statistical concept to assess the uncertainty in Bayesian model weights and its impact on model ranking. Water Resources Research. 51 (2015) 7524-46, doi: 10.1002/2015WR016918.
M
[57] Tsai FTC, X Li. Inverse groundwater modeling for hydraulic conductivity
ED
estimation using Bayesian model averaging and variance window. Water Resources Research. 44 (2008) W09434, doi: 10.1029/2007WR006576.
PT
[58] Volpi E, G Schoups, G Firmani, JA Vrugt. Sworn testimony of the model evidence: Gaussian Mixture Importance (GAME) sampling. Water Resources Research. 53
CE
(2017) 6133-58, doi: 10.1002/2016WR020167. [59] Strauch M, C Bernhofer, S Koide, M Volk, C Lorz, F Makeschin. Using
AC
precipitation data ensemble for uncertainty analysis in SWAT streamflow simulation. Journal of Hydrology. 414-415 (2012) 413-24, doi: 10.1016/j.jhydrol.2011.11.014.
[60] Vrugt JA, CJF ter Braak, CGH Diks, BA Robinson, JM Hyman, D Higdon. Accelerating Markov Chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling. Int J Nonlin Sci Num. 10 (2009) 273-90, doi: 10.1515/ijnsns.2009.10.3.273. 48
ACCEPTED MANUSCRIPT
[61] Laloy E, JA Vrugt. High-dimensional posterior exploration of hydrologic models using multiple-try DREAM(ZS) and high-performance computing. Water Resources Research. 48 (2012), doi: 10.1029/2011wr010608. [62] Zheng Y, F Han. Markov Chain Monte Carlo (MCMC) uncertainty analysis for watershed water quality modeling and management. Stochastic Environmental
CR IP T
Research and Risk Assessment. 30 (2015) 293-308, doi: 10.1007/s00477-015-1091-8.
[63] ter Braak CJF, JA Vrugt. Differential Evolution Markov Chain with snooker updater and fewer chains. Statistics and Computing. 18 (2008) 435-46, doi: 10.1007/s11222-008-9104-9.
AN US
[64] Gelman A, DB Rubin. Inference from iterative simulation using multiple
sequences. Statistical Science. 7 (1992) 457-72, doi: 10.1214/ss/1177011136. [65] Dempster AP, NM Laird, DB Rubin. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B
M
(Methodological)1977. pp. 1-38.
ED
[66] Vrugt JA, CGH Diks, MP Clark. Ensemble Bayesian model averaging using Markov Chain Monte Carlo sampling. Environmental Fluid Mechanics. 8 (2008)
PT
579-595, doi: 10.1007/s10652-008-9106-3. [67] U.S. Environmental Protection Agency. Total Maximum Daily Loads for Nutrients
CE
San Diego Creek and Newport Bay, California. 1998. http://www. epa.gov/owow/tmdl/examples/nutrients/ca_sdnbay.pdf. Accessed Sep. 10, 2014.
AC
[68] U.S. Environmental Protection Agency. 1998 California 303(d) list and TMDL priority schedule. 1999. http://www.krisweb.com/biblio/ ncc_swrcb_ncrwqcb_2001_wqcontrolplanapdx.pdf. Accessed Oct. 14, 2017.
[69] Zheng Y, AA Keller. Stochastic watershed water quality simulation for TMDL development - a case study in the Newport Bay watershed. Journal of the American Water Resources Association. 44 (2008) 1397-410, doi: 10.1111/j.1752-1688.2008.00232.x. 49
ACCEPTED MANUSCRIPT
[70] Balin Talamba D, E Parent, A Musy. Bayesian multiresponse calibration of TOPMODEL: Application to the Haute-Mentue catchment, Switzerland. Water Resources Research. 46 (2010), doi: 10.1029/2007wr006449. [71] Campolongo F, J Cariboni, A Saltelli. An effective screening design for sensitivity analysis of large models. Environmental Modelling & Software. 22 (2007)
CR IP T
1509-18, doi: 10.1016/j.envsoft.2006.10.004. [72] Wu B, Y Zheng, Y Tian, X Wu, Y Yao, F Han, et al. Systematic assessment of the uncertainty in integrated surface water-groundwater modeling based on the
probabilistic collocation method. Water Resources Research. 50 (2014) 5848-65,
AC
CE
PT
ED
M
AN US
doi: 10.1002/2014wr015366.
50
ACCEPTED MANUSCRIPT
Tables Table 1. Different SWAT models in the synthetic and real-situation modeling cases True model
Model 1
Model 2
Model 3
Model 4
Model 5
/
RM1
RM2
RM3
RM4
RM5
Point source A a
/
48.8 (C)
61 (C)
73.2 (C)
61 (SF)
61 (BF)
a
/
24.8 (C)
31 (C)
37.2 (C)
31 (SF)
31 (BF)
Fertilizer use b
/
216 (HU)
240 (HU)
264 (HU)
240 (HU)
240 (HU)
Stormwater model
/
USGS
USGS
USGS
USGS
USGS
Model name
Point source B
Synthetic case TM
IM1
IM2
IM3
IM4
IM5
Point source A
a
132 (D)
100 (C)
125 (C)
150 (C)
175 (C)
200 (C)
Point source B
a
68 (D)
50 (C)
62.5 (C)
75 (C)
87.5 (C)
100 (C)
240 (Date)
165 (HU)
206.3 (HU)
247.5 (HU)
288.8 (HU)
330 (HU)
buildup-washoff
USGS
USGS
USGS
USGS
USGS
Fertilizer use
b
Stormwater model a
AN US
Model ID
CR IP T
Real-situation case
The average PS loading of nitrate pollution (kg N/day). "D" means the true PS loading in the synthetic case is Dynamic or time-variant. "C" represents a Constant loading input. "SF" and Fluctuations, respectively.
b
M
"BF" in the real-situation case indicate artificial time-variant loading inputs with Small and Big Fertilization rate in the NBW (kg N/ha/year). "Date" means the fertilizer applications are
AC
CE
PT
ED
scheduled on specific dates, and "HU" indicates the heat unit scheduling approach.
51
ACCEPTED MANUSCRIPT
Table 2. Uncertain SWAT parameters considered in the Bayesian inference Definition
Type
Range
ESCO
Soil evaporation compensation factor (/)
I
[0.01, 1]
Synthetic True value 0.02312
CN2
SCS runoff curve number for moisture condition II (/)
II
[-20, 15]
-9.770
SURLAG
Surface runoff lag coefficient (/)
I
[0.1, 12]
2.235
SOL_AWC
Available water capacity of the soil layer (mm H2O/mm soil)
II
[-0.2, 0.2]
0.1103
CH_K1
Hydraulic conductivity in tributary channel alluvium (mm/h)
I
[1, 25]
1.458
CH_N1
Manning's "n" value for the tributary channels (/)
I
[0.01, 0.2]
0.0511
CH_K2
Hydraulic conductivity in main channel alluvium (mm/h)
I
[1, 25]
7.094
ALPHA_BNK
Baseflow alpha factor for bank storage (day)
I
[0.001, 1]
0.0047
SHALLST
Initial depth of water in the shallow aquifer (mm H2O)
I
[500, 2000]
644.60
GW_DELAY
Groundwater delay time (day)
I
[10, 800]
723.69
NPERCO
Nitrogen percolation coefficient (/)
I
[0.01, 1]
0.0843
SDNCO
Denitrification threshold water content (/)
I
[0.2, 1.2]
0.2556
AC
CE
PT
ED
M
AN US
CR IP T
Parameters
52
ACCEPTED MANUSCRIPT
Table 3. Hypothetical future loading scenarios for predicting potential water quality responses Loading scenario
Synthetic case
Real-situation case
Point source A a
210 kg N/day
90 kg N/day
a
105 kg N/day
45 kg N/day
372 kg N/ha/year
356 kg N/ha/year
Point source B Fertilizer use
b
Scenario 2 (steady loading) Point source A a Point source B Fertilizer use
a
b
140 kg N/day
60 kg N/day
70 kg N/day
30 kg N/day
248 kg N/ha/year
238 kg N/ha/year
Scenario 3 (decreased loading) 70 kg N/day
30 kg N/day
a
35 kg N/day
15 kg N/day
Point source B Fertilizer use
b
124 kg N/ha/year
AN US
Point source A a
CR IP T
Scenario 1 (increased loading)
119 kg N/ha/year
Constant loading.
b
Fertilizer application is scheduled on specific dates.
AC
CE
PT
ED
M
a
53
ACCEPTED MANUSCRIPT
Table 4. Efficiency of DREAM(ZS) in both the synthetic and real-situation cases AR a
NC b
[θ, φ, w]
5.17%
489,840
w
62.83%
2640
MCMC_IM1
[θ, φ]
2.99%
350,160
MCMC_IM2
[θ, φ]
2.96%
328,320
MCMC_IM3
[θ, φ]
4.02%
288,000
MCMC_IM4
[θ, φ]
3.97%
357,600
MCMC_IM5
[θ, φ]
1.88%
325,680
BAIPU_RMs
[θ, φ, w]
23.19%
143,040
w
66.25%
2,400
MCMC_RM1
[θ, φ]
4.69%
141,840
MCMC_RM2
[θ, φ]
6.27%
143,760
MCMC_RM3
[θ, φ]
5.47%
107,760
MCMC_RM4
[θ, φ]
5.39%
173,040
MCMC_RM5
[θ, φ]
6.84%
100,080
BMA_IMs Synthetic case
BMA_RMs Real-situation case
AN US
BAIPU_IMs
CR IP T
Parameters
Approach
AR is the acceptance rate.
b
NC is the number of samples for the Markov chains to achieve convergence.
AC
CE
PT
ED
M
a
54
ACCEPTED MANUSCRIPT Table 5. Correlation coefficients a between model parameters and the weights derived by BAIPU in the synthetic case w2
w3
w4
w5
ESCO
0.08
0.06
0.07
-0.12
0.02
CN2
0.08
0.06
0.00
0.01
-0.08
SURLAG
0.02
0.04
-0.08
0.09
-0.07
SOL_AWC
0.06
0.03
0.11
0.03
-0.17
CH_K1
-0.05
-0.13
0.00
0.00
0.10
CH_N1
0.02
0.05
-0.08
0.09
-0.08
CH_K2
-0.02
0.04
-0.04
0.05
ALPHA_BNK
0.03
-0.03
0.00
0.04
SHALLST
-0.27
-0.12
-0.14
0.02
GW_DELAY
-0.12
-0.07
-0.07
0.11
NPERCO
-0.08
-0.07
-0.09
-0.04
SDNCO
0.17
0.07
0.12
0.07
-0.03 -0.05 0.27 0.01 0.18
-0.29
AN US
a
CR IP T
w1
5
Because the sample size is huge (above 10 ), all the correlations are statistically significant at the
confidence level of 99% (i.e., p<0.01), except those with 0.00 coefficients. The coefficients larger
AC
CE
PT
ED
M
than 0.1 or smaller than -0.1 are shaded in the table to indicate relatively strong correlations.
55
ACCEPTED MANUSCRIPT Table 6. Correlation coefficients a between model parameters and the weights derived by the BAIPU in the real-situation case w2
w3
w4
w5
ESCO
0.04
-0.01
-0.05
-0.03
0.03
CN2
0.03
-0.02
0.01
0.00
-0.02
SURLAG
-0.05
0.00
0.01
0.00
0.04
SOL_AWC
0.09
0.00
-0.07
0.04
-0.09
CH_K1
-0.04
0.05
-0.02
-0.01
0.02
CH_N1
-0.03
0.02
0.00
0.00
0.02
CH_K2
-0.02
0.02
0.02
0.01
ALPHA_BNK
0.00
-0.01
-0.02
0.03
SHALLST
-0.09
0.02
-0.01
0.06
GW_DELAY
-0.11
-0.01
0.01
0.09
NPERCO
-0.23
-0.04
0.07
0.03
SDNCO
0.13
0.00
-0.10
-0.03
-0.02 -0.02 -0.01 -0.02 0.14
-0.03
AN US
a
CR IP T
w1
5
Because the sample size is huge (above 10 ), all the correlations are statistically significant at the
confidence level of 99% (i.e., p<0.01), except those with 0.00 coefficients. The coefficients larger
AC
CE
PT
ED
M
than 0.1 or smaller than -0.1 are shaded in the table to indicate relatively strong correlations.
56
ACCEPTED MANUSCRIPT
Figure captions Figure 1. The ―returning‖ operation proposed for sampling model weights using DREAM(ZS). Figure 2. The framework of the Bayesian Analysis of Input and Parametric
CR IP T
Uncertainties (BAIPU).
Figure 3. The study area and modeling domain. (a) The Newport Bay Watershed (NBW) and (b) location of the NBW.
AN US
Figure 4. Cumulative distribution functions (CDFs) of model parameters inferred by the BAIPU and MCMC analyses for the individual models in the synthetic case. (a)-(l) are for the twelve different SWAT parameters.
M
Figure 5. Marginal probability distribution functions of model weights derived based
ED
on BAIPU (red bars) and the sequential BMA approach (grey bars) in the synthetic case. Red points and black circles indicate the respective mean values. (a)-(e)
PT
correspond to the five models.
CE
Figure 6. Results of the uncertainty analysis in the synthetic case. (a) and (b) illustrate
AC
the uncertainty bands (50% and 95% confidence intervals) of nitrate concentration based on BAIPU and the sequential BMA approach, respectively. (c) and (d) compare the means and standard deviations of the predictive distributions based on the two approaches. Each circle in the subplots (c) and (d) represents a simulation day. Figure 7. Quantile-Quantile (QQ) plots of nitrate concentration based on the BAIPU, sequential BMA and individual MCMC analyses. (a) and (b) are for the synthetic and 57
ACCEPTED MANUSCRIPT
real-situation cases, respectively. The solid black lines are the 45-degree lines. Figure 8. True and simulated values of the two management variables under the three future loading scenarios in the synthetic modeling case. (a), (c) and (e) are the results
exceeding (FOE) 6 mg N/L.
CR IP T
of average nitrate concentration (AVC). (b), (d) and (f) are the results of frequency of
Figure 9. Cumulative distribution functions (CDFs) of uncertain model parameters
inferred by the BAIPU and individual MCMC analyses in the real-situation modeling
AN US
case. (a)-(i) are for the twelve model parameters.
Figure 10. Marginal probability distribution functions of model weights derived based on BAIPU (red bars) and the sequential BMA approach (grey bars) in the
M
real-situation modeling case. Red points and black circles indicate the respective
ED
mean values. (a)-(e) correspond to the five models. Figure 11. Results of the uncertainty analysis in the real-situation modeling case. (a)
PT
and (b) illustrate the uncertainty bands (50% and 95% confidence intervals) of nitrate
CE
concentration based on the BAIPU and the sequential BMA approach, respectively. (c)
AC
and (d) compare the means and standard deviations of the predictive distributions based on the two approaches. Each circle in the subplots (c) and (d) represents a simulation day. Figure 12. Simulated values of the two management-relevant variables under the three future loading scenarios in the real-situation modeling case. (a), (c) and (e) are the results of average nitrate concentration (AVC). (b), (d) and (f) are the results of 58
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
frequency of exceeding (FOE) 6 mg N/L.
59
ACCEPTED MANUSCRIPT
-∆o
∆o
θ*,i
θb
θ**,i ∆
boundary
boundary
CR IP T
θi
Figure 1. The ―returning‖ operation proposed for sampling model weights using
AC
CE
PT
ED
M
AN US
DREAM(ZS).
60
ACCEPTED MANUSCRIPT
Input X2
. . .
Input XK
Model M1
Model M2
. . .
Model MK
Hyperparameters φ
MLE: EM algorithm
The residual error model
p(Y|X1,θ,φ)
CR IP T
Model parameters θ
Input X1
p(Y|X2,θ,φ)
w1(θ,φ)
w2(θ,φ)
MCMC sampling
. . .
p(Y|XK,θ,φ)
. . . wK(θ,φ)
AN US
BMA prediction: p(Y|X1,...,K,θ,φ)
No
ED
MCMC sampling
Likelihood function: l(θ,φ,w1,...K,X1,...K,Z)
M
Observations Z
Has the Markov chain converged ? Yes
Stop the sampling
PT
Predictive distribution
CE
Figure 2. The framework of the Bayesian Analysis of Input and Parametric
AC
Uncertainties (BAIPU).
61
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 3. The study area and modeling domain. (a) The Newport Bay Watershed
AC
CE
PT
ED
M
(NBW) and (b) location of the NBW.
62
ACCEPTED MANUSCRIPT
1
1
1
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
(a) ESCO
0 0
0.1
0.2
0.3
0.4
0.5
CDF
0.8
(d) SOL_AWC
0.2
0.3
0.4
0.8
0.8
0.4
0.4
0.2
0.2 0.5
0.7
0.9 1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
(e) CH_K1
0.3
0.4
0.5
1 (j) GW_DELAY
0.6
0.1
0
0.4
PT
0
0.6
0.8
0.2
0.8
11
0
1
1
0.8
0.8
CE
0
0.4
0.6
0.8
1
0.6
0.8
1
(l) SDNCO 0.2 0
0
1 0
0.2
0.4 (k) NPERCO
0.05
0.1
0.15
0.2
0
MCMC_IM1
MCMC_IM2
MCMC_IM5
True value
0.04
AC
MCMC_IM4
(i) SHALLST
0.6
0
BAIPU_IMs
0.4
0
0.0070.6
0.2
1
0.3
0.4
0.4
0.2
0.2
0.6
0.6
0.4
0.1
0.8
(h) ALPHA_BNK
ED
0.8
0.025 0.05 0.075
M
0.2
(f) CH_N1
1
0 0.003
0
0.25
0
0.2
(g) CH_K2
0.2
0.4
0
1
0.15
0.6
0 0.3
0.1 1
0.6
0.1
CDF
0.1
1
0.6
0
CDF
0 0
1
(c) SURLAG
CR IP T
0
0.2
(b) CN2
AN US
CDF
0.8
0.2
0.4
MCMC_IM3
0.08
Figure 4. Cumulative distribution functions (CDFs) of model parameters inferred by the BAIPU and MCMC analyses for the individual models in the synthetic case. (a)-(l) are for the twelve different SWAT parameters.
63
ACCEPTED MANUSCRIPT
0.2
0.1
0.0 0.0
0.5
0.5
1.0
(c) w3 0.2
0.0 0.0
0.5
1.0
0.6
(d) w4 0.3
0.2
0.4
0.6
(e) w5
Sequential BMA BAIPU
0.3
CR IP T
Posterior density
Posterior density
0.1
0.0 0.0
1.0
0.6
0.0 0.0
0.4
(b) w2
Posterior density
(a) w1
Posterior density
Posterior density
0.2
Synthetic case
0.0 0.0
0.2
0.4
AN US
Figure 5. Marginal probability distribution functions of model weights derived based on BAIPU (red bars) and the sequential BMA approach (grey bars) in the synthetic case. Red points and black circles indicate the respective mean values. (a)-(e)
AC
CE
PT
ED
M
correspond to the five models.
64
ACCEPTED MANUSCRIPT
60
95% conf. int.
-2 -12 Median
NO3-N (mg N/L)
50
50% conf. int.
(a) BAIPU
Observations
40 30 20 10
CR IP T
0
60
(b) Sequential BMA 40
30
AN US
20 10 0
16
(c) Predictive Mean
M
20
ED
15 10 5 0 0
PT
By sequential BMA (mg N/L)
25
5
10
15
By sequential BMA (mg N/L)
NO3-N (mg N/L)
50
(d) Predictive std. dev.
12
8
4
0 20
0
25
4
8
12
16
By BAIPU (mg N/L)
CE
By BAIPU (mg N/L)
AC
Figure 6. Results of the uncertainty analysis in the synthetic case. (a) and (b) illustrate the uncertainty bands (50% and 95% confidence intervals) of nitrate concentration based on BAIPU and the sequential BMA approach, respectively. (c) and (d) compare the means and standard deviations of the predictive distributions based on the two approaches. Each circle in the subplots (c) and (d) represents a simulation day.
65
ACCEPTED MANUSCRIPT
1
0.8
0.6 0.4 0.2
0
(b) NO3; Real-situation case
0.8
0.6 0.4 0.2
0 0
0.2 0.4 0.6 0.8 Theoretical Quantile of U[0,1]
1
0
0.2 0.4 0.6 0.8 Theoretical Quantile of U[0,1]
1 BAIPU_IMs(RMs)
BMA_IMs(RMs)
MCMC_IM2(RM2) 0
MCMC_IM3(RM3)
0
CR IP T
(a) NO3; Synthetic case
Quantile of observed p-values
Quantile of observed p-values
1
1
MCMC_IM1(RM1) MCMC_IM4(RM4) 0.1
AN US
MCMC_IM5(RM5)
Figure 7. Quantile-Quantile (QQ) plots of nitrate concentration based on the BAIPU, sequential BMA and individual MCMC analyses. (a) and (b) are for the synthetic and
AC
CE
PT
ED
M
real-situation cases, respectively. The solid black lines are the 45-degree lines.
66
ACCEPTED MANUSCRIPT
1
(a) Scenario 1; AVC
0.8
Cumulative frequency
0.6 0.4 0.2 0 12
14
16
18
0.6 0.4 0.2
1
0.8
Cumulative frequency
(c) Scenario 2; AVC
0.6 0.4
0 8
10
0.8
0.94
0.95
12
0.6 0.4 0.2 0 0.75
14
1
0.96
(d) Scenario 2; FOE
AN US
0.2
0.8
0.85
0.9
0.95
0.5
0.6
1
0.8
M
0.6 0.4 0.2
ED
(e) Scenario 3; AVC
0 4
5
6
7
PT
0
(f) Scenario 3; FOE
0.6 0.4 0.2
0.2
0.3
0.4
Frequency of exceeding
BMA_IMs
MCMC_IM1
MCMC_IM2
MCMC_IM4
MCMC_IM50.1
True Model
CE
MCMC_IM3
0
0.8
0 8
Average NO3-N concentration BAIPU_IMs 1
Cumulative frequency
Cumulative frequency
(b) Scenario 1: FOE
0 0.93
20
1
Cumulative frequency
0.8
CR IP T
Cumulative frequency
1
AC
Figure 8. True and simulated values of the two management variables under the three future loading scenarios in the synthetic modeling case. (a), (c) and (e) are the results of average nitrate concentration (AVC). (b), (d) and (f) are the results of frequency of exceeding (FOE) 6 mg N/L.
67
1
1
1
0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
(a) ESCO
0
CDF
0
0.1
0.2
0.3
0
1
1
0.8
0.8
0.6
0.6
0.4
0.4 (d) SOL_AWC
0.2 0.3
0.4
0.5
0.6
0.4 0.2
(e) CH_K1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.25
0.3
0
0.2
(h) ALPHA_BNK
0
0.005
0.01
0.4
0.8
0.8
0.8
0.6
0.6
0.4
0.4
ED
1
PT
0.2
0.2
0
0.6
0.8
1 0
0 1
0.8
1
(l) SDNCO
0 0
0.1
0.2
BAIPU_RMs
MCMC_RM1
MCMC_RM3 0
MCMC_RM4
AC
CE
0.4
0.6
0.2
(k) NPERCO
(j) GW_DELAY
(i) SHALLST
0
1
0.4
0.1
0.4
1
0.6
0.05
0.6
0 0.2
0.01
0.8
M
0
0.005
1
0.2
(g) CH_K2
(f) CH_N1
0
0
1
0.2
0.8
0 0.2
CDF
0.1
1
0.2
0
CDF
0.05
(c) SURLAG
0 0.15
CR IP T
0
0.2
(b) CN2
AN US
CDF
ACCEPTED MANUSCRIPT
0.3
0
0.5
1
MCMC_RM2 0.1
MCMC_RM5
Figure 9. Cumulative distribution functions (CDFs) of uncertain model parameters inferred by the BAIPU and individual MCMC analyses in the real-situation modeling case. (a)-(i) are for the twelve model parameters.
68
ACCEPTED MANUSCRIPT
(a) w1 0.2
0.4
0.8
(d) w4 0.15
0.00 0.0
0.25
0.00 0.0
0.3
0.6
(c) w3 0.3
0.0 0.00
0.25
0.50
0.30
Posterior density
Posterior density
0.30
(b) w2
0.5
(e) w5 Sequential BMA BAIPU
0.15
CR IP T
0.0 0.0
0.6
Posterior density
0.50
Posterior density
Posterior density
0.4
Real-situation case
0.00 1.0 0.0
0.5
1.0
AN US
Figure 10. Marginal probability distribution functions of model weights derived based on BAIPU (red bars) and the sequential BMA approach (grey bars) in the
real-situation modeling case. Red points and black circles indicate the respective
AC
CE
PT
ED
M
mean values. (a)-(e) correspond to the five models.
69
ACCEPTED MANUSCRIPT
30
95% conf. int.
-2 -12 Median
NO3-N (mg N/L)
25
50% conf. int.
(a) BAIPU
Observations
20 15 10 5
CR IP T
0
30
(b) Sequential BMA 20
15
AN US
10 5 0
7
(c) Predictive Mean
M
20
ED
15 10 5 0 0
PT
By sequential BMA (mg N/L)
25
5
10
15
By sequential BMA (mg N/L)
NO3-N (mg N/L)
25
(d) Predictive std. dev.
6 5 4 3 2
20
2
25
3
4
5
6
7
By BAIPU (mg N/L)
CE
By BAIPU (mg N/L)
AC
Figure 11. Results of the uncertainty analysis in the real-situation modeling case. (a) and (b) illustrate the uncertainty bands (50% and 95% confidence intervals) of nitrate concentration based on the BAIPU and the sequential BMA approach, respectively. (c) and (d) compare the means and standard deviations of the predictive distributions based on the two approaches. Each circle in the subplots (c) and (d) represents a simulation day. 70
ACCEPTED MANUSCRIPT
1
(a) Scenario 1; AVC
0.8
Cumulative frequency
0.6 0.4 0.2 0 10
11
12
13
0.6 0.4 0.2
Cumulative frequency
1
0.8 0.6 0.4
0
0.9
0.8 0.6 0.4
0.95
(d) Scenario 2; FOE
0.2
AN US
(c) Scenario 2; AVC
0.2
0
7
7.5
8
8.5
1
9
0.6
0.65
0.7
0.75
0.8
1
0.8
M
0.6 0.4
(e) Scenario 3; AVC
0 4
ED
0.2
4.5
Cumulative frequency
Cumulative frequency
0.8
0 0.85
14
1
Cumulative frequency
(b) Scenario 1; FOE
CR IP T
Cumulative frequency
1
(f) Scenario 3; FOE
0.8 0.6 0.4 0.2
0 5
0
PT
Average NO3-N concentration BAIPU_RMs 1
0
0.2
0.3
BMA_RMs
MCMC_RM1
MCMC_RM4
MCMC_RM5
MCMC_RM2
0.08
CE
MCMC_RM3 0
0.1
Frequency of exceeding
AC
Figure 12. Simulated values of the two management-relevant variables under the three future loading scenarios in the real-situation modeling case. (a), (c) and (e) are the results of average nitrate concentration (AVC). (b), (d) and (f) are the results of frequency of exceeding (FOE) 6 mg N/L.
71