9 The authors explore simple concepts of persistence in streamflow forecasting based on the 10 real-time streamflow observations from the years 2002 to 2018 at 140 U.S. Geological Survey 11 (USGS) streamflow gauges in Iowa. The spatial scale of the basins ranges from about 7 km2 to 12 37,000 km2. Motivated by the need for evaluating the skill of real-time streamflow forecasting 13 systems, the authors perform quantitative skill assessment of different persistence schemes 14 across spatial scales and lead-times. They show that skill in temporal persistence forecasting has 15 a strong dependence on basin size, and a weaker, but non-negligible, dependence on geometric 16 properties of the river networks in the basins. Building on results from this temporal persistence, 17 they extend the streamflow persistence forecasting to space through flow-connected river 18 networks. The approach simply assumes that streamflow at a station in space will persist to 19 another station which is flow-connected; these are referred to as pure spatial persistence 20 forecasts (PSPF). The authors show that skill of PSPF of streamflow is strongly dependent on 21 the monitored vs. predicted basin area-ratio and lead-times, and weakly related to the 22 downstream flow distance between stations. River network topology shows some effect on the 23 hydrograph timing and timing of the peaks, depending on the stream gauge configuration. The 24 study shows that the skill depicted in terms of Kling-Gupta efficiency (KGE) > 0.5 can be 25 achieved for basin area ratio > 0.6 and lead-time up to three days. The authors discuss the 26 implications of their findings for assessment and improvements of rainfall-runoff models, data 27 assimilation schemes, and stream gauging network design. 28 1 Introduction 29 In this paper we demonstrate, in a systematic way, that persistence is a hard-to-beat 30 streamflow forecasting method (Palash et al. 2018), a fact well-known to operational forecasters. 31 We limit our considerations to real-time forecasting in space and time in a river network. River 32 networks aggregate water flow that originates from the transformation of rainfall and/or 33 snowmelt to runoff. Additional recognition of the importance of river networks stems from the 34 fact that many communities reside by streams and rivers and thus require streamflow forecasts 35 on a regular basis for a number of applications ranging from flood protection to water resource 36 management. 37 Our study is motivated by ten years of experience of the Iowa Flood Center (IFC) providing 38 operational real-time streamflow forecasts for the citizens of Iowa (Demir and Krajewski 2013; 39 Krajewski et al. 2017). The IFC has accomplished this in a traditional way, i.e. using a state-40 wide implementation of a distributed rainfall-runoff model. The IFC model is referred to as the 41 Hillslope-Link Model (HLM) and has its origins in landscape decomposition discussed by 42 Mantilla and Gupta (2005). The HLM has been presented and its performance discussed in 43 several prior studies, including Cunha et al. (2012), Ayalew et al. (2014), ElSaadani et al. (2018), 44 and Quintero et al. (2018). In addition, the National Weather Service, via its North Central and 45 Central River Forecasts Centers issues “official” streamflow forecasts for a number of locations 46 throughout the state of Iowa (e.g. Zalenski et al. 2017). 47 In this study, we take a step back and ask the question: how well could someone forecast 48 streamflow without a rainfall-runoff model? What would be the skill and limitations of such 49 forecasts? We address these questions using the approach of simple persistence applied both in 50 space and in time. We explore a range of forecast lead times, from hours to days ahead, and 51 characterize their skill according to several popular metrics used in evaluation of hydrologic 52 models. Our study domain is the entire state of Iowa but without considering the border rivers, 53 i.e. the Mississippi and the Missouri Rivers. Iowa has a dense drainage network of streams and 54 rivers that feed its border rivers. The state is home to some 1,000 communities located near 55 streams and rivers and is frequented by severe storms and the resultant floods. The U.S. 56 Geological Survey (USGS) operates and collects streamflow data in real time at some 140 stream 57 gauges in Iowa and the IFC has deployed 250 stage-only sensors (Kruger et al. 2016; Krajewski 58 et al. 2017). 59 Streamflow forecasting in space without using hydrologic models is little studied, and 60 limited efforts have been documented in the literature. Some recent probabilistic approaches 61 (Botter et al. 2008; Müller et al. 2014) include a physically based mechanism for hydrologic 62 predictions in sparsely gauged basins under changing climatic conditions (Doulatyari et al. 2015; 63 Müller and Thompson 2016; Dralle et al. 2016). However, the spatial cross-correlations among 64 multiple basin outlets have not been studied in this framework (Betterle et al. 2017). 65 Specifically, geo-statistical methods have been effectively employed to explore spatial patterns 66 in streamflows accounting for the river network topology (e.g. Betterle et al., 2017; Castiglioni et 67 al., 2011; Laaha et al., 2014; Skøien et al., 2006). Spatial streamflow correlation has been 68 proven to be a better indicator in transferring the flow attributes to ungauged outlets from gauged 69 outlets in comparison to proximity (Betterle et al. 2017). Notwithstanding the progress achieved 70 over the last few years in spatial streamflow predictions, further efforts on exploiting the spatial 71 and temporal persistence of streamflow in the context of streamflow forecasting is needed and 72 this study follows in that direction. 73 We organize this paper as follows. In the next section, Experimental Domain and Data, we 74 present the data used in our investigations and describe Iowa and its major river basins. In the 75 Methods section, we describe the experimental setup for our analysis and discuss the 76 performance metrics used to characterize the skill of the persistence forecasting. The Results 77 section documents our findings reported as a function of spatial scale characteristics. We also 78 discuss the connection to the river network geometry depicted in terms of width function. We 79 end with a discussion of implications of our study for operational forecasting systems, 80 hydrologic and hydraulic models, river basin gauging network design as well as data assimilation 81 and Artificial Intelligence (AI) methods. 82 2 Experimental Domain and Data 83 We explore the forecasting skill of simple persistence in both time and space. Our domain 84 of interest is the state of Iowa and its interior rivers (Figure 1). Iowa is bordered by the 85 Mississippi River on the east and the Missouri River on the west side. The interior rivers all 86 discharge to these two large rivers with about two-thirds of the state’s 150,000 km2 area draining 87 to the Mississippi. Most drainage basins of those interior rivers are contained within the 88 administrative borders of the state, but a small portion drains parts of southern Minnesota and 89 South Dakota. 90 The climate of Iowa is characterized by marked seasonal variations due to its latitude and 91 interior continental location. During the warm season, the moist air coming from the Gulf of 92 Mexico produces a summer rainfall maximum; during the cold season, flow of dry Canadian air 93 causes cold and dry winters. The air masses moving across the western United States throughout 94 the year produce mild and dry weather. 95 The landscape is predominantly agricultural, characterized by low relief, fertile soils, high 96 drainage density and thus numerous rivers (e.g., Prior 1991). Much of the state is covered by tile 97 drains, thus changing the hydrologic regime of the region (e.g., Schilling et al. 2019). The 98 combination of the characteristics of its landscape, together with significant rainfall during the 99 summer months makes the state prone to the occurrence of floods (e.g., Villarini et al. 2011a, 100 2013; Mallakpour and Villarini 2015; Villarini et al. 2011b). 101 102 Figure 1: USGS stream gauges over the state of Iowa used in this study. Numbered labels show 103 the river basins with multiple gauges thus enabling assessment of skill of persistence forecasting 104 in space. 105 The U.S. Geological Survey currently operates about 140 stream gauges in the state (Figure 106 1), which provide 15-minute resolution stage data that can be converted to discharge using well-107 maintained and frequently updated rating curves. The USGS record of 15-minute data has only 108 17 years but we claim that this is sufficient to illustrate the objectives of our study. Though the 109 Iowa Flood Center has added 250 stage-only gauges (Kruger et al. 2016), we will not use the 110 data in this study. Figure 2 shows the basin size histogram at the outlet of these two sets of 111 gauges. It is clear that the IFC emphasizes monitoring smaller basins than the USGS. 112 113 Figure 2: Distribution of the IFC (blue) and the USGS (green) stream gauges over the state of 114 Iowa. It is apparent that most USGS gauges monitor streamflow at larger basin scales. 115 In addition to the discharge data, we will make use of the water flow velocity data collected 116 by the USGS. The velocity data are synthetized into a simple power-law model documented by 117 Ghimire et al. (2018). In the later part of the paper, we will refer to the model simulation results 118 for the Iowa domain. The model, called Hillslope Link Model (HLM) is briefly presented in 119 Krajewski et al. (2017) and several other publications (e.g., Quintero et al. 2016; ElSaadani et al. 120 2018; Cunha et al. 2012; Ayalew et al. 2014b; Small et al. 2013). A notable point of the model’s 121 use in the operation of the IFC is that the model has not been calibrated in the classic meaning of 122 the practice. 123 3 Methods 124 3.1 Experimental Design 125 We investigated the streamflow forecasting performance of a simple persistence approach in 126 time and in space. We considered a range of spatial scales (~10-10,000 km2) as well as a range 127 of the forecast lead times (15 minutes to 5 days). Let us denote locations where we would like to 128 forecast streamflow with index i and the stations which data we use to make the forecast as j. In 129 the time domain, the streamflow Q forecasting equation becomes simply 130 𝑄௜(𝑖,𝑡଴+Δ𝑡)=𝑄௜(𝑖,𝑡଴) (1) 131 where, the time of issuing the forecast is 𝑡଴ and Δ𝑡 is the forecast lead time, Δ𝑡=[0, Δ𝑡max]. Here, 132 we simply repeat the main findings documented in Ghimire and Krajewski (2019). These 133 authors investigated three different persistence approaches for Δ𝑡 ranging between 15 minutes 134 and five days, but here we present the results of the simplest one only, i.e. that the future will be 135 the same as the present, as given by (1). In Figure 3, we illustrate the result from this approach 136 with 15-minutes of streamflow data from 2002 to 2018 at about 140 USGS stations. It is clear 137 that for small scale basins, where the streamflow variability is large and the basin response time 138 short, the persistence forecasting quickly loses skill. For large basins, the skill persists for 139 considerable forecast lead time. 140 141 Figure 3: Variability of Kling-Gupta Efficiency based on simple persistence forecasting 142 approach across basin scales for different lead-times (after Ghimire and Krajewski 2019). The 143 forecast skill in terms of KGE shows clear dependence on basin scales and forecast time horizon. 144 In space, to be fully parallel to the investigation in the time domain, we first investigated 145 scenarios where time is not involved. Here, however, the 140 USGS stream gauges must be 146 stratified not only by basin size but also by basin membership. We only considered gauges that 147 belong to the same basin and thus are flow-connected in either downstream or upstream direction 148 (but not both). Therefore, we excluded consideration of hydroclimatic similarity and land use or 149 ecosystem classification (e.g., Schilling et al. 2015). For example, a basin in western Iowa is 150 excluded from considerations for persistence forecasting in eastern Iowa. 151 In Figure 1, we show the major basins in the state and in Table 1 we provide details on 152 USGS gauges in each. Now, within a basin we can consider all potential forecasting vs. 153 evaluation scenarios. For example, consider the situation of a supposed basin with five gauges 154 as numbered in Figure 4. First, consider downstream direction only. Hypothetically, one could 155 apply the persistence forecasting paradigm using data from Gauge 1 to forecast streamflow at 156 Gauge 5. The results might not be very good as Gauge 1 misses the contributions from a large 157 portion of the basin and major tributaries upstream from Gauge 5. Similarly, one could evaluate 158 the forecasting performance at Gauge 5 using data from Gauges 2, 3, and 4. Or, one could 159 evaluate the spatial persistence at Gauge 3 using data from Gauges 1 or 2 but not 4 or 5. Gauge 160 5 is flow-connected with Gauge 3 but in the upstream direction. Gauge 4 is not flow-connected 161 with Gauge 3 as the connection would require downstream and upstream direction for water 162 travel. 163 Table 1: Basin characteristics with linear regression-based model parameters for KGE 164 165 ID Basin Area [km2] Number of Stations Number of Stations/km2 α β R2 RMSE 1 Wapsipinicon River 6,000 4 0.67 -0.45 1.44 0.94 0.11 2 Iowa River 7,300 11 1.51 -0.39 1.41 0.97 0.09 3 Skunk River 11,300 11 0.97 -0.43 1.44 0.97 0.08 4 Cedar River 20,100 20 1.00 -0.44 1.43 0.98 0.07 5 Des Moines River 14,300 9 0.63 -0.61 1.60 0.96 0.10 6 Raccoon River 9,400 9 0.96 -0.41 1.39 0.96 0.09 7 Turkey River 4,100 5 1.22 -0.38 1.39 0.97 0.08 8 Little Sioux River 6,300 4 0.63 -0.41 1.40 0.98 0.08 * 𝐾𝐺𝐸=𝛼+𝛽𝐴𝑟 166 Figure 4: Demonstration of experimental setup for the spatial persistence forecasting. Arrows 167 show the direction for streamflow forecast in a flow connected river network. The schematic 168 illustrates Single Point Forecasts and Multiple Point Forecasts. In the upstream direction, only 169 SPF is possible. 170 Similarly, if we consider the upstream direction, one could make forecast at Gauge 1 using 171 data from Gauge 5. Full analogy follows for all other pairing. In summary, for this example 172 basin, we have seven pairs to consider in each direction. These pairs of the forecast and 173 evaluation points will allow us to evaluate persistence forecasting in space as function of spatial 174 scale. We will refer to these forecasting scenarios as single point forecasts (SPF). 175 The above discussion is valid for the forecast lead time equal to zero, i.e. streamflow 176 observed now at location i is the same as the forecast at location j. We can extend the same 177 construction to other lead times Δ𝑡= 0,…,Δ 𝑡௠௔௫. The forecast at location i for the lead time Δ𝑡 178 is the same as the flow now at the stream gauge j. In terms of the spatial arrangement between 179 the forecast and evaluation points, nothing changes i.e. 180 𝑄௜(𝑖,𝑡଴+Δ𝑡)=𝑄௝(𝑗,𝑡଴) for Δ 𝑡= 0,…,Δ 𝑡௠௔௫ (2) 181 Now, let us expand the forecasting setup. Following the same example in Figure 4, it is 182 clear that if one adds data from Gauge 3 and Gauge 4, the forecasts at Gauge 5 should improve, 183 perhaps considerably. However, forecasting at Gauge 3 does not benefit from considering data 184 from Gauge 1 and Gauge 2 together. For the example basin in Figure 4, there is only one multi 185 point forecast (MPF) combination possible. However, for larger basins with rich river networks 186 and dense system of gauges, there could be several to consider. Note, that in the upstream 187 direction there are no MPFs possible. As earlier, considering different forecast lead times does 188 not change the spatial arrangements between the forecast and evaluation points. Mathematically, 189 𝑄௜(𝑖,𝑡଴+Δ𝑡)=∑𝑄௝(𝑗,𝑡଴)௡ ௝ୀଵ for Δ 𝑡= 0,…,Δ 𝑡௠௔௫ (3) 190 where n is the number of stream gauges that constitute an MPF. 191 Next, let us add a simple travel time consideration. Since it is well established that in most 192 places the direction of water movement is easy to predict, we can add water flow velocity and 193 river distance information to the persistence forecasts. The simplest choice for water flow 194 velocity is a constant value. Mathematically, for a SPF 195 𝑄௜(𝑖,𝑡଴+Δ𝑡)=𝑄௝ቀ𝑗,𝑡଴+ௗ ௩ቁ for Δ 𝑡= 0,…,Δ 𝑡௠௔௫ (4) 196 where v is a constant water flow velocity and d is the distance along the river network between 197 locations i and j. Recognizing that water flow velocity is a non-linear function of discharge and 198 other channel parameters (e.g., Ghimire et al. 2018; Mantilla 2007; Ayalew et al. 2014a,b), we 199 limit our consideration to the constant velocity case to keep the setup simple. 200 As river networks display a wide range of shapes (topologies) that may affect hydrologic 201 response (e.g., Ayalew and Krajewski 2017; Perez et al. 2018), it is interesting to check their 202 effect on the forecasting skill of persistence (also see Ghimire and Krajewski 2019). To 203 accomplish this, we consider the width function of the basins defined by the existing USGS 204 gauges. Basin width function can be interpreted as a distribution of stream segments (distances) 205 with respect to the outlet of the basin (e.g., Rodriguez-Iturbe and Rinaldo 1997). If one assumes 206 constant water flow velocity, the distribution of distances is equivalent to the distribution of 207 travel times. Consider the width functions of the Cedar River at Cedar Rapids and the Iowa 208 River at Marengo (Figure 5). We can see that the maximum of the travel time for the Cedar 209 River basin is located some 200 km away from the outlet at Cedar Rapids. For the Iowa River, 210 the maximum is just upstream from Marengo, some 50 km from the outlet. Note that in Figure 211 5, the travel time between the forecast point and the evaluation point (basin outlet) is the same. 212 In addition, the base of the width function is also similar (~350 km). However, the forecast point 213 of the Cedar basin entails not only the peak of the width function, but also a larger upstream 214 drainage area. 215 These differences in the shape of the basin may influence the collective performance of 216 persistence-based forecasting, especially for the MPF scenarios. However, since considering the 217 entire distribution is inconvenient, we examined several simple geometric descriptors of width 218 functions (e.g., Perez et al. 2018) such as base of the width function( DB), centroid of the width 219 function ( 𝑑̅), mass of the width function ( M), and distance to the peak of the width function 220 (dWmax). Perez et al. (2018) and Ghimire and Krajewski (2019) pointed out that normalized form 221 of DB explains the most variability in performance metrics. Whether the forecasting scenario is 222 SPF or MPF, we define relationships between forecast and evaluation points in pairs. The ratio 223 of DB explains the same variability as normalized DB, which we discuss later in detail. 224 225 Figure 5: Width function depicting normalized histogram (density) of distance of each link in the 226 river network to the basin outlet. (a) Width function for the Cedar River at Cedar Rapids and (b) 227 Width function for the Iowa River at Marengo. The red dot represents the Single Forecast Point 228 with the evaluation site being the basin outlet. 229 3.2 Evaluation Metrics 230 We evaluated persistence derived streamflow forecasts in space and time against USGS 231 streamflow observations at evaluation points. We computed forecast skills in terms of several 232 statistical metrics. The major performance skill measures we report here are Kling-Gupta 233 Efficiency ( KGE), non-parametric KGE, mean absolute error ( MAE), normalized MAE (nMAE), 234 hydrograph timing ( TH), timing of the annual peak ( TP), and percent peak difference ( PPD). The 235 KGE metric has three components: Pearson’s correlation ( r), variance ratio ( ) and mean ratio 236 (β) (see Gupta et al. 2009). The ideal value of KGE is equal to 1. 237 𝐾𝐺𝐸= 1−ඥ(𝑟−1)𝟐+(𝛼−1)𝟐+(𝛽−1)𝟐 (5) 238 where 𝛼=ఙೞ ఙ೚, and 𝛽=ఓೞ ఓ೚, σs is the variance of forecasts, σo is the variance of observations, μs is 239 the mean of forecasts, and μo is the mean of observations. As an example, if the variance ratio  240 and the mean ration  are close to one, the KGE is dominated by the correlation coefficient r. 241 The frontier of KGE for good/bad forecasts relative to the mean of the observations could be 242 considered as -0.41 (Knoben et al. 2019). A modification was proposed by Pool et al. (2018) to 243 get away with the implicit assumption of linearity of data, normality of data and the absence of 244 outliers in the formulation of KGE by Gupta et al. (2009). It comprises non-parametric 245 components, i.e. Spearman rank correlation and normalized flow duration curve while the mean 246 ration (multiplicative bias) component stays the same. The modified measure is given by 247 𝐾𝐺𝐸ே௉= 1−ඥ(𝑟௦−1)𝟐+(𝛼ே௉−1)𝟐+(𝛽−1)𝟐 (6) 248 where rs is the Spearman rank correlation and αNP is the normalized flow duration curve. For 249 further details, readers are referred to Pool et al. (2018). 250 Likewise, we report MAE for stage forecasts associated with streamflow forecasts derived 251 using rating curves and is given by, 252 𝑀𝐴𝐸=ଵ ே|ℎ𝒔−ℎ௢| (7) 253 where hs and ho are forecast and observed stage, respectively, and N is the total number of time- 254 series nodes. Also, it is interesting for streamflow forecasters to consider stage forecast error in 255 terms of actual values depicted in terms of median (Δ hmedian) i.e. 256 ∆ℎ௠௘ௗ௜௔௡=𝑚𝑒𝑑𝑖𝑎𝑛(ℎ௦−ℎ𝒐) (8) 257 The normalized MAE (nMAE) for streamflow forecasts is computed similar to (7) but 258 normalized by upstream drainage area. Hydrograph timing ( TH) is associated with the cross- 259 correlation between the entire observed and forecast time-series. We computed it as the number 260 of hours a forecast time-series needs to be shifted so that the cross-correlation coefficient is 261 maximized. A positive value indicates delay in the timing of hydrograph while a negative value 262 indicates early timing of forecast relative to the observed hydrograph. Timing of the peak TP is 263 simply the difference in time at which the peak of forecast occurs relative to the peak of 264 observations. 265 Next, we computed PPD as, 266 𝑃𝑃𝐷=௣௘௔௞𝒔𝒊𝒎ି௣௘௔௞ 𝒐𝒃𝒔 ௣௘௔௞𝒐𝒃𝒔×100 (9) 267 where, peaksim and peakobs are the peaks of forecasts and observation respectively. We limit the 268 consideration only to annual peaks (based on the observations). 269 4 Results 270 4.1 Purely Temporal Persistence Forecasts 271 While Ghimire and Krajewski (2019) provide a comprehensive characterization of three 272 persistence schemes, here we present their key result as it will be useful for our discussion of the 273 spatio-temporal persistence. In Figure 3, we show the KGE statistic for different forecast lead 274 times up to three days ahead organized as a function of basin drainage area. It is interesting to 275 note that for one-day ahead forecast for basins of the size of about 1000 km2, this simple 276 persistence method skill is about 0.8, equivalent to a very good model (for a discussion of what 277 is very good, see the analysis in Knoben et al. (2019). As Ghimire and Krajewski (2019) show, 278 this skill can be even somewhat increased by incorporating simple climatology information. 279 4.2 Purely Spatial Persistence Forecasts (PSPF): Downstream Direction 280 In this scenario, the point of interest, e.g. a community, is downstream from the forecast 281 point. In this study, we assume we have a stream gauge at the forecast location and can use the 282 data to evaluate the quality of the forecast that used data from a point upstream. For example, 283 consider the Cedar River basin and take streamflow data at Waverly, Iowa (Gauge ‘2’ in Figure 284 6) and use them as the forecast for Cedar Rapids (Gauge ‘5’ in Figure 6). The performance is 285 not very good: the KGE is -0.1 with components of correlation of 0.49, the mean ratio of 0.25, 286 and the variance ratio of 0.38. However, if we now use data at Waterloo, Iowa (Gauge ‘3’ in 287 Figure 6) and use them to make forecasts for Cedar Rapids, the skill increases dramatically: the 288 KGE goes up to 0.7. Now, we summarize the performance results for all possible combinations 289 of the Single Point Forecasts (SPFs) in Figure 6 (see red dots) organizing them as function of the 290 upstream drainage area ratio, hereafter referred to as monitored area fraction, Ar (e.g. gauge at 291 Waterloo monitors drainage of about 0.8 of the basin area at Cedar Rapids). It is evident that 292 using data from stations that capture over 70% of the basin of interest, i.e. Ar = 0.7, persistence- 293 based forecasting performs very well (KGE > 0.5). The results shown here for the Cedar River 294 basin, indicate strong relationship of the KGE measure and the monitored area fraction. 295 296 Figure 6: Relationship of the KGE metric with upstream drainage area ratio, Ar for the Cedar 297 River basin for the year 2016. Every dot in this figure implicitly represents a pair of forecast and 298 evaluation points. Red dots show SPF performance while blue dots represent MPF forecast 299 performance. As a demonstration of improvement in performance of MPF based forecast from 300 SPF based forecast consider the green and magenta points. The green dots are evaluation points 301 3 to 7 using SPF at forecast point 1. The magenta dots are the same evaluation points upon 302 addition of forecast point 2 to 1 showing a dramatic increase in the KGE metric. 303 The performance increases if we combine data coming from different upstream branches of 304 the river network (MPFs). A simple demonstration of how MPFs contribute to performance 305 increase is shown by blue dots in Figure 6. Here, green dots with forecast issued from Gauge 1 306 transforms to magenta dots upon addition of forecast from Gauge 2. The benefit does not apply 307 to all stations as there are many locations that have only one station upstream. However, it is 308 interesting to note that the performance based on MPFs also follows the strong relationship with 309 Ar. 310 The above results are based on data for 2016. A natural question arises whether the same 311 performance stays in other years. To illustrate the year-to-year variability, we have explored 312 data for the 17-year period, 2002-2018. In Figure 7, we show violin plots depicting distribution 313 of KGE over this period for eight major basins (refer Figure 1 and Table 1) in Iowa. It is clear 314 that both the median KGE (red line) and overall distribution display the same signature as in 315 Figure 6. A linear-regression based predictive model with goodness-of-fit measures are 316 presented in Table 1. Clearly, R2 and RMSE measures corroborates the strong relationship of 317 KGE with Ar. Note that the intercept, α and slope parameter, β of the predictive models show 318 very small variability across basins except the Des Moines River. We speculate that this 319 difference could be attributed to the distinct hydrogeology of this basin much of which is in the 320 Des Moines Lobe characterized with pothole connectivity augmented by an extensive tile 321 drainage network. 322 323 Figure 7: Violin plots showing inter-annual variability of KGE metric with upstream drainage 324 across basins in Iowa. Each violin plot shows distribution of KGE over the period of 2002–325 2018. The red line and black line show median and mean, respectively, of the KGE over the 326 years. 327 In Figure 8, we show relationship of KGE and its decomposition components with Ar for the 328 pooled data from all basins for the year 2016. As discussed above, KGE statistic is explained 329 strongly by Ar. Among its components, both bias and variance ratio show similar signature with 330 Ar while correlation does not show as strong relationship. It is interesting to explore if river 331 network geometry in terms of width function can explain the variability in KGE and its 332 components. It is clear from Figure 9 that a width function descriptor of DB ratio does not 333 explain the variability of KGE and its components as strongly as Ar. However, note that 334 correlation is explained almost equally well by the DB ratio. If this result is interpreted in 335 conjunction with Figure 5, larger variability in correlation with Ar can also be explained by DB 336 ratio. This, in turn, has particular implication for TH and TP, which we discuss later. 337 338 Figure 8: Relationships of KGE metric decomposition components with upstream drainage area 339 ratio for pooled data from all basins in Iowa for the year 2016. KGE (a) shows strong linear 340 dependence with Ar. The mean ratio (b), variance ratio (c) and correlation coefficient (d) show 341 decreasing degree of dependence on Ar. 342 343 Figure 9: Relationships of KGE metric decomposition components with base of the width 344 function, DB ratio, for pooled data from all basins in Iowa for the year 2016. The dependence of 345 KGE and its components with DB ratio is not as strong as with Ar (see Figure 8). Note that, the 346 correlation is explained similarly by both DB ratio and Ar. 347 While KGE may not be the ultimate performance criterion for the operational forecasters, its 348 regular behavior with the monitored area fraction provides high level of confidence in the 349 streamflow forecasting. However, for many applications, including flood forecasting, the 350 variable of most interest is river stage. In Figure 10, we evaluate our persistence forecasting in 351 terms of stage errors. We use the rating curves provided by the USGS to convert the propagated 352 discharge to stage. As the relationship between discharge and stage is influenced by the local 353 hydraulic properties of the river channel, and because stage is measured with respect to an 354 arbitrary datum, we expressed the results in terms of stage error in units of length [m]. In Figure 355 10, we show the distribution of median of the stage errors obtained over a period of seventeen 356 years across eight basins in Iowa. The year-to-year variability is high, in particular, for cases of 357 small Ar. Somewhat surprisingly, the median stage forecast errors are small for Ar larger than 358 0.7. For Ar > 0.7, median stage forecast errors are within 30 cm (one foot) for all basins 359 indicating the potential of PSPF, which model-based forecasts should improve upon. 360 361 Figure 10: Violin plots showing inter-annual variability of median stage forecast error with 362 monitored area fraction, Ar across basins in Iowa. Each violin plot shows distribution of stage 363 forecast error over the period of 2002–2018. The red line shows the median of stage errors over 364 the years. 365 In addition to the skill measures discussed above, operational forecasters are interested in 366 other statistics while making streamflow forecasting decisions. In Figure 11, we show variability 367 of nMAE, PPD, TH and TP over the seventeen-year period with Ar for the Cedar River basin. The 368 variability and median values of nMAE and PPD show similar signature as KGE with Ar. These 369 relationships of TH and TP with Ar are similar but weaker. This weaker relationship is partly 370 explained by the variability of correlation explained by the DB ratio (also see Figure 9). 371 372 Figure 11: Box plots showing inter-annual variability of forecast skill measures normalized MAE 373 (nMAE), percent peak difference ( PPD), hydrograph timing ( TH) and timing of the peaks ( TP) 374 with monitored area fraction, Ar for the Cedar River basin. The nMAE and PPD show strong 375 dependence on Ar while the signature of TH and TP with Ar is not as strong. 376 4.3 Purely Spatial Persistence Forecasts (PSPF): Upstream Direction 377 In this scenario, the point of interest, e.g. a community or a water intake or a waste water 378 discharge, is upstream from the forecast point. The results are similar but, surprisingly, worse 379 than in the downstream direction, at least in terms of the KGE statistic (Figure 12). This lack of 380 symmetry between the downstream and upstream direction can be easily explained by the 381 construction of the KGE metric. While the correlation coefficient between the predicted and 382 observed series is immune to the upstream/downstream translation, the other terms are ratios. In 383 the downstream direction the predicted values are systematically lower than the observed (by 384 construction of the forecasting scheme), the mean ratio is always less than one. For the upstream 385 direction, the ratio has the inverse value. The variance ratio behaves similarly. As a result, the 386 KGE values are generally lower for the upstream case and not better than the mean of the 387 observations, especially for the Ar below 0.5. Note that we kept the ratio constraint by the larger 388 scale gauge, although it is inconsistent with its original definition (forecasted vs. monitored 389 area). 390 391 Figure 12: Effect of forecast direction in the river network on the relationship between KGE and 392 Ar. With upstream direction of persistence forecast, the persistence forecasts lose skill faster 393 with Ar. 394 4.4 Spatio-temporal Persistence Forecasts 395 Here, we consider forecasting scenarios with lead-time, Δ t > 0 to the purely spatial 396 persistence forecasts (PSPF) (see equation 2 and 3). In the same framework discussed in 4.3, we 397 investigate the two-dimensional variability of performance metrics, i.e. with monitored area 398 fraction, Ar and the lead-time, Δ t. In Figure 13, we show contour plots depicting the spatio- 399 temporal persistence forecast skill in terms of non-parametric KGE for the Cedar River basin. 400 Since these contour plots are empirically produced based on a limited sample size, lack of 401 smoothness of the contours is expected. As a demonstration, we show here three basin scales in 402 top to bottom panels (Figures 13a-13c) which represent the SPF points used to issue forecast at 403 locations downstream. It is interesting to note that, for all scales at Δ t = 3 days and Ar = 0.6, 404 KGENP is larger than 0.5. Following from our discussion in 4.3, one can fit regression models 405 with Ar and Δt as explanatory variables giving rise to a contour plot that works for all basin 406 scales and can be used by forecasters for model evaluation. We have developed similar results 407 (not shown) in terms of the KGE statistic. 408 409 Figure 13: Contour plots depicting the spatio-temporal persistence forecast skill in terms of non-410 parametric KGE (see contour labels) for the Cedar River basin. Three basin scales from top to 411 bottom panels (a-c) represent the SPF points used to issue forecast at locations downstream. For 412 example, for all scales at lead time of three days and Ar of 0.6, KGENP is larger than 0.5. 413 4.5 Adding Velocity Effect 414 Staying true to the principle of extreme simplicity of this investigation, we would like to 415 explore the effect of downstream water transport to improve the forecasting skill. Intuitively, it 416 seems reasonable to expect that since water flows downstream, we should be able to improve the 417 forecasting by applying the proper time offset. Since the distance along the river is known, the 418 other piece of information we need is the water flow velocity. Assuming a constant velocity, the 419 question becomes what value should we use? We explored the empirical distribution of 420 streamflow velocity based on USGS gauges in Iowa. We investigated over 100,000 data points 421 for the velocity distribution. The observed velocity in Iowa’s streams and rivers ranges from 0.2 422 to over 2 m/s. We investigated the velocities that correspond to each 10th quantile of the 423 distribution. It turns out that the best choice, in terms of maximizing the KGE, is a velocity 424 significantly higher than the median or the mean but not the very highest velocity. 425 In Figure 14, we show the dependence of the KGE and its components on the chosen 426 constant velocity (each 10th quantile of the distribution) for the Cedar River basin. We show the 427 components of the KGE because the metric itself shows little sensitivity with respect to velocity. 428 In the Cedar River basin, we have a large number of USGS gauges, some added following the 429 2008 extreme flood (Smith et al. 2013). This results in a large number of pairs of gauges and 430 thus a wide range of the Ar. In Figure 14, we show (in red dots) all possible SPF pairs. Those 431 with small values of Ar result in poor performance in terms of KGE. To facilitate the analysis, 432 we selected two pairs: one with a small value of Ar (~0.04) shown in green, and second one with 433 high (but not highest) Ar (~0.81) shown in blue. The corresponding color horizontal lines show 434 the KGE for the purely spatial persistence forecasting (PSPF). Recall that “purely spatial” refers 435 to the case of Δ𝑡 = 0. 436 437 Figure 14: Effect of constant velocity routing on Pure SPF forecasts skill for the Cedar River 438 basin. The panels show effect of constant velocity on each component of KGE. The velocity 439 axis shows the constant velocity used for routing PSPF forecasts while vertical axis shows KGE 440 components conditional on a given constant velocity. Each grey line connecting the red dots 441 represents a value of Ar. Green and blue horizontal lines show the corresponding PSPF forecast 442 skills for smaller ( Ar=0.04) and larger ( Ar=0.81) area ratios. The green and blue lines with dots 443 show forecast skills for respective area ratios with the addition of constant routing velocity. 444 Several aspects of the problem become clear from the plots. First, the sensitivity of the 445 KGE is dominated by the mean and the variance ratios. These two components show no 446 sensitivity to the velocity as the same hydrographs are passed between the two gauges, just some 447 faster and some slower. The sensitivity is in the correlation coefficient. Increasing the velocity 448 improves the correlation but only to a point. Further increase of the velocity will compromise 449 the correlation. Also, the sensitivity is higher for the small Ar. The gauge pairs characterized by 450 high Ar, show negligible impact of the velocity. This implies that the forecasting performance is 451 dominated by the Ar. However, for some values of Ar, including velocity from the mid-range of 452 the observed values, improves the forecasting performance, again, mostly in terms of the 453 correlation. 454 Deeming the results found for the spatio-temporal persistence forecasts in 4.4, we expect 455 that similar results will hold upon addition of constant velocity. The effect of constant velocity 456 on other statistics, such as nMAE, PPD, TH and TP, results directly from the above results. The 457 values of TH and TP depend mainly on the correlation. 458 4.6 Closing Comments 459 We hope that readers find the results of this simply constructed study interesting and useful. 460 We demonstrate that data-only based real-time streamflow forecasting has a high level of skill 461 for well-instrumented basins. Our results could be useful in several ways. We establish a 462 performance reference level that more elaborate approaches, such as distributed hydrologic or AI 463 based models should be expected to improve on. This reference is location-dependent based on 464 the configuration of the stream gauges. One simple demonstration of this is in Figure 15. Here, 465 we show the KGE along the river network in the Cedar River basin obtained by using the linear 466 relationship from Figure 6. The skill changes as the distance from the gauges increases (as the 467 monitored area fraction changes as well). For each segment of the river network (or a link in the 468 HLM nomenclature), we find the maximum skill for upstream or downstream (SPF or MPF) 469 options. Adding IFC sensors would expand the useful skill of the persistence (assuming 470 “perfect” rating curves). 471 472 Figure 15: Demonstration of Pure SPF skill in terms of KGE in the river network using the 473 USGS only (a), and USGS and IFC bridge sensors (b) as forecasting locations for each link. 474 Addition of IFC bridge sensor in the network shows the potential for enhanced coverage of 475 forecasts and improved the corresponding skill in the river network. The plot assumes perfect 476 rating curves for both USGS and IFC gauge location. 477 The above results could serve as a basis for streamflow gauging network design as they 478 clearly guide placing additional stream gauges to increase the Ar. They can also help with the 479 assessment of the quality of synthetic rating curves developed for stage-only gauges (Kruger et 480 al. 2016). 481 We foresee several extensions of this work in the direction of systematic assessment of the 482 forecasting skill of hydraulic and hydrologic models. For example, we have demonstrated that it 483 is difficult to find a constant velocity that improves upon a purely spatial approach. Would 484 hydraulic models of open channel flow lead to better results? They require additional 485 information in the form of channel geometric and hydraulic (roughness) properties. While the 486 cross sections and slopes are straightforward to obtain, parametric description of roughness is 487 subject to considerable uncertainty. 488 Finally, while comparing the PSPF skill to that of highly distributed (physically based) 489 hydrologic models is outside the scope of this paper, we include an example in Figure 16. It 490 shows the KGE of an open loop HLM model as well as the same model with a simple 491 streamflow data assimilation scheme. The scheme is a replacement of the initial conditions (the 492 HML is a system of ordinary differential equations) and thus directly comparable to the PSPF. 493 The scheme does not modify system states or parameters upstream of the streamflow gauge. The 494 results show better forecast skill of the HLM in terms of KGE for smaller basin scales and for 495 longer lead times. For example, from lead-time of one day to five days the forecast skill of the 496 HLM model starts bettering the performance of PSPF. With the addition of the simple data 497 assimilation framework, the HLM forecast skill is much enhanced over PSPF, as expected. The 498 HLM simulation assumes “perfect” rainfall forecast (i.e. limited only to the observational and 499 estimation errors). In reality, the rainfall forecasts are subject to considerable uncertainty, 500 especially for the longer lead times (e.g., Seo et al. 2018). 501 502 Figure 16: One-to-one comparison of a hydrologic model forecast skill with the PSPF. Red dots 503 show the comparison of KGE of the open loop HLM model with PSPF while blue dots show the 504 comparison of KGE of simple data assimilated open loop model with PSPF. This scatterplot 505 shows a potential of enhanced forecast skill upon adding modeling complexity to PSPF. 506 Acknowledgements 507 This study was funded by the Iowa Flood Center of the University of Iowa. The first author 508 also acknowledges partial support from the Rose & Joseph Summers endowment. The authors 509 are grateful to many colleagues at the IFC who facilitated the study by providing observational 510 and computational support. 511