Uncertainty analysis of model inputs in riverine water temperature simulations
In this study, the HFLUX model was coupled with the SCEM-UA algorithm for analyzing the uncertainties of the model inputs. The specific procedures started with selecting the inputs of the HFLUX model. With the linked HFLUX and SCEM-UA model and implementation of an iteration scheme, the uncertainty of each of the selected inputs was obtained based on the ranges (minimum and maximum values) of the input data/parameters and the Latin hypercube sampling. The simulations were then compared against the observed data to evaluate the performance of the SCEM-UA algorithm. These steps are depicted in Fig. 1.Figure 1Flowchart for the uncertainty analysis.Full size imageRiver water temperatures simulated by the HFLUX modelRiver water temperature affects the water quality and the ecosystem health, and hence control of river water temperature is important to mitigation of its adverse effects1. The HFLUX model was used to simulate the streamflow temperatures at different locations and times. The model is highly flexible in terms of choosing the solution methods for solving the governing equations and selecting the energy budget terms such as shortwave solar radiation, latent heat flux, and sensible heat transfer flux. The model input data include the initial spatial and temporal temperature conditions, stream geometry data, discharge data, and meteorological data8. The water balance and energy balance equations are respectively given by8:$$frac{partial A}{{partial t}} + frac{partial Q}{{partial x}} = mathop qnolimits_{L}$$
(1)
$$frac{{partial left( {Amathop Tnolimits_{w} } right)}}{partial t} + frac{{partial left( {Qmathop Tnolimits_{w} } right)}}{partial x} = mathop qnolimits_{L} mathop Tnolimits_{L} + R$$
(2)
$$R = frac{{Bmathop varphi nolimits_{total} }}{{mathop rho nolimits_{w} mathop Cnolimits_{w} }}$$
(3)
where A is the cross section area of the stream (m2), x is the distance along the stream (m), t is the time (s), Q is the discharge of the stream (m3/s), qL is the lateral inflow per unit stream length (m2/s), Tw is the stream temperature ((^circ C)), TL is the temperature of the lateral inflow ((^circ C)), R is the energy flux (source or sink) per unit stream length ((^circ C) m2/s), B is the width of the stream (m), (mathop varphi nolimits_{total}) is the total energy flux to the stream per surface area (W/m2), (mathop rho nolimits_{w}) is the density of water (kg/m3), and (mathop Cnolimits_{w}) is the specific heat of water (J/kg (^circ C)). Equation (3) is based on a thermal datum of 0 (^circ C) and the impact on the absolute value of the advective heat flux term. In Eq. (2), if qL is negative, the first term on the right-hand side of the equation becomes a loss of qLTw. Also, dispersive heat transport that is omitted in Eq. 2 is negligible when the longitudinal change in water temperature is small in comparison to the temporal changes8.SCEM-UA algorithmThe SCEM-UA algorithm provides posterior distribution functions for the model parameters and input data by generating an initial sample from the parameter space. First, the indicators of n, q, and s that are respectively dimension (the number of investigate inputs), number of complexes (the population to be divided), and population (the number of sample points) are determined for the algorithm. Then, the algorithm searches the sampling points in the feasible space and sorts the points according to the density. The algorithm determines the sequence and complexes based on those points. The sequence is the first q points of the population and complexes are a collection of m points from the population. Note that m = s/q. In the next step, the points of each complex are sorted based on the density, which can be mathematically expressed as20:$$left{ {begin{array}{*{20}c} {mathop alpha nolimits^{k} le T,,,,,,,,,mathop theta nolimits^{t + 1} = Nleft( {mathop theta nolimits^{t} ,,mathop Cnolimits_{n}^{2} mathop Sigma nolimits^{k} } right)} \ {mathop alpha nolimits^{k} > T,,,,,,,,mathop theta nolimits^{t + 1} = Nleft( {mathop mu nolimits^{k} ,,mathop Cnolimits_{n}^{2} mathop Sigma nolimits^{k} } right)} \ end{array} } right.$$
(4)
where k = 1,2,…,q, α is the ratio of the mean posterior density of the m points of complexes to the mean posterior density of the last m generated points of sequences, (theta) is the points of complexes, ({c}_{n}=frac{2.4}{sqrt{n}}) , (T={10}^{6}), (mu) is the mean, and ∑ denotes the covariance. To investigate the new points created by the algorithm, the points of complexes are replaced by20:$$left{ {begin{array}{*{20}l} {Omega ge Zquad replace,best,member,of,mathop Cnolimits^{k} ,with,mathop theta nolimits^{t + 1} } \ {Omega < Zquad mathop theta nolimits^{t + 1} = mathop theta nolimits^{t} ,,,,,,,,,,,,,,,,,,,,,} \ end{array} } right.$$
(5)
where (mathop Cnolimits^{k}) is the Kth complex, Z is drawn from the uniform distribution in the range of 0–1, and Ω is calculated by20:$$Omega = frac{{Pleft( {left. {mathop theta nolimits^{t + 1} } right|y} right)}}{{Pleft( {left. {mathop theta nolimits^{t} } right|y} right)}}$$
(6)
where (Pleft( {left. {mathop theta nolimits^{t + 1} } right|y} right)) and (Pleft( {left. {mathop theta nolimits^{t} } right|y} right)) are the posterior probability distributions for (mathop theta nolimits^{t + 1}) and (mathop theta nolimits^{t}), respectively. Then, the algorithm examines the following condition for each complex. If it is rejected, the algorithm replaces the worst member ({c}^{k})(the point with the lowest density) with ({theta }^{t+1}) 20.$$mathop Gamma nolimits^{k} le T,,and,,Pleft( {{{mathop theta nolimits^{t + 1} } mathord{left/ {vphantom {{mathop theta nolimits^{t + 1} } y}} right. kern-nulldelimiterspace} y}} right) < ,Pleft( {{{mathop Cnolimits_{m}^{k} } mathord{left/ {vphantom {{mathop Cnolimits_{m}^{k} } y}} right. kern-nulldelimiterspace} y}} right)$$
(7)
where ({Gamma }^{k}) is the ratio of the posterior density of the best (the point with the highest density) to the posterior density of the worst member of ({c}^{k}). The last step is to examine (beta) and L. Note that (beta) = 1 and L = m/10. If (beta < L), (beta = beta + 1) and the algorithm returns to sort complex points. Otherwise, the algorithm examines the Gelman and Rubin convergence6, and eventually provides the posterior distribution functions20. The value of the Gelman and Rubin convergence should be less than 1.2. The Gelman and Rubin convergence is examined by:$$R = sqrt {frac{g - 1}{g} + frac{q + 1}{{q.g}}frac{B}{W}}$$
(8)
where g is the number of iterations within each sequence, B is the variance between the q sequence means, and W is the average of the q within-sequence variances for the parameter under consideration20.Study AREAMeadowbrook Creek was selected to test the methods proposed in this study8. The creek flows through the City of Syracuse in New York. Thus, this catchment consists of high residential and industrial land covers, which contribute runoff to the main channel. The creek is about 4 km long. A portion of this creek (475 m long) was selected for the modeling for a period of June 13–19, 2012 in this study. The upstream boundary condition in the HFLUX model was set based on the water temperature of the creek observed at the upstream station8. The uncertainty of the model inputs was examined at three selected points as shown in Fig. 2. Note that the input values at these three points had greater relative changes than the changes at other locations, which provided the possibility to improve the evaluation of the algorithm performance. In addition, these three locations had the same sampling of the selected input data. During the simulation period, the streamflow velocity varied within a range of 0.06–0.63 (m/s). The daily temperature changed between 8.9 and 28.2 °C. The relative humidity, used to calculate the total energy flux to the stream per surface area, changed from 36 to 93%. The creek bed mainly consisted of clay, cobbles, sand, and gravel materials. The basic statistics of the data/variables used in the HFLUX model are presented in Table 1. Figure 2 shows the study area, the creek, and the three selected points for analysis.Figure 2Study area and the locations of three evaluation sections (the gray enlarged map shows the State of New York), the map in this Figure is created by Google Earth 7.0.2.8415 (https://google.com/earth/versions).Full size imageTable 1 Basic statistics of the data/variables used in the HFLUX model.Full size tableEthical approvalAll authors accept all ethical approvals.Consent to participateAll authors consent to participate.Consent to publishAll authors consent to publish. More
