- Email: [email protected]
Sediment concentration prediction and statistical evaluation for annual load estimation
Journal of Hydrology 482 (2013) 69–78
Contents lists available at SciVerse ScienceDirect
Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol
Sediment concentration prediction and statistical evaluation for annual load estimation You-Gan Wang a,⇑, Ting Tian b a b
Centre for Applications in Natural Resource Mathematics (CARM), School of Mathematics and Physics, The University of Queensland, St. Lucia 4072, Australia School of Agriculture and Food Science, The University of Queensland, St. Lucia 4072, Australia
a r t i c l e
i n f o
Article history: Received 22 September 2012 Received in revised form 19 December 2012 Accepted 19 December 2012 Available online 8 January 2013 This manuscript was handled by Geoff Syme, Editor-in-Chief Keywords: Bootstrap Load estimation Nash–Sutcliffe model efficiency coefficient Rating curve Suspended sediment Uncertainty
s u m m a r y We consider the development of statistical models for prediction of constituent concentration of riverine pollutants, which is a key step in load estimation from frequent flow rate data and less frequently collected concentration data. We consider how to capture the impacts of past flow patterns via the average discounted flow (ADF) which discounts the past flux based on the time lapsed – more recent fluxes are given more weight. However, the effectiveness of ADF depends critically on the choice of the discount factor which reflects the unknown environmental cumulating process of the concentration compounds. We propose to choose the discount factor by maximizing the adjusted R2 values or the Nash–Sutcliffe model efficiency coefficient. The R2 values are also adjusted to take account of the number of parameters in the model fit. The resulting optimal discount factor can be interpreted as a measure of constituent exhaustion rate during flood events. To evaluate the performance of the proposed regression estimators, we examine two different sampling scenarios by resampling fortnightly and opportunistically from two real daily datasets, which come from two United States Geological Survey (USGS) gaging stations located in Des Plaines River and Illinois River basin. The generalized rating-curve approach produces biased estimates of the total sediment loads by 30% to 83%, whereas the new approaches produce relatively much lower biases, ranging from 24% to 35%. This substantial improvement in the estimates of the total load is due to the fact that predictability of concentration is greatly improved by the additional predictors. Crown Copyright Ó 2013 Published by Elsevier B.V. All rights reserved.
1. Introduction The estimation of fluvial constituent loads is an essential for helping researchers and water managers understand water quality and its effect on environmental resources. Much work on estimation of constituent loads focuses on the ‘‘rating curve’’ approach on the assumption that there is a linear relationship between concentration of sediment and flow on log transformed space (Brown et al., 1970; Colby, 1956; Miller, 1951; Thomas, 1985, 1988; Verhoff et al., 1980; Walling, 1977). Generally, the linear regression analysis uses the log-transformed flow rate as the predictor. However, constituent concentrations often depend on many other environmental changes such as ‘‘first flush’’ or ‘‘depletion’’ effects. The first flush is ‘‘the first significant channelized flow of the wet season which is generally accompanied by relatively high sediment and nutrient concentrations (Wang et al., 2011). Such dynamics is also well discussed in Davis et al. (2012). There are higher flows in wet seasons, and the start of a wet season may change from year to year which makes stratification difficult (Wallace et al., 2008). High nutrient concentrations
⇑ Corresponding author. Tel.: +61 7 3365 2311. E-mail address: [email protected] (Y.-G. Wang).
often come with high flux initially, but not sustainable due to the exhaustion and slow accumulation process of sediments. Simple linear or periodic mathematical functions therefore often yield poor prediction. Therefore, seeking additional hydrologically meaningful predictors or directly modeling the physical process is fundamental in improving load estimation (Wang et al., 2011). Limitations of the rating curve method were well studied by Ferguson (1986). Since then, new statistical methods have been developed expressing the relationship between constituent concentration and flow. One is the seven-parameter log linear model for concentration, which generalizes the rating curve by incorporating five additional predictors: a quadratic term of the logarithm of flow rate, and four parametric functions for seasonality (sine and cosine functions for annual and semi-annual cycles) (Cohn et al., 1992). This pioneering work has provided fundamental guidelines for constructing an appropriate regression model for load estimation based on frequently collected flow data but less-frequently collected concentration data. Additional predictive covariates to the original 7-parameter model can better account for hysteresis in sediment transport (Wang and Linker, 2008). Integration of flow history was proposed by Singer and Dunne (2001) using ARIMA models. Such time series models can extract information from temporal correlations. This is useful if the exhaustion is very fast leading to short lag effects (a few days). Note that one parameter
0022-1694/$ - see front matter Crown Copyright Ó 2013 Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jhydrol.2012.12.043
70
Y.-G. Wang, T. Tian / Journal of Hydrology 482 (2013) 69–78
is needed for each lag and overfit can become a concern if too many parameters are used. Stratification by season or region can also improve predictions. However, stratification can be quite subjective, and specific knowledge is required to determine the strata (Thomas and Lewis, 1993). One can also consider incorporating strata heterogeneity into a regression framework by creating indicator/dummy variables. As pointed out by Wang et al. (2011), simple linear or periodic time functions will not adequately capture information such as a ‘‘flush’’, and more novel covariates characterizing these phenomena will need to be considered. Wang et al. (2011) proposed a new predictive variable, namely, the average discounted flow (ADF). This aims to provide appropriate weights for all the past flow records in a discounted way because more recent fluxes may exhaust more sediment, whereas flux occurring a long time ago is given less weight. There are increasing concerns over quality of water, leading to a demand for accurate estimates of total sediment loads. There are
(a)
1200
Concentration (mg/L)
Concentratation (mg/L)
300
several methods available to estimate total loads from frequent flow rate data and concentration data; however, frequency of concentration data is often very low, and hence concentration prediction becomes a key step in load estimation (Darnell et al., 2012). Estimation of total sediment loads therefore becomes challenging. Existing rating curve and 7-parameter methods may produce biased load estimates, and multiple regression approach appears to be generally more versatile and easy to incorporate additional covariates. Both Koch and Smillie (1986) and Ferguson (1986) found that the regression model (7-parameter) under-predicts concentrations. However, such bias due to nonlinear transformation is often small relative to biases from other sources. Recently, Verma et al. (2012) also investigated bias and standard errors for three estimators and concluded importance of accurate estimators for cost-effective monitoring designs. In this paper we focus on all models for predicting concentration and evaluate their performance in terms of R2 and the resulting load estimates. In particular, we define the new predictive variable (i.e. ADF), which can mimic
250 200 150 100 50
(b)
1000 800 600 400 200 0
0 03 Apr Jul Oct 04 Apr Jul Oct 05 Apr Jul
03 Apr Jul Oct 04 Apr Jul Oct 05 Apr Jul
Date
Date
(c)
80
Flow rate (CM/s)
Flow rate (CM/s)
4 3 2
(d)
60 40
1
20
0
0 03 Apr Jul Oct 04 Apr Jul Oct 05 Apr Jul
03 Apr Jul Oct 04 Apr Jul Oct 05 Apr Jul
Date
Date
(e)
140
Load (1000 tonnes/day)
Load (tonnes/day)
2500 2000 1500 1000 500
(f)
120 100 80 60 40 20 0
0 03 Apr Jul Oct 04 Apr Jul Oct 05 Apr Jul
Date
03 Apr Jul Oct 04 Apr Jul Oct 05 Apr Jul
Date
Fig. 1. The trends of concentration, flow rate, and load for two stations between 2003 and 2005. (a), (c), (e) Station 05532500 (Des Plaines River), (b), (d), (f) Station 05586100 (Illinois River).
71
Y.-G. Wang, T. Tian / Journal of Hydrology 482 (2013) 69–78
models in predicting the annual cumulative load, each regression model using two different sampling designs (uniform sampling and opportunistic sampling) will be described in Section 3. We subsample from the 3 years’ daily records for sediment, using uniform (once every fortnight) and opportunistic sampling (more frequent in wet periods than in dry periods) and use the adjusted R2 as a criterion to determine the most appropriate concentration predictive model. In Section 4, we present our findings; then, a discussion of the methods and possible generalization will be presented.
the recovery and exhaustion process commonly seen in temporal measures of concentration and flow; also, an optimal discount factor is sought from the data. Unlike the time series models, the number of additional parameters does not change even if exhaustion process is slow (but require discount factor close to 1). The paper is organized as follows: in Section 2 we will briefly introduce the two datasets from USGS (see http://il.water.usgs.gov/ for details about the profiler). The new models we propose are presented and discussed. To evaluate the performance of regression
2. Materials and methods Table 1 The mean values of flow (CM/s), concentration (C, mg/L) and load (tonnes/day) for the period of January 2003 to September 2005. Station 05532500 Des Plaines River
2.1. Materials
Station 05586100 Illinois River
Period
Flow
C
Load
Flow
C
Load
2003.01–2003.04 2003.05–2003.08 2003.09–2004.04 2004.05–2004.08 2004.09–2005.04 2005.05–2005.08 2005.09–2005.09
10 21 17 33 17 7 5
37 47 48 56 57 34 37
51 122 104 234 136 35 23
273 634 479 737 898 230 178
128 256 159 205 165 48 33
3963 20046 9626 16523 15310 1242 557
The data considered in this study were collected daily from 2003 to 2005 by US Geological Survey. Two stations will be considered: Station (05532500) of the Des Plaines River at Riverside, located at latitude (41°4900 200 ) and longitude (87°4900 150 ), and Station (05586100) of the Illinois River at Valley City, located at latitude (39°4200 120 ) and longitude (90°3800 430 ). Drainage area of the former Station is 1632 square kilometers (km2) and that of the latter is 69,264 square kilometers (km2). Fig. 1 shows the three time series plots (daily measurements of flow rate, concentration and load from January 1, 2003 to Septem-
Station 05532500 5.5 5.0
(a)
(b)
True Concentration
4.5 4.0 3.5 3.0 2.5 5.5 5.0
3.0
3.5
4.0
4.5
5.0
5.5
6.0
(c)
2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
3.0
3.5
4.0
4.5
5.0
5.5
6.0
(d)
4.5 4.0 3.5 3.0 2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
2.5
Predicted Concentration
True Concentration
Station 05586100 6.5 6.0 5.5 5.0 4.5 4.0 3.5 3.0
(b)
(a)
3 6.5 6.0 5.5 5.0 4.5 4.0 3.5 3.0
4
5
6
3
4
5
6
3
4
5
6
(d)
(c)
3
4
5
6
Predicted Concentration Fig. 2. Plot of predicted concentrations based on four regression models (uniform sampling) versus the true concentration for two stations during 2005 (Log Scale Space). (a) Rating curve: black solid circle ( ), (b) seven-parameter model: red diamond ( ), (c) M1: green triangle point-up ( ), (d) M2: blue square ( ). Solid line: predicted concentration equal to true concentration.
72
Y.-G. Wang, T. Tian / Journal of Hydrology 482 (2013) 69–78
ber 30, 2005) for the two stations (05532500 and 05586100). The concentration data in Station 05532500 shows a periodic pattern, whereas the concentration in Station 05586100 shows two global peaks in the middles of year 2003 and year 2004. It seems concentration may be seasonal, which can be modeled by periodic function. The flow rate in each station reaches annual peaks in the middle of 2003 and 2004, but reaches the annual peak at the beginning of 2005. Similarly, daily suspended sediment loads reach a local peak in the middles of 2003 and 2004, but reaches the annual peak at the beginning of 2005 in each station. 2.2. Methods Suppose the flow rate Q(t) (e.g., in L/s) and concentration C(t) (e.g., in mg/L) are given, then the total load L can be written as the following equation:
Z
L¼K
T
CðtÞQðtÞdt;
ð1Þ
0
where K is a unit-conversion constant depending on the units applied, T is time measured in years. In practice, the samples are only collected at discrete time intervals (daily, weekly or even longer intervals). As a result, the total load L in Eq. (1) may be approximated by
LK
M X
C m Q m d;
ð2Þ
models only make use of the current flow rate and ignore the past rates. To this end, Wang et al. (2011) introduced ADF. This new variable is constructed using weighted (discounted) mean of the previous flux. The discount is based on the time lapsed since the flow. The discount factor can be so chosen to optimize the model fit (such as coefficient of determination or R2 value). Note that the R2 is also known as the Nash–Sutcliffe model efficiency in hydrology. This modeling approach greatly improved the concentration prediction model for the Tully River in the Great Barrier Reef catchments, Australia. The idea is to employ exponential discounting to all the flow historic data. Suppose (Q1, Q2, . . ., Qj) are the flow rates from day 1 to day j. The ADF with discount factor d at day j is defined by the following equation:
ADFðdÞ ¼
b L¼K
M X
^ m Q d; C m
j1
j1
1 þ d þ þ d
Q1
ð7Þ
;
where discount factor d can be measured from 0.1 to 1. This Eq. (7) summarizes past flow rates (Q1, Q2, . . ., Qj) as an index by giving more weight to more recent events. Roughly, d can be interpreted as the exhaustion rate. We propose to find the optimal d in terms of adjusted R2 values by taking account of the number of parameters used. The ADF will be used to replace the current flow rate Q(t). We first introduce our proposed new regression model, M1, within the extended seven-parameter model:
lnðCÞ ¼ b0 þ b1 lnðADFÞ þ b2 lnðQ Þ þ b3 ½lnðQ Þ2 þ b4 sinð2pTÞ þ b5 cosð2pTÞ þ b6 sinð4pTÞ þ b7 cosð4pTÞ þ e
m¼1
where d is the time interval and dM = T, Qm and Cm are the true values of flow and concentration, respectively, at time tm (the beginning of the mth time interval). However, the concentration data is not always available, so the concentration has to be predicted (at small intervals, e.g. every 10 min.) from a generalized rating curve (an empirical relationship between C and Q); the total load L can be estimated by
Q j þ dQ j1 þ þ d
ð8Þ
The optimal d is the one that produces the highest adjusted R2 values. It is possible that an improved model would employ two ADF variables with two different discount rates, d1 and d2:
lnðCÞ ¼ b0 þ b1 lnðADF 1 Þ þ b2 lnðADF 2 Þ þ b3 ½lnðADF 1 Þ2 þ b4 ½lnðADF 2 Þ2 þ b5 sinð2pTÞ þ b6 cosð2pTÞ þ b7 sinð4pTÞ þ b8 cosð4pTÞ þ e
ð9Þ
ð3Þ
m¼1
b L¼K
2 b i Q exp s ; C i 2 i¼1
M X
10
rating curve seven−parameter model M1 M2
8
Density
where the unknown concentration of Cm in Eq. (2) is substituted b m , which is often predicted from a statistical using estimates of C model such as rating curve or other regression model. Generally, b m is predicted on log scale, and the bias corrected total load estiC mator is
6 4
ð4Þ
2
where s2 is the estimated variance of the model error on log scale (see Ferguson, 1986; Koch and Smillie, 1986; Cohn, 1995).
0 0.2
0.4
0.6
0.8
2
1.0 2
Adjusted R for Station05532500 (1632km ) 2.2.1. Regression models We now briefly introduce two widely used regression models. The rating curve approach is based on the following regression model:
ð5Þ
where ln() denotes the natural logarithm function, b0 and b1 are parameters to be estimated, and e is random error, which is identical independent distributed (i.i.d). A total of k = 7 parameters are supposed to be required to formulate the relationship between concentration and flow rate:
6 4 2 0 0.2
0.4
0.6
0.8
1.0
Adjusted R2 for Station05586100 (69264km2)
2
lnðCÞ ¼ b0 þ b1 lnðQ Þ þ b2 ½lnðQ Þ þ b3 sinð2pTÞ þ b4 cosð2pTÞ þ b5 sinð4pTÞ þ b6 cosð4pTÞ þ e;
8
Density
lnðCÞ ¼ b0 þ b1 ln Q þ e;
rating curve seven−parameter model M1 M2
10
ð6Þ
where the errors, e Nð0; r2e Þ, and b are parameters of the model that can be estimated from the data. The above two regression
Fig. 3. Comparisons of adjusted R2 from four different regression models (uniform sampling) based on the relationship between concentration and flow rate. Rating curve: black dashed line ( ), seven-parameter model: red solid line ( ), M1: green dot dash line ( ), M2: blue dotted line ( ).
73
Y.-G. Wang, T. Tian / Journal of Hydrology 482 (2013) 69–78
If we restrict d1 close to 0, ADF1 will become Q. Using multiple ADF variables will capture exhaustion processes with different cycles. Again, we will choose the optimal d values to maximize the adjusted R2. We will refer to this model as M2. The adjusted R2 values can be obtained for a range of discount values: 0, 0.1, 0.2, . . ., 0.9, 0.95, 0.96, . . ., 0.99, and hence the optimal d can be chosen. 2.2.2. Resampling In order to compare the four load estimation methods based on above regression models, we apply two sampling designs, i.e. uniform sampling and opportunistic sampling. Uniform sampling randomly selects one sample from every fixed period of two weeks (sampling interval is not fixed but has a mean of 2 weeks), while opportunistic sampling randomly selects one sample from a shorter (1 week) or longer (1 month) period depending on wet or dry seasons. The full data set incorporates 1004 daily flow rate observations and corresponding concentration data, covering from 1st of January for 2003 to 30th of September for 2005. To achieve a moderate annual sample size (20–30), we will examine fortnightly and monthly sampling schemes. Table 1 shows the mean values for the selected wet and dry seasons. For illustrative purpose, we will consider random uniform sampling and opportunistic sampling using a subjectively defined ‘‘wet/dry’’ seasons. The uniform random samples are selected every fortnight (resulting in 66 samples during these 3 years). For opportunistic sampling, samples are se-
lected once a month in the dry season (from September to next May) and once a week in the wet season (between June and August). As a result, the total sample size using opportunistic sampling is the same as that using the uniform sampling. These two sampling designs have a total sample size of 66 in each simulation. Note that the resampling process is repeated for many times and overall the resampled data will preserve both the variability of original flow data (albeit at different frequency) and their relationship with corresponding concentration. The two stations are located on different rivers, where dry and wet seasons varies over years as shown in Fig. 1. This will allow us to demonstrate how the two sampling designs might differ. In terms of sampling frequency, weekly samples can also be examined. The estimators discussed here are all regression-based, unlike other imputation or composite methods, relative performance, will not so depend on the sample frequency. 3. Results 3.1. Uniform sampling 3.1.1. Concentration prediction based on four regression models The performance of predicted daily concentration where the unit for daily concentration is mg/L in the log-space in the two stations is depicted in Fig. 2. For each predictive regression model, the
Density
Station 05532500 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0
(a)
1.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0
2.0
2.5
3.0
(c)
1.5
(b)
μ=2.10 σ=0.13
3.5
4.0
1.5
2.5
3.0
3.5
2.0
2.5
3.0
(d)
μ=2.27 σ=0.15
2.0
μ=2.25 σ=0.16
4.0
1.5
3.5
4.0
μ=2.45 σ=0.20
2.0
2.5
3.0
3.5
4.0
Estimated Load (104 tonnes)
Density
Station 05586100 1.2 1.0 0.8 0.6 0.4 0.2 0.0
(a)
1 1.0
2
3
4
(c)
0.8
(b)
μ=4.00 σ=0.39
5
6
1
μ=3.22 σ=0.46
2
3
4
(d)
μ=3.14 σ=0.52
5
6
μ=2.65 σ=0.44
0.6 0.4 0.2 0.0 1
2
3
4
5
6
1
2
3
4
5
6
Estimated Load (106 tonnes) Fig. 4. Density plots of the four estimators (uniform sampling) based on rating curve (a), seven-parameter model (b), M1 (c), and M2 (d) of total loads for the two stations in 2005. Dash vertical line: True total sediment loads. Top panel: Station 05532500 (Des Plaines River), Bottom panel: Station 05586100 (Illinois River).
74
Y.-G. Wang, T. Tian / Journal of Hydrology 482 (2013) 69–78
performance of predictive logarithm concentration seems to be close to the true logarithm concentration in Station 05532500. Typically, when the values of logarithm concentration are below 4, predictive logarithm concentrations significantly match their corresponding logarithm true concentrations. However, the values of predictive logarithm concentration between 3.5 and 4.5 underestimate true logarithm concentration, while they overestimate true logarithm concentration over 5.0 for four regression models in Station 05586100. The reasons or causes for such patterns are unknown implying some other important factors affecting the concentration are to be found. In particular, the rating curve provides larger predictive logarithm concentration at large values of true logarithm concentration. Nevertheless, the values of predictive logarithm concentration for M1 and M2 are relatively closer to the true logarithm concentration compared to with the rating curve and seven-parameter model. Overall, four regression models provide more accurate predictive logarithm concentration at relatively lower true logarithm concentration (below 4.5), and M1 and M2 yield comparatively better results. 3.1.2. Adjusted R2 of prediction concentration based on four regression models Fig. 3 illustrates the effect of average discount flow on the prediction concentration based on log linear flow rate. We randomly
resample the flow rate over the timeframe and the corresponding concentration of two stations 1000 times, and compare the density of each adjusted R2 of four regression models. It is apparent that the adjusted R2 of new regression model one (M1) and new regression model two (M2) can be improved in each station. Choosing the optimal discount factor, d, in M1 (Eq. (8)) and M2 (Eq. (9)) for 66 random subsets of data samples associated with the highest adjusted R2, it is clear that M1 and M2 based on the optimal discount factor, d by adjusted R2, provide higher adjusted R2 than that generalized from the rating curve and seven-parameter model. The adjusted R2 of Station 05532500 ranges for all regression models as follow: rating curve approach (0.18, 0.7), seven-parameter model (0.2, 0.76), M1 (0.25,0.82) and M2 (0.3,0.84), in addition, the peaks of the adjusted R2 for M1 and M2 are higher than those for the rating curve approach (0.43) and seven-parameter model (0.52), increasing to 0.58 and 0.62, respectively. Meanwhile, the adjusted R2 of Station 05586100 ranges from 0.4 to 0.8 for M1 and M2, respectively, which is larger than that obtained from the rating curve approach (0.2, 0.6) and seven-parameter model (0.22, 0.74). Also, the peaks of the adjusted R2 for M1 and M2 are higher than those for the rating curve approach (0.4) and seven-parameter model (0.48), increasing to 0.61 and 0.65, respectively Thus, above results of these two stations indicate that M1 and M2 yield higher adjusted R2.
Station 05532500 5.5 5.0
(a)
(b)
True Concentration
4.5 4.0 3.5 3.0 2.5 5.5 5.0
3.0
3.5
4.0
4.5
5.0
5.5
6.0
(c)
2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
3.0
3.5
4.0
4.5
5.0
5.5
6.0
(d)
4.5 4.0 3.5 3.0 2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
2.5
Predicted Concentration
True Concentration
Station 05586100 6.5 6.0 5.5 5.0 4.5 4.0 3.5 3.0
(b)
(a)
3 6.5 6.0 5.5 5.0 4.5 4.0 3.5 3.0
4
5
6
3
4
5
6
3
4
5
6
(d)
(c)
3
4
5
6
Predicted Concentration Fig. 5. Plotof predicted concentration based on four regression models (opportunistic sampling) versus the true concentration for the two stations during 2005 (Log Space). (a) Rating curve: black solid circle ( ), (b) seven-parameter model: red diamond ( ), (c) M1: green triangle point-up ( ), M2: (d) blue square ( ).Solid line: predicted concentration equal to true concentration.
75
Y.-G. Wang, T. Tian / Journal of Hydrology 482 (2013) 69–78
3.2. Opportunistic sampling 3.2.1. Concentration prediction based on the four regression models As shown in Fig. 5, it is clear that, by means of opportunistic sampling, there is a similar performance of predictive logarithm concentration in Station 05532500 compared to the uniform sampling approach (upper panel of Fig. 2). However, the performance of predictive logarithm concentration at Station 05586100 is very different from that based on uniform sampling (lower panel of Fig. 2). There are more variations when the values of logarithm concentration exceed 5.0 mg/L. Typically, rating curve and sevenparameter models provide larger values of predictive logarithm concentration than M1 and M2 do at high concentrations. Moreover, the ranges of predicted concentration (on logarithm) by the rating curve and the seven-parameter model are wider, in which the upper bounds of predicted concentration are as high as 5.5 and 6.0, respectively. The values of predictive logarithm concentration overestimate true values for all regression models at their lower bounds. 3.2.2. Adjusted R2 of prediction concentration based on four regression models Compared with Fig. 3, the lower bounds of adjusted R2 for M1 and M2 increase to 0.3 and 0.4, respectively, while the upper bound of adjusted R2 for M1 and M2 increases to 0.85 in Station 05532500 (upper panel of Fig. 6). Corresponding peak values of adjusted R2 for M1 and M2 increase to 0.63 and 0.7, respectively. Moreover, the lower bounds of adjusted R2 of the four regression models slightly increase by 0.2 (the minimum values of adjusted R2 of rating curve and seven-parameter model from 0.2 to 0.4 and the minimum values of adjusted R2 of M1 and M2 from 0.4 to 0.6) in Station 05586100 (lower panel of Fig. 6). The peaks of adjusted R2 values of M1 and M2 increase to around 0.72. Overall, the values of adjusted R2 of the four regression models increase for opportunistic sampling, and there is an increase of adjusted R2 in M1 and M2.
rating curve seven−parameter model M1 M2
10
Density
8 6 4 2 0 0.2
0.4
0.6
0.8
2
1.0 2
Adjusted R for Station05532500 (1632km ) 10
rating curve seven−parameter model M1 M2
8
Density
3.1.3. Estimation of loads Estimates of total sediment loads for two stations in 2005 include the estimator based on rating curve ðb L rc Þ, the estimator based on seven-parameter model ðb L 7p Þ, the estimator based on M1 ðb L M1 Þ, and the estimator based on M2 ðb L M2 Þ. The estimates of total sediment load (Eq. (4)) based on four regression models of the two stations are shown in Fig. 4. We randomly resample 1000 times 66 subset samples to estimate total sediment loads. Although all four estimators of total sediment load underestimate the true load for Station 05532500 (upper panel of Fig. 4.), the averages of b L M1 and b L M1 are closer to true total sediment loads (the biases are 20% and 15%, while the bias in b L rc and b L 7p are 27% and 21%, respectively). On the other hand, b L rc and b L 7 produce overestimates, whereas b L M1 and b L M2 produce much less biased estimates especially for the second Station (Fig. 4). The rating curve produces larger estimates than the seven-parameter model, M1, and M2. Comparing four estimators in two stations, M1 and M2, regardless of whether the year is relatively dry or wet, provide better results on estimation of total sediment loads. In reality, we often intend to sample more often when flux is large. This opportunistic sampling makes sense for more reliable load estimation because the load is a product of flux and concentration (and there is no point sampling when there is no flux). Moreover, in order to justify whether the results above remain consistent for the two stations within each analysis, we randomly resample 1000 times of 66 subset samples again at different time intervals (i.e., opportunistic sampling).
6 4 2 0 0.2
0.4
0.6
0.8
1.0
2 2 Adjusted R for Station05586100 (69264km )
Fig. 6. Adjusted R2 for four different regression models (opportunistic sampling) based on the relationship between concentration and flow rate. . Rating curve: black dashed line ( ), seven-parameter model: red solid line ( ), M1: green dot dash line ( ), M2: blue dotted line ( ).
3.2.3. Estimation of loads Both M1 and M2 appear to have improved the traditional regression methods (Rating curve and seven-parameter model) for both stations. (Fig. 7). For Station 05532500, all regression models underestimate the true total sediment loads, the averages of b L rc , b L 7p , b L M1 , and b L M2 are smaller but their standard errors are larger than those using uniform sampling are. Meanwhile, for Station 05586100, rating curve and seven-parameter model, which overestimates true total sediment loads, yield larger estimation results (average and standard error) than those using uniorm sampling. M1 and M2, which slightly overestimate the total sediment loads, produce larger average estimates of total sediment loads and larger standard errors than those using uniform sampling do. A summary of the bias and standard error of estimation of four regression models for two sampling designs in the two stations is given (Table 2). The biases of two estimators (b L rc and b L 7p ) are relatively larger than estimator ðb L M1 Þ but much larger than estimator ðb L M2 Þ at Station 05532500 regardless of sampling designs. On the other hand, at Station 05586100, the bias of estimator ðb L rc Þ is substantially larger in uniform sampling (65%), and even more so in opportunistic sampling (83%), whereas the biases of estimators b L M1 and b L M2 are much smaller in two sampling designs. Therefore, M1 and M2 provide improvement in estimation of total sediment loads based on two different sampling designs. Employing adjusted R2 as the criterion to choose the appropriate predictive model, it is imperative to construct the relationship between adjusted R2 and biases of estimation of total sediment loads (Fig. 8). As illustrated in Fig. 8, it is clear that the averages of adjusted R2 in M1 and M2 are much higher, with corresponding lower average of biases in two stations regardless of sampling designs. When considering which regression model performs best overall, we note that, for all four regression models, the biases of estimating total load and standard error values obtained by uniform sampling are smaller than those obtained by opportunistic sampling are, whereas the adjusted R2 values obtained by uniform
76
Y.-G. Wang, T. Tian / Journal of Hydrology 482 (2013) 69–78
Density
Station 05532500 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0
(a)
1.5
2.0
2.5
3.0
3.5
(c)
2.0
(b)
μ=2.01 σ=0.12
4.0
4.5
1.5
μ=2.16 σ=0.16
2.0
2.5
3.0
3.5
(d)
μ=2.19 σ=0.16
4.0
4.5
μ=2.40 σ=0.23
1.5 1.0 0.5 0.0 1.5
2.0
2.5
3.0
3.5
4.0
4.5
1.5
2.0
2.5
3.0
3.5
4.0
4.5
Estimated Load (104 tonnes)
Station 05586100
(a)
0.8
(b)
μ=4.45 σ=0.44
0.6
μ=3.76 σ=0.61
0.4
Density
0.2 0.0 1
2
3
4
5
(c)
0.6
6
7
1
2
3
4
5
(d)
μ=3.27 σ=0.56
6
7
μ=2.78 σ=0.58
0.4 0.2 0.0 1
2
3
4
5
7
6
1
2
3
4
5
6
7
Estimated Load (106 tonnes) Fig. 7. Density plot of the estimators (opportunistic sampling) based on rating curve (a), seven-parameter model (b), M1 (c), and M2 (d) of total loads for the two stations in 2005. Dash vertical line: True total sediment loads. Top panel: Station 05532500 (Des Plaines River), Bottom panel: Station 05586100 (Illinois River).
Table 2 Comparison of uniform sampling with opportunistic sampling for the two stations for average bias and SE in a total of 66 subset samples. Estimators
a b
Station 05532500 Des Plaines River
Station 05586100 Illinois River
Uniform Sampling
Opportunistic Sampling
Uniform Sampling
Opportunistic Sampling
biasa
SEb
bias
SE
bias
SE
bias
SE
b L rc b L 7p
–27%
1284
–30%
1183
65%
391135
83%
439861
–21%
1594
–25%
1634
33%
460208
55%
606043
b L M1 b L M2
–20%
1518
–24%
1559
30%
515477
35%
560574
–15%
2205
–16%
2321
9%
443282
14%
582739
Average bias in percentage. Standard error.
sampling are smaller than those obtained by opportunistic sampling are. Also, the flow rate is heteroskedastic and thus using the opportunistic sampling is reasonable and practically desirable. However, more work is needed for efficient sampling design. With the opportunistic sampling, M2 provides the best predictive concentration and the best estimates of total loads, followed by M1. All analysis was performed using statistical package R (free download from http://www.r-project.org/). Our computing codes are available upon request.
4. Discussion Regression approaches to estimation of total nutrient loads based on regression models is the most widely used method for predicting concentration data, which is required when there are less frequent measurements. Thus, better understanding of the accuracies of loads based on curve-fitting models has important implications in setting and evaluating targets in management plans. Incorporating more predictive variables in the regression
77
Y.-G. Wang, T. Tian / Journal of Hydrology 482 (2013) 69–78
Uniform Sampling
Average Bias
−0.10
(a)
0.8
−0.15
M2
0.6
−0.20
(b) rating curve
0.4
M1
7 parameter model
7 parameter model
M1
0.2
−0.25
M2
0.0
rating curve
−0.30 0.4
0.5
0.6
0.7
Average Adjusted R
0.4
0.8
2
0.5
0.6
0.7
Average Adjusted R
0.8
2
Opportunistic Sampling
Average Bias
−0.10
(a)
0.8
−0.15
(b)
7 parameter model
0.6
M2
−0.20
rating curve
0.4 M1
−0.25
M1
0.2
M2
7 parameter model
0.0
rating curve
−0.30 0.4
0.5
0.6
0.7
Average Adjusted R
2
0.8
0.4
0.5
0.6
0.7
Average Adjusted R
0.8
2
2
Fig. 8. Plot of average adjusted R and average bias of load estimates (opportunistic sampling) based on rating curve, seven-parameter model, M1, and M2 for the two stations in 2005 ((a) Des Plaines River and (b) Illinois River). Top panel: uniform sampling, Bottom panel: opportunistic sampling. Rating curve : black solid circle ( ), sevenparameter model: red diamond ( ), M1: green triangle point-up ( ), M2: blue square ( ). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
will produce more accurate concentration prediction. For example, recently, Kuhnert et al. (2012) considered vegetation cover together with discounted flow terms suggested by Wang et al. (2011) and obtained much improved model performance. Load estimates using any regression models will produce biases because many important factors are unknown and real data in practice are never generated from a regression model. Various methods are developed in the literature to correct for biases (Wang et al., 2011; Verma et al., 2012). In this paper, we have emphasized the development of regression models that produce more accurate estimation of both concentration and sediment loads. In particularly, adding a meaningful local predictive variable (ADF) reduces the model errors. As proposed by Kulasova et al. (2012), instantaneous load needs to be predicted that rely on continuous estimates of concentration over time as well as the continuous flux to obtain accurate estimate of annual load and the associated uncertainty. The proposed average discounted flow capturing the past flow patterns appears to be useful in improving the concentration prediction and hence the total annual load. The two basins studied here are immensely different in sizes. Our limited experience is consistent to the observed pattern here – the ADF works more effectively for large catchment basins. This may be due to more stable exhaustion/accumulation process of the concentration compounds. In our data analysis we have intensive flow samples spanning 3 years, and concentration sampling times for total suspended sediment (TSS) random selected 66 samples throughout that period. It is apparent that such selected flow rate data are biased. The time
period over which flow rate and corresponding concentration are collected does influence the accuracy of estimation of load values. Consequently, we employ two different sampling designs (uniform sampling and opportunistic sampling) to compare regression models’ performance. For both Station 05532500 and Station 05586100, M1 and M2 provide more accurate estimation of total load with. Similarly, both M1 and M2 produce higher adjusted R2 compared to other estimators. Moreover, there seems a relationship between adjusted R2 and estimation bias in our illustration, where the higher the adjusted R2 the more accurate the estimation of total load is. Thus, the adjusted R2, as the traditional model selection criteria in statistics or the Nash–Sutcliffe model efficiency in hydrology, is useful to choose appropriate models. We have presented results based on fortnightly and monthly sampling frequencies (with more samples in wet season). The Monte Carlo results using other sampling frequencies can also be easily obtained. For example, we also investigated weekly sampling schemes and have obtained similar conclusions. This is as one might expect because all the estimators here are regressionbased. Based on our studies, we can conclude that the new proposed procedures can significantly improve the predictability especially for large rivers. How much improvement by using the ADF will be site-dependent and the optimal discount factor should be site-specific reflecting the local environmental conditions. It would be exciting to see how the procedures work under different environmental conditions. Another fundamental challenge encountered in estimating sediment total load is the proper uncertainty analysis, such as model
78
Y.-G. Wang, T. Tian / Journal of Hydrology 482 (2013) 69–78
error, of the predictive model of concentration. To effectively investigate uncertainty in loads over a year, as described by Rustomji and Wilkinson (2008), using the bootstrap technique to resample the residuals produces reliable estimates of total sediment loads. Also, to compare the accuracies of sediment total loads, employing adjusted R2 as a reflection of regression model’s variation (i.e., how well a model describes the data) is a basis for developing the predictive model of concentration. On the other hand, improved surrogate technologies, such as optics and acoustics, are available to collect highly resolved time series of sediment concentration which will lead to more accurate estimates of annual load. Improved methods incorporating spatial–temporal dynamics will lead to better governmental decisions in managing costal and ecosystem health (Brodie et al., 2012). Acknowledgements We wish to thank two referees for their very constructive comments and detailed suggestions which led to a much improved paper. References Brodie, J.E., Kroon, F.J., Schaffelke, B., Wolanski, E.C., Lewis, S.E., Devlin, M.J., Bohnet, I.C., Bainbridge, Z.T., Waterhouse, J., Davis, A.M., 2012. Terrestrial pollutant runoff to the Great Barrier Reef: an update of issues, priorities and management responses. Mar. Pollut. Bull. 65, 81–100. Brown, H.E., Hansen, E.A., Champagne Jr., N.E., 1970. A system for measuring total sediment yield from small watersheds. Water Resour. Res. 6 (3), 818–826. Cohn, T.A., 1995. Recent advances in statistical methods for the estimation of sediment and nutrient transport in rivers. Rev. Geophys. 33, 1117–1123. Cohn, T.A., Caulder, D.L., Gilroy, E.J., Zynjuk, L.D., Summers, R.M., 1992. The validity of a simple statistical model for estimating fluvial constituent loads: an empirical study involving nutrient loads entering Chesapeake bay. Water Resour. Res. 28 (9), 2353–2363. Colby, B.R., 1956. The Relationship of Sediment Discharge to Streamflow, Rep., U.S. Geol. Surv., Reston, Va, p. 170. Darnell, R., Henderson, B., Kroon, F.J., Kuhnert, P., 2012. Statistical power of detecting trends in total suspended sediment loads to the Great Barrier Reef. Mar. Pollut. Bull. 65, 203–209. Davis, A.M., Lewis, S.E., Bainbridge, Z.T., Glendenning, L., Turner, R.D., Brodie, J.E., 2012. Dynamics of herbicide transport and partitioning under event flow
conditions in the lower Burdekin region, Australia. Mar. Pollut. Bull. 65, 182– 193. Ferguson, R.I., 1986. River loads underestimated by rating curves. Water Resour. Res. 22 (1), 74–76. Koch, R.W., Smillie, G.M., 1986. Bias in hydrologic prediction using log-transformed regression models. Water Resour. Bull. 22 (5), 717–723. Kuhnert, P.M., Henderson, B.L., Lewis, S.E., Bainbridge, Z.T., Wilkinson, S.N., Brodie, J.E., 2012. Quantifying total suspended sediment export from the Burdekin River catchment using the loads regression estimator tool. Water Resour. Res. 48, W04533. Kulasova, A., Smith, P.J., Beven, K.J., Blazkova, S.D., 2012. A method of computing uncertain nitrogen and phosphorus loads in a small stream from an agricultural catchment using continuous monitoring data. J. Hydrol. 458–459, 1–8. Miller, C.R., 1951. Analysis of Flow-duration, Sediment-Rating Curve Method of Computing Sediment Yield, Rep., U.S. Bur. of Reclam., Denver, Colo, p. 15. Rustomji, P., Wilkinson, S., 2008. Applying bootstrap resampling to quantify uncertainty in fluvial suspended sediment loads estimated using rating curve. Water Resour. Res. 44, W09434. http://dx.doi.org/10.1029/2007WR006088. Singer, M.B., Dunne, T., 2001. Identifying eroding and depositional reaches of valley by analysis of suspended sediment Transport in the Sacramento River, California. Water Resour. Res. 37 (12), 3371–3381 (Identifying eroding and depositional reaches of valley by analysis of suspended sediment Transport in the Sacramento River, California). Thomas, R.B., 1985. Estimating total suspended sediment yield with probability sampling. Water Resour. Res. 21 (9), 1381–1388. Thomas, R.B., 1988. Monitoring baseline suspended sediment in forested basins: the effects of sampling of suspended sediment rating curves. Hydrol. Sci. 33 (5), 499–514. Thomas, R.B., Lewis, J., 1993. A comparison of selection at list time and timestratified sampling for estimating suspended sediment loads. Water Resour. Res. 29, 1247–1256. Verhoff, F.H., Yaksich, S.M., Melfi, D.A., 1980. River nutrient and chemical transport estimation. J. Environ. Eng. Div. Am. Soc. Civ. Eng. 106, 591–608. Verma, S., Markus, M., Cooke, R.A., 2012. Development of error correction techniques for nitrate-N load estimation. J. Hydrol. 432, 12–25. Wallace, J., Stewart, L., Hawdon, A., Keen. R., 2008. The role of coastal floodplains in generating sediment and nutrient fluxes to the great barrier reef lagoon in Australia. In: Ecohydrological Processes and Sustainable Floodplain Management Opportunities and Concepts for Water Hazard Mitigation and Ecological and Socioeconomic Sustainability, Lodz, Poland. Walling, D.E., 1977. Assessing the accuracy of suspended sediment rating curves for a small basin. Water Resour. Res. 13 (3), 531–538. Wang, P., Linker, L.C., 2008. Improvement of regression simulation in fluvial sediment loads. J. Hydraul. Eng. 134 (10), 1527–1531. Wang, Y.-G., Kuhnert, P., Henderson, B., 2011. Load estimation with uncertainties from opportunistic sampling data – a semiparametric approach. J. Hydrol. 396, 148–157.
Contents lists available at SciVerse ScienceDirect
Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol
Sediment concentration prediction and statistical evaluation for annual load estimation You-Gan Wang a,⇑, Ting Tian b a b
Centre for Applications in Natural Resource Mathematics (CARM), School of Mathematics and Physics, The University of Queensland, St. Lucia 4072, Australia School of Agriculture and Food Science, The University of Queensland, St. Lucia 4072, Australia
a r t i c l e
i n f o
Article history: Received 22 September 2012 Received in revised form 19 December 2012 Accepted 19 December 2012 Available online 8 January 2013 This manuscript was handled by Geoff Syme, Editor-in-Chief Keywords: Bootstrap Load estimation Nash–Sutcliffe model efficiency coefficient Rating curve Suspended sediment Uncertainty
s u m m a r y We consider the development of statistical models for prediction of constituent concentration of riverine pollutants, which is a key step in load estimation from frequent flow rate data and less frequently collected concentration data. We consider how to capture the impacts of past flow patterns via the average discounted flow (ADF) which discounts the past flux based on the time lapsed – more recent fluxes are given more weight. However, the effectiveness of ADF depends critically on the choice of the discount factor which reflects the unknown environmental cumulating process of the concentration compounds. We propose to choose the discount factor by maximizing the adjusted R2 values or the Nash–Sutcliffe model efficiency coefficient. The R2 values are also adjusted to take account of the number of parameters in the model fit. The resulting optimal discount factor can be interpreted as a measure of constituent exhaustion rate during flood events. To evaluate the performance of the proposed regression estimators, we examine two different sampling scenarios by resampling fortnightly and opportunistically from two real daily datasets, which come from two United States Geological Survey (USGS) gaging stations located in Des Plaines River and Illinois River basin. The generalized rating-curve approach produces biased estimates of the total sediment loads by 30% to 83%, whereas the new approaches produce relatively much lower biases, ranging from 24% to 35%. This substantial improvement in the estimates of the total load is due to the fact that predictability of concentration is greatly improved by the additional predictors. Crown Copyright Ó 2013 Published by Elsevier B.V. All rights reserved.
1. Introduction The estimation of fluvial constituent loads is an essential for helping researchers and water managers understand water quality and its effect on environmental resources. Much work on estimation of constituent loads focuses on the ‘‘rating curve’’ approach on the assumption that there is a linear relationship between concentration of sediment and flow on log transformed space (Brown et al., 1970; Colby, 1956; Miller, 1951; Thomas, 1985, 1988; Verhoff et al., 1980; Walling, 1977). Generally, the linear regression analysis uses the log-transformed flow rate as the predictor. However, constituent concentrations often depend on many other environmental changes such as ‘‘first flush’’ or ‘‘depletion’’ effects. The first flush is ‘‘the first significant channelized flow of the wet season which is generally accompanied by relatively high sediment and nutrient concentrations (Wang et al., 2011). Such dynamics is also well discussed in Davis et al. (2012). There are higher flows in wet seasons, and the start of a wet season may change from year to year which makes stratification difficult (Wallace et al., 2008). High nutrient concentrations
⇑ Corresponding author. Tel.: +61 7 3365 2311. E-mail address: [email protected] (Y.-G. Wang).
often come with high flux initially, but not sustainable due to the exhaustion and slow accumulation process of sediments. Simple linear or periodic mathematical functions therefore often yield poor prediction. Therefore, seeking additional hydrologically meaningful predictors or directly modeling the physical process is fundamental in improving load estimation (Wang et al., 2011). Limitations of the rating curve method were well studied by Ferguson (1986). Since then, new statistical methods have been developed expressing the relationship between constituent concentration and flow. One is the seven-parameter log linear model for concentration, which generalizes the rating curve by incorporating five additional predictors: a quadratic term of the logarithm of flow rate, and four parametric functions for seasonality (sine and cosine functions for annual and semi-annual cycles) (Cohn et al., 1992). This pioneering work has provided fundamental guidelines for constructing an appropriate regression model for load estimation based on frequently collected flow data but less-frequently collected concentration data. Additional predictive covariates to the original 7-parameter model can better account for hysteresis in sediment transport (Wang and Linker, 2008). Integration of flow history was proposed by Singer and Dunne (2001) using ARIMA models. Such time series models can extract information from temporal correlations. This is useful if the exhaustion is very fast leading to short lag effects (a few days). Note that one parameter
0022-1694/$ - see front matter Crown Copyright Ó 2013 Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jhydrol.2012.12.043
70
Y.-G. Wang, T. Tian / Journal of Hydrology 482 (2013) 69–78
is needed for each lag and overfit can become a concern if too many parameters are used. Stratification by season or region can also improve predictions. However, stratification can be quite subjective, and specific knowledge is required to determine the strata (Thomas and Lewis, 1993). One can also consider incorporating strata heterogeneity into a regression framework by creating indicator/dummy variables. As pointed out by Wang et al. (2011), simple linear or periodic time functions will not adequately capture information such as a ‘‘flush’’, and more novel covariates characterizing these phenomena will need to be considered. Wang et al. (2011) proposed a new predictive variable, namely, the average discounted flow (ADF). This aims to provide appropriate weights for all the past flow records in a discounted way because more recent fluxes may exhaust more sediment, whereas flux occurring a long time ago is given less weight. There are increasing concerns over quality of water, leading to a demand for accurate estimates of total sediment loads. There are
(a)
1200
Concentration (mg/L)
Concentratation (mg/L)
300
several methods available to estimate total loads from frequent flow rate data and concentration data; however, frequency of concentration data is often very low, and hence concentration prediction becomes a key step in load estimation (Darnell et al., 2012). Estimation of total sediment loads therefore becomes challenging. Existing rating curve and 7-parameter methods may produce biased load estimates, and multiple regression approach appears to be generally more versatile and easy to incorporate additional covariates. Both Koch and Smillie (1986) and Ferguson (1986) found that the regression model (7-parameter) under-predicts concentrations. However, such bias due to nonlinear transformation is often small relative to biases from other sources. Recently, Verma et al. (2012) also investigated bias and standard errors for three estimators and concluded importance of accurate estimators for cost-effective monitoring designs. In this paper we focus on all models for predicting concentration and evaluate their performance in terms of R2 and the resulting load estimates. In particular, we define the new predictive variable (i.e. ADF), which can mimic
250 200 150 100 50
(b)
1000 800 600 400 200 0
0 03 Apr Jul Oct 04 Apr Jul Oct 05 Apr Jul
03 Apr Jul Oct 04 Apr Jul Oct 05 Apr Jul
Date
Date
(c)
80
Flow rate (CM/s)
Flow rate (CM/s)
4 3 2
(d)
60 40
1
20
0
0 03 Apr Jul Oct 04 Apr Jul Oct 05 Apr Jul
03 Apr Jul Oct 04 Apr Jul Oct 05 Apr Jul
Date
Date
(e)
140
Load (1000 tonnes/day)
Load (tonnes/day)
2500 2000 1500 1000 500
(f)
120 100 80 60 40 20 0
0 03 Apr Jul Oct 04 Apr Jul Oct 05 Apr Jul
Date
03 Apr Jul Oct 04 Apr Jul Oct 05 Apr Jul
Date
Fig. 1. The trends of concentration, flow rate, and load for two stations between 2003 and 2005. (a), (c), (e) Station 05532500 (Des Plaines River), (b), (d), (f) Station 05586100 (Illinois River).
71
Y.-G. Wang, T. Tian / Journal of Hydrology 482 (2013) 69–78
models in predicting the annual cumulative load, each regression model using two different sampling designs (uniform sampling and opportunistic sampling) will be described in Section 3. We subsample from the 3 years’ daily records for sediment, using uniform (once every fortnight) and opportunistic sampling (more frequent in wet periods than in dry periods) and use the adjusted R2 as a criterion to determine the most appropriate concentration predictive model. In Section 4, we present our findings; then, a discussion of the methods and possible generalization will be presented.
the recovery and exhaustion process commonly seen in temporal measures of concentration and flow; also, an optimal discount factor is sought from the data. Unlike the time series models, the number of additional parameters does not change even if exhaustion process is slow (but require discount factor close to 1). The paper is organized as follows: in Section 2 we will briefly introduce the two datasets from USGS (see http://il.water.usgs.gov/ for details about the profiler). The new models we propose are presented and discussed. To evaluate the performance of regression
2. Materials and methods Table 1 The mean values of flow (CM/s), concentration (C, mg/L) and load (tonnes/day) for the period of January 2003 to September 2005. Station 05532500 Des Plaines River
2.1. Materials
Station 05586100 Illinois River
Period
Flow
C
Load
Flow
C
Load
2003.01–2003.04 2003.05–2003.08 2003.09–2004.04 2004.05–2004.08 2004.09–2005.04 2005.05–2005.08 2005.09–2005.09
10 21 17 33 17 7 5
37 47 48 56 57 34 37
51 122 104 234 136 35 23
273 634 479 737 898 230 178
128 256 159 205 165 48 33
3963 20046 9626 16523 15310 1242 557
The data considered in this study were collected daily from 2003 to 2005 by US Geological Survey. Two stations will be considered: Station (05532500) of the Des Plaines River at Riverside, located at latitude (41°4900 200 ) and longitude (87°4900 150 ), and Station (05586100) of the Illinois River at Valley City, located at latitude (39°4200 120 ) and longitude (90°3800 430 ). Drainage area of the former Station is 1632 square kilometers (km2) and that of the latter is 69,264 square kilometers (km2). Fig. 1 shows the three time series plots (daily measurements of flow rate, concentration and load from January 1, 2003 to Septem-
Station 05532500 5.5 5.0
(a)
(b)
True Concentration
4.5 4.0 3.5 3.0 2.5 5.5 5.0
3.0
3.5
4.0
4.5
5.0
5.5
6.0
(c)
2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
3.0
3.5
4.0
4.5
5.0
5.5
6.0
(d)
4.5 4.0 3.5 3.0 2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
2.5
Predicted Concentration
True Concentration
Station 05586100 6.5 6.0 5.5 5.0 4.5 4.0 3.5 3.0
(b)
(a)
3 6.5 6.0 5.5 5.0 4.5 4.0 3.5 3.0
4
5
6
3
4
5
6
3
4
5
6
(d)
(c)
3
4
5
6
Predicted Concentration Fig. 2. Plot of predicted concentrations based on four regression models (uniform sampling) versus the true concentration for two stations during 2005 (Log Scale Space). (a) Rating curve: black solid circle ( ), (b) seven-parameter model: red diamond ( ), (c) M1: green triangle point-up ( ), (d) M2: blue square ( ). Solid line: predicted concentration equal to true concentration.
72
Y.-G. Wang, T. Tian / Journal of Hydrology 482 (2013) 69–78
ber 30, 2005) for the two stations (05532500 and 05586100). The concentration data in Station 05532500 shows a periodic pattern, whereas the concentration in Station 05586100 shows two global peaks in the middles of year 2003 and year 2004. It seems concentration may be seasonal, which can be modeled by periodic function. The flow rate in each station reaches annual peaks in the middle of 2003 and 2004, but reaches the annual peak at the beginning of 2005. Similarly, daily suspended sediment loads reach a local peak in the middles of 2003 and 2004, but reaches the annual peak at the beginning of 2005 in each station. 2.2. Methods Suppose the flow rate Q(t) (e.g., in L/s) and concentration C(t) (e.g., in mg/L) are given, then the total load L can be written as the following equation:
Z
L¼K
T
CðtÞQðtÞdt;
ð1Þ
0
where K is a unit-conversion constant depending on the units applied, T is time measured in years. In practice, the samples are only collected at discrete time intervals (daily, weekly or even longer intervals). As a result, the total load L in Eq. (1) may be approximated by
LK
M X
C m Q m d;
ð2Þ
models only make use of the current flow rate and ignore the past rates. To this end, Wang et al. (2011) introduced ADF. This new variable is constructed using weighted (discounted) mean of the previous flux. The discount is based on the time lapsed since the flow. The discount factor can be so chosen to optimize the model fit (such as coefficient of determination or R2 value). Note that the R2 is also known as the Nash–Sutcliffe model efficiency in hydrology. This modeling approach greatly improved the concentration prediction model for the Tully River in the Great Barrier Reef catchments, Australia. The idea is to employ exponential discounting to all the flow historic data. Suppose (Q1, Q2, . . ., Qj) are the flow rates from day 1 to day j. The ADF with discount factor d at day j is defined by the following equation:
ADFðdÞ ¼
b L¼K
M X
^ m Q d; C m
j1
j1
1 þ d þ þ d
Q1
ð7Þ
;
where discount factor d can be measured from 0.1 to 1. This Eq. (7) summarizes past flow rates (Q1, Q2, . . ., Qj) as an index by giving more weight to more recent events. Roughly, d can be interpreted as the exhaustion rate. We propose to find the optimal d in terms of adjusted R2 values by taking account of the number of parameters used. The ADF will be used to replace the current flow rate Q(t). We first introduce our proposed new regression model, M1, within the extended seven-parameter model:
lnðCÞ ¼ b0 þ b1 lnðADFÞ þ b2 lnðQ Þ þ b3 ½lnðQ Þ2 þ b4 sinð2pTÞ þ b5 cosð2pTÞ þ b6 sinð4pTÞ þ b7 cosð4pTÞ þ e
m¼1
where d is the time interval and dM = T, Qm and Cm are the true values of flow and concentration, respectively, at time tm (the beginning of the mth time interval). However, the concentration data is not always available, so the concentration has to be predicted (at small intervals, e.g. every 10 min.) from a generalized rating curve (an empirical relationship between C and Q); the total load L can be estimated by
Q j þ dQ j1 þ þ d
ð8Þ
The optimal d is the one that produces the highest adjusted R2 values. It is possible that an improved model would employ two ADF variables with two different discount rates, d1 and d2:
lnðCÞ ¼ b0 þ b1 lnðADF 1 Þ þ b2 lnðADF 2 Þ þ b3 ½lnðADF 1 Þ2 þ b4 ½lnðADF 2 Þ2 þ b5 sinð2pTÞ þ b6 cosð2pTÞ þ b7 sinð4pTÞ þ b8 cosð4pTÞ þ e
ð9Þ
ð3Þ
m¼1
b L¼K
2 b i Q exp s ; C i 2 i¼1
M X
10
rating curve seven−parameter model M1 M2
8
Density
where the unknown concentration of Cm in Eq. (2) is substituted b m , which is often predicted from a statistical using estimates of C model such as rating curve or other regression model. Generally, b m is predicted on log scale, and the bias corrected total load estiC mator is
6 4
ð4Þ
2
where s2 is the estimated variance of the model error on log scale (see Ferguson, 1986; Koch and Smillie, 1986; Cohn, 1995).
0 0.2
0.4
0.6
0.8
2
1.0 2
Adjusted R for Station05532500 (1632km ) 2.2.1. Regression models We now briefly introduce two widely used regression models. The rating curve approach is based on the following regression model:
ð5Þ
where ln() denotes the natural logarithm function, b0 and b1 are parameters to be estimated, and e is random error, which is identical independent distributed (i.i.d). A total of k = 7 parameters are supposed to be required to formulate the relationship between concentration and flow rate:
6 4 2 0 0.2
0.4
0.6
0.8
1.0
Adjusted R2 for Station05586100 (69264km2)
2
lnðCÞ ¼ b0 þ b1 lnðQ Þ þ b2 ½lnðQ Þ þ b3 sinð2pTÞ þ b4 cosð2pTÞ þ b5 sinð4pTÞ þ b6 cosð4pTÞ þ e;
8
Density
lnðCÞ ¼ b0 þ b1 ln Q þ e;
rating curve seven−parameter model M1 M2
10
ð6Þ
where the errors, e Nð0; r2e Þ, and b are parameters of the model that can be estimated from the data. The above two regression
Fig. 3. Comparisons of adjusted R2 from four different regression models (uniform sampling) based on the relationship between concentration and flow rate. Rating curve: black dashed line ( ), seven-parameter model: red solid line ( ), M1: green dot dash line ( ), M2: blue dotted line ( ).
73
Y.-G. Wang, T. Tian / Journal of Hydrology 482 (2013) 69–78
If we restrict d1 close to 0, ADF1 will become Q. Using multiple ADF variables will capture exhaustion processes with different cycles. Again, we will choose the optimal d values to maximize the adjusted R2. We will refer to this model as M2. The adjusted R2 values can be obtained for a range of discount values: 0, 0.1, 0.2, . . ., 0.9, 0.95, 0.96, . . ., 0.99, and hence the optimal d can be chosen. 2.2.2. Resampling In order to compare the four load estimation methods based on above regression models, we apply two sampling designs, i.e. uniform sampling and opportunistic sampling. Uniform sampling randomly selects one sample from every fixed period of two weeks (sampling interval is not fixed but has a mean of 2 weeks), while opportunistic sampling randomly selects one sample from a shorter (1 week) or longer (1 month) period depending on wet or dry seasons. The full data set incorporates 1004 daily flow rate observations and corresponding concentration data, covering from 1st of January for 2003 to 30th of September for 2005. To achieve a moderate annual sample size (20–30), we will examine fortnightly and monthly sampling schemes. Table 1 shows the mean values for the selected wet and dry seasons. For illustrative purpose, we will consider random uniform sampling and opportunistic sampling using a subjectively defined ‘‘wet/dry’’ seasons. The uniform random samples are selected every fortnight (resulting in 66 samples during these 3 years). For opportunistic sampling, samples are se-
lected once a month in the dry season (from September to next May) and once a week in the wet season (between June and August). As a result, the total sample size using opportunistic sampling is the same as that using the uniform sampling. These two sampling designs have a total sample size of 66 in each simulation. Note that the resampling process is repeated for many times and overall the resampled data will preserve both the variability of original flow data (albeit at different frequency) and their relationship with corresponding concentration. The two stations are located on different rivers, where dry and wet seasons varies over years as shown in Fig. 1. This will allow us to demonstrate how the two sampling designs might differ. In terms of sampling frequency, weekly samples can also be examined. The estimators discussed here are all regression-based, unlike other imputation or composite methods, relative performance, will not so depend on the sample frequency. 3. Results 3.1. Uniform sampling 3.1.1. Concentration prediction based on four regression models The performance of predicted daily concentration where the unit for daily concentration is mg/L in the log-space in the two stations is depicted in Fig. 2. For each predictive regression model, the
Density
Station 05532500 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0
(a)
1.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0
2.0
2.5
3.0
(c)
1.5
(b)
μ=2.10 σ=0.13
3.5
4.0
1.5
2.5
3.0
3.5
2.0
2.5
3.0
(d)
μ=2.27 σ=0.15
2.0
μ=2.25 σ=0.16
4.0
1.5
3.5
4.0
μ=2.45 σ=0.20
2.0
2.5
3.0
3.5
4.0
Estimated Load (104 tonnes)
Density
Station 05586100 1.2 1.0 0.8 0.6 0.4 0.2 0.0
(a)
1 1.0
2
3
4
(c)
0.8
(b)
μ=4.00 σ=0.39
5
6
1
μ=3.22 σ=0.46
2
3
4
(d)
μ=3.14 σ=0.52
5
6
μ=2.65 σ=0.44
0.6 0.4 0.2 0.0 1
2
3
4
5
6
1
2
3
4
5
6
Estimated Load (106 tonnes) Fig. 4. Density plots of the four estimators (uniform sampling) based on rating curve (a), seven-parameter model (b), M1 (c), and M2 (d) of total loads for the two stations in 2005. Dash vertical line: True total sediment loads. Top panel: Station 05532500 (Des Plaines River), Bottom panel: Station 05586100 (Illinois River).
74
Y.-G. Wang, T. Tian / Journal of Hydrology 482 (2013) 69–78
performance of predictive logarithm concentration seems to be close to the true logarithm concentration in Station 05532500. Typically, when the values of logarithm concentration are below 4, predictive logarithm concentrations significantly match their corresponding logarithm true concentrations. However, the values of predictive logarithm concentration between 3.5 and 4.5 underestimate true logarithm concentration, while they overestimate true logarithm concentration over 5.0 for four regression models in Station 05586100. The reasons or causes for such patterns are unknown implying some other important factors affecting the concentration are to be found. In particular, the rating curve provides larger predictive logarithm concentration at large values of true logarithm concentration. Nevertheless, the values of predictive logarithm concentration for M1 and M2 are relatively closer to the true logarithm concentration compared to with the rating curve and seven-parameter model. Overall, four regression models provide more accurate predictive logarithm concentration at relatively lower true logarithm concentration (below 4.5), and M1 and M2 yield comparatively better results. 3.1.2. Adjusted R2 of prediction concentration based on four regression models Fig. 3 illustrates the effect of average discount flow on the prediction concentration based on log linear flow rate. We randomly
resample the flow rate over the timeframe and the corresponding concentration of two stations 1000 times, and compare the density of each adjusted R2 of four regression models. It is apparent that the adjusted R2 of new regression model one (M1) and new regression model two (M2) can be improved in each station. Choosing the optimal discount factor, d, in M1 (Eq. (8)) and M2 (Eq. (9)) for 66 random subsets of data samples associated with the highest adjusted R2, it is clear that M1 and M2 based on the optimal discount factor, d by adjusted R2, provide higher adjusted R2 than that generalized from the rating curve and seven-parameter model. The adjusted R2 of Station 05532500 ranges for all regression models as follow: rating curve approach (0.18, 0.7), seven-parameter model (0.2, 0.76), M1 (0.25,0.82) and M2 (0.3,0.84), in addition, the peaks of the adjusted R2 for M1 and M2 are higher than those for the rating curve approach (0.43) and seven-parameter model (0.52), increasing to 0.58 and 0.62, respectively. Meanwhile, the adjusted R2 of Station 05586100 ranges from 0.4 to 0.8 for M1 and M2, respectively, which is larger than that obtained from the rating curve approach (0.2, 0.6) and seven-parameter model (0.22, 0.74). Also, the peaks of the adjusted R2 for M1 and M2 are higher than those for the rating curve approach (0.4) and seven-parameter model (0.48), increasing to 0.61 and 0.65, respectively Thus, above results of these two stations indicate that M1 and M2 yield higher adjusted R2.
Station 05532500 5.5 5.0
(a)
(b)
True Concentration
4.5 4.0 3.5 3.0 2.5 5.5 5.0
3.0
3.5
4.0
4.5
5.0
5.5
6.0
(c)
2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
3.0
3.5
4.0
4.5
5.0
5.5
6.0
(d)
4.5 4.0 3.5 3.0 2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
2.5
Predicted Concentration
True Concentration
Station 05586100 6.5 6.0 5.5 5.0 4.5 4.0 3.5 3.0
(b)
(a)
3 6.5 6.0 5.5 5.0 4.5 4.0 3.5 3.0
4
5
6
3
4
5
6
3
4
5
6
(d)
(c)
3
4
5
6
Predicted Concentration Fig. 5. Plotof predicted concentration based on four regression models (opportunistic sampling) versus the true concentration for the two stations during 2005 (Log Space). (a) Rating curve: black solid circle ( ), (b) seven-parameter model: red diamond ( ), (c) M1: green triangle point-up ( ), M2: (d) blue square ( ).Solid line: predicted concentration equal to true concentration.
75
Y.-G. Wang, T. Tian / Journal of Hydrology 482 (2013) 69–78
3.2. Opportunistic sampling 3.2.1. Concentration prediction based on the four regression models As shown in Fig. 5, it is clear that, by means of opportunistic sampling, there is a similar performance of predictive logarithm concentration in Station 05532500 compared to the uniform sampling approach (upper panel of Fig. 2). However, the performance of predictive logarithm concentration at Station 05586100 is very different from that based on uniform sampling (lower panel of Fig. 2). There are more variations when the values of logarithm concentration exceed 5.0 mg/L. Typically, rating curve and sevenparameter models provide larger values of predictive logarithm concentration than M1 and M2 do at high concentrations. Moreover, the ranges of predicted concentration (on logarithm) by the rating curve and the seven-parameter model are wider, in which the upper bounds of predicted concentration are as high as 5.5 and 6.0, respectively. The values of predictive logarithm concentration overestimate true values for all regression models at their lower bounds. 3.2.2. Adjusted R2 of prediction concentration based on four regression models Compared with Fig. 3, the lower bounds of adjusted R2 for M1 and M2 increase to 0.3 and 0.4, respectively, while the upper bound of adjusted R2 for M1 and M2 increases to 0.85 in Station 05532500 (upper panel of Fig. 6). Corresponding peak values of adjusted R2 for M1 and M2 increase to 0.63 and 0.7, respectively. Moreover, the lower bounds of adjusted R2 of the four regression models slightly increase by 0.2 (the minimum values of adjusted R2 of rating curve and seven-parameter model from 0.2 to 0.4 and the minimum values of adjusted R2 of M1 and M2 from 0.4 to 0.6) in Station 05586100 (lower panel of Fig. 6). The peaks of adjusted R2 values of M1 and M2 increase to around 0.72. Overall, the values of adjusted R2 of the four regression models increase for opportunistic sampling, and there is an increase of adjusted R2 in M1 and M2.
rating curve seven−parameter model M1 M2
10
Density
8 6 4 2 0 0.2
0.4
0.6
0.8
2
1.0 2
Adjusted R for Station05532500 (1632km ) 10
rating curve seven−parameter model M1 M2
8
Density
3.1.3. Estimation of loads Estimates of total sediment loads for two stations in 2005 include the estimator based on rating curve ðb L rc Þ, the estimator based on seven-parameter model ðb L 7p Þ, the estimator based on M1 ðb L M1 Þ, and the estimator based on M2 ðb L M2 Þ. The estimates of total sediment load (Eq. (4)) based on four regression models of the two stations are shown in Fig. 4. We randomly resample 1000 times 66 subset samples to estimate total sediment loads. Although all four estimators of total sediment load underestimate the true load for Station 05532500 (upper panel of Fig. 4.), the averages of b L M1 and b L M1 are closer to true total sediment loads (the biases are 20% and 15%, while the bias in b L rc and b L 7p are 27% and 21%, respectively). On the other hand, b L rc and b L 7 produce overestimates, whereas b L M1 and b L M2 produce much less biased estimates especially for the second Station (Fig. 4). The rating curve produces larger estimates than the seven-parameter model, M1, and M2. Comparing four estimators in two stations, M1 and M2, regardless of whether the year is relatively dry or wet, provide better results on estimation of total sediment loads. In reality, we often intend to sample more often when flux is large. This opportunistic sampling makes sense for more reliable load estimation because the load is a product of flux and concentration (and there is no point sampling when there is no flux). Moreover, in order to justify whether the results above remain consistent for the two stations within each analysis, we randomly resample 1000 times of 66 subset samples again at different time intervals (i.e., opportunistic sampling).
6 4 2 0 0.2
0.4
0.6
0.8
1.0
2 2 Adjusted R for Station05586100 (69264km )
Fig. 6. Adjusted R2 for four different regression models (opportunistic sampling) based on the relationship between concentration and flow rate. . Rating curve: black dashed line ( ), seven-parameter model: red solid line ( ), M1: green dot dash line ( ), M2: blue dotted line ( ).
3.2.3. Estimation of loads Both M1 and M2 appear to have improved the traditional regression methods (Rating curve and seven-parameter model) for both stations. (Fig. 7). For Station 05532500, all regression models underestimate the true total sediment loads, the averages of b L rc , b L 7p , b L M1 , and b L M2 are smaller but their standard errors are larger than those using uniform sampling are. Meanwhile, for Station 05586100, rating curve and seven-parameter model, which overestimates true total sediment loads, yield larger estimation results (average and standard error) than those using uniorm sampling. M1 and M2, which slightly overestimate the total sediment loads, produce larger average estimates of total sediment loads and larger standard errors than those using uniform sampling do. A summary of the bias and standard error of estimation of four regression models for two sampling designs in the two stations is given (Table 2). The biases of two estimators (b L rc and b L 7p ) are relatively larger than estimator ðb L M1 Þ but much larger than estimator ðb L M2 Þ at Station 05532500 regardless of sampling designs. On the other hand, at Station 05586100, the bias of estimator ðb L rc Þ is substantially larger in uniform sampling (65%), and even more so in opportunistic sampling (83%), whereas the biases of estimators b L M1 and b L M2 are much smaller in two sampling designs. Therefore, M1 and M2 provide improvement in estimation of total sediment loads based on two different sampling designs. Employing adjusted R2 as the criterion to choose the appropriate predictive model, it is imperative to construct the relationship between adjusted R2 and biases of estimation of total sediment loads (Fig. 8). As illustrated in Fig. 8, it is clear that the averages of adjusted R2 in M1 and M2 are much higher, with corresponding lower average of biases in two stations regardless of sampling designs. When considering which regression model performs best overall, we note that, for all four regression models, the biases of estimating total load and standard error values obtained by uniform sampling are smaller than those obtained by opportunistic sampling are, whereas the adjusted R2 values obtained by uniform
76
Y.-G. Wang, T. Tian / Journal of Hydrology 482 (2013) 69–78
Density
Station 05532500 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0
(a)
1.5
2.0
2.5
3.0
3.5
(c)
2.0
(b)
μ=2.01 σ=0.12
4.0
4.5
1.5
μ=2.16 σ=0.16
2.0
2.5
3.0
3.5
(d)
μ=2.19 σ=0.16
4.0
4.5
μ=2.40 σ=0.23
1.5 1.0 0.5 0.0 1.5
2.0
2.5
3.0
3.5
4.0
4.5
1.5
2.0
2.5
3.0
3.5
4.0
4.5
Estimated Load (104 tonnes)
Station 05586100
(a)
0.8
(b)
μ=4.45 σ=0.44
0.6
μ=3.76 σ=0.61
0.4
Density
0.2 0.0 1
2
3
4
5
(c)
0.6
6
7
1
2
3
4
5
(d)
μ=3.27 σ=0.56
6
7
μ=2.78 σ=0.58
0.4 0.2 0.0 1
2
3
4
5
7
6
1
2
3
4
5
6
7
Estimated Load (106 tonnes) Fig. 7. Density plot of the estimators (opportunistic sampling) based on rating curve (a), seven-parameter model (b), M1 (c), and M2 (d) of total loads for the two stations in 2005. Dash vertical line: True total sediment loads. Top panel: Station 05532500 (Des Plaines River), Bottom panel: Station 05586100 (Illinois River).
Table 2 Comparison of uniform sampling with opportunistic sampling for the two stations for average bias and SE in a total of 66 subset samples. Estimators
a b
Station 05532500 Des Plaines River
Station 05586100 Illinois River
Uniform Sampling
Opportunistic Sampling
Uniform Sampling
Opportunistic Sampling
biasa
SEb
bias
SE
bias
SE
bias
SE
b L rc b L 7p
–27%
1284
–30%
1183
65%
391135
83%
439861
–21%
1594
–25%
1634
33%
460208
55%
606043
b L M1 b L M2
–20%
1518
–24%
1559
30%
515477
35%
560574
–15%
2205
–16%
2321
9%
443282
14%
582739
Average bias in percentage. Standard error.
sampling are smaller than those obtained by opportunistic sampling are. Also, the flow rate is heteroskedastic and thus using the opportunistic sampling is reasonable and practically desirable. However, more work is needed for efficient sampling design. With the opportunistic sampling, M2 provides the best predictive concentration and the best estimates of total loads, followed by M1. All analysis was performed using statistical package R (free download from http://www.r-project.org/). Our computing codes are available upon request.
4. Discussion Regression approaches to estimation of total nutrient loads based on regression models is the most widely used method for predicting concentration data, which is required when there are less frequent measurements. Thus, better understanding of the accuracies of loads based on curve-fitting models has important implications in setting and evaluating targets in management plans. Incorporating more predictive variables in the regression
77
Y.-G. Wang, T. Tian / Journal of Hydrology 482 (2013) 69–78
Uniform Sampling
Average Bias
−0.10
(a)
0.8
−0.15
M2
0.6
−0.20
(b) rating curve
0.4
M1
7 parameter model
7 parameter model
M1
0.2
−0.25
M2
0.0
rating curve
−0.30 0.4
0.5
0.6
0.7
Average Adjusted R
0.4
0.8
2
0.5
0.6
0.7
Average Adjusted R
0.8
2
Opportunistic Sampling
Average Bias
−0.10
(a)
0.8
−0.15
(b)
7 parameter model
0.6
M2
−0.20
rating curve
0.4 M1
−0.25
M1
0.2
M2
7 parameter model
0.0
rating curve
−0.30 0.4
0.5
0.6
0.7
Average Adjusted R
2
0.8
0.4
0.5
0.6
0.7
Average Adjusted R
0.8
2
2
Fig. 8. Plot of average adjusted R and average bias of load estimates (opportunistic sampling) based on rating curve, seven-parameter model, M1, and M2 for the two stations in 2005 ((a) Des Plaines River and (b) Illinois River). Top panel: uniform sampling, Bottom panel: opportunistic sampling. Rating curve : black solid circle ( ), sevenparameter model: red diamond ( ), M1: green triangle point-up ( ), M2: blue square ( ). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
will produce more accurate concentration prediction. For example, recently, Kuhnert et al. (2012) considered vegetation cover together with discounted flow terms suggested by Wang et al. (2011) and obtained much improved model performance. Load estimates using any regression models will produce biases because many important factors are unknown and real data in practice are never generated from a regression model. Various methods are developed in the literature to correct for biases (Wang et al., 2011; Verma et al., 2012). In this paper, we have emphasized the development of regression models that produce more accurate estimation of both concentration and sediment loads. In particularly, adding a meaningful local predictive variable (ADF) reduces the model errors. As proposed by Kulasova et al. (2012), instantaneous load needs to be predicted that rely on continuous estimates of concentration over time as well as the continuous flux to obtain accurate estimate of annual load and the associated uncertainty. The proposed average discounted flow capturing the past flow patterns appears to be useful in improving the concentration prediction and hence the total annual load. The two basins studied here are immensely different in sizes. Our limited experience is consistent to the observed pattern here – the ADF works more effectively for large catchment basins. This may be due to more stable exhaustion/accumulation process of the concentration compounds. In our data analysis we have intensive flow samples spanning 3 years, and concentration sampling times for total suspended sediment (TSS) random selected 66 samples throughout that period. It is apparent that such selected flow rate data are biased. The time
period over which flow rate and corresponding concentration are collected does influence the accuracy of estimation of load values. Consequently, we employ two different sampling designs (uniform sampling and opportunistic sampling) to compare regression models’ performance. For both Station 05532500 and Station 05586100, M1 and M2 provide more accurate estimation of total load with. Similarly, both M1 and M2 produce higher adjusted R2 compared to other estimators. Moreover, there seems a relationship between adjusted R2 and estimation bias in our illustration, where the higher the adjusted R2 the more accurate the estimation of total load is. Thus, the adjusted R2, as the traditional model selection criteria in statistics or the Nash–Sutcliffe model efficiency in hydrology, is useful to choose appropriate models. We have presented results based on fortnightly and monthly sampling frequencies (with more samples in wet season). The Monte Carlo results using other sampling frequencies can also be easily obtained. For example, we also investigated weekly sampling schemes and have obtained similar conclusions. This is as one might expect because all the estimators here are regressionbased. Based on our studies, we can conclude that the new proposed procedures can significantly improve the predictability especially for large rivers. How much improvement by using the ADF will be site-dependent and the optimal discount factor should be site-specific reflecting the local environmental conditions. It would be exciting to see how the procedures work under different environmental conditions. Another fundamental challenge encountered in estimating sediment total load is the proper uncertainty analysis, such as model
78
Y.-G. Wang, T. Tian / Journal of Hydrology 482 (2013) 69–78
error, of the predictive model of concentration. To effectively investigate uncertainty in loads over a year, as described by Rustomji and Wilkinson (2008), using the bootstrap technique to resample the residuals produces reliable estimates of total sediment loads. Also, to compare the accuracies of sediment total loads, employing adjusted R2 as a reflection of regression model’s variation (i.e., how well a model describes the data) is a basis for developing the predictive model of concentration. On the other hand, improved surrogate technologies, such as optics and acoustics, are available to collect highly resolved time series of sediment concentration which will lead to more accurate estimates of annual load. Improved methods incorporating spatial–temporal dynamics will lead to better governmental decisions in managing costal and ecosystem health (Brodie et al., 2012). Acknowledgements We wish to thank two referees for their very constructive comments and detailed suggestions which led to a much improved paper. References Brodie, J.E., Kroon, F.J., Schaffelke, B., Wolanski, E.C., Lewis, S.E., Devlin, M.J., Bohnet, I.C., Bainbridge, Z.T., Waterhouse, J., Davis, A.M., 2012. Terrestrial pollutant runoff to the Great Barrier Reef: an update of issues, priorities and management responses. Mar. Pollut. Bull. 65, 81–100. Brown, H.E., Hansen, E.A., Champagne Jr., N.E., 1970. A system for measuring total sediment yield from small watersheds. Water Resour. Res. 6 (3), 818–826. Cohn, T.A., 1995. Recent advances in statistical methods for the estimation of sediment and nutrient transport in rivers. Rev. Geophys. 33, 1117–1123. Cohn, T.A., Caulder, D.L., Gilroy, E.J., Zynjuk, L.D., Summers, R.M., 1992. The validity of a simple statistical model for estimating fluvial constituent loads: an empirical study involving nutrient loads entering Chesapeake bay. Water Resour. Res. 28 (9), 2353–2363. Colby, B.R., 1956. The Relationship of Sediment Discharge to Streamflow, Rep., U.S. Geol. Surv., Reston, Va, p. 170. Darnell, R., Henderson, B., Kroon, F.J., Kuhnert, P., 2012. Statistical power of detecting trends in total suspended sediment loads to the Great Barrier Reef. Mar. Pollut. Bull. 65, 203–209. Davis, A.M., Lewis, S.E., Bainbridge, Z.T., Glendenning, L., Turner, R.D., Brodie, J.E., 2012. Dynamics of herbicide transport and partitioning under event flow
conditions in the lower Burdekin region, Australia. Mar. Pollut. Bull. 65, 182– 193. Ferguson, R.I., 1986. River loads underestimated by rating curves. Water Resour. Res. 22 (1), 74–76. Koch, R.W., Smillie, G.M., 1986. Bias in hydrologic prediction using log-transformed regression models. Water Resour. Bull. 22 (5), 717–723. Kuhnert, P.M., Henderson, B.L., Lewis, S.E., Bainbridge, Z.T., Wilkinson, S.N., Brodie, J.E., 2012. Quantifying total suspended sediment export from the Burdekin River catchment using the loads regression estimator tool. Water Resour. Res. 48, W04533. Kulasova, A., Smith, P.J., Beven, K.J., Blazkova, S.D., 2012. A method of computing uncertain nitrogen and phosphorus loads in a small stream from an agricultural catchment using continuous monitoring data. J. Hydrol. 458–459, 1–8. Miller, C.R., 1951. Analysis of Flow-duration, Sediment-Rating Curve Method of Computing Sediment Yield, Rep., U.S. Bur. of Reclam., Denver, Colo, p. 15. Rustomji, P., Wilkinson, S., 2008. Applying bootstrap resampling to quantify uncertainty in fluvial suspended sediment loads estimated using rating curve. Water Resour. Res. 44, W09434. http://dx.doi.org/10.1029/2007WR006088. Singer, M.B., Dunne, T., 2001. Identifying eroding and depositional reaches of valley by analysis of suspended sediment Transport in the Sacramento River, California. Water Resour. Res. 37 (12), 3371–3381 (Identifying eroding and depositional reaches of valley by analysis of suspended sediment Transport in the Sacramento River, California). Thomas, R.B., 1985. Estimating total suspended sediment yield with probability sampling. Water Resour. Res. 21 (9), 1381–1388. Thomas, R.B., 1988. Monitoring baseline suspended sediment in forested basins: the effects of sampling of suspended sediment rating curves. Hydrol. Sci. 33 (5), 499–514. Thomas, R.B., Lewis, J., 1993. A comparison of selection at list time and timestratified sampling for estimating suspended sediment loads. Water Resour. Res. 29, 1247–1256. Verhoff, F.H., Yaksich, S.M., Melfi, D.A., 1980. River nutrient and chemical transport estimation. J. Environ. Eng. Div. Am. Soc. Civ. Eng. 106, 591–608. Verma, S., Markus, M., Cooke, R.A., 2012. Development of error correction techniques for nitrate-N load estimation. J. Hydrol. 432, 12–25. Wallace, J., Stewart, L., Hawdon, A., Keen. R., 2008. The role of coastal floodplains in generating sediment and nutrient fluxes to the great barrier reef lagoon in Australia. In: Ecohydrological Processes and Sustainable Floodplain Management Opportunities and Concepts for Water Hazard Mitigation and Ecological and Socioeconomic Sustainability, Lodz, Poland. Walling, D.E., 1977. Assessing the accuracy of suspended sediment rating curves for a small basin. Water Resour. Res. 13 (3), 531–538. Wang, P., Linker, L.C., 2008. Improvement of regression simulation in fluvial sediment loads. J. Hydraul. Eng. 134 (10), 1527–1531. Wang, Y.-G., Kuhnert, P., Henderson, B., 2011. Load estimation with uncertainties from opportunistic sampling data – a semiparametric approach. J. Hydrol. 396, 148–157.