Estimating watershed-scale storage changes from hourly discharge data in mountainous humid watersheds: toward a new way of dominant process modeling

The present study demonstrates and suggests a methodology for estimating watershed-scale storage changes from hourly discharge data in mountainous watersheds under a humid climate in Japan as the basis for watershed characterizations by combining a hydrograph separation and a storage-discharge formulizing method. Firstly, we separated hydrographs into different sub-components with respect to the time constants of the flow recession periods of the hydrograph by the filter-separation autoregressive method. Then we applied a storage-discharge formulizing method to the relationship between watershed-scale storage and discharge sub-components independently. As the result, we obtained a realistic estimate of annual watershed-scale storage change in the upper Abukuma River watershed. In addition, we compared our methodology with the non-linear reservoir modeling approach to explore the potential for extensive applications such as dominant process modeling, watershed classifications, and practical watershed management.


INTRODUCTION
Estimating watershed-scale storage is one of the fundamental questions of hydrology. Its importance becomes explicit when we recognize that it is the source of our water resources everywhere on the earth. Recently, hydrologists have been revisiting the estimation of storage (Spence, 2007(Spence, , 2010Kirchner, 2009;Soulsby et al., 2011;McNamura et al., 2011;Sayama et al., 2011;Peters and Aulenbach, 2011). However, because of our inability to know the bottom boundary conditions of watershed-scale storage, it is almost impossible to quantify the watershed-scale storage itself.
In practice, storage change, rather than storage itself, is more critical for us because our practical interests require knowledge of the increase or decrease of storages under flooding or drought events rather than storage amount in many cases. Recently, estimation of storage change became more realistic through the work of Kirchner (2009). He suggested that recession analysis in the nighttime rainless period after rainfall events would allow us to formulate the storage-discharge relationship by introducing traditional recession analysis by Brutsaert and Neiber (1977). Kirchner (2009) demonstrated his methodology in a UK watershed that can be modeled by a single nonlinear reservoir model, and hence modifications are necessary in watersheds with more complex hydrological processes as shown in Birkel et al. (2011), Krakauer and Temimi (2011), and . Teuling et al. (2010) extended Kirchner's approach by introducing variable recession rate as a function of discharge in a Swiss pre-alpine watershed, but their modification essentially does not fit to a watershed with seasonal storage and discharge as pointed out by McMillan et al. (2011).
Storage change estimations in a watershed with more complex hydrological processes are also going to be practical. Okazaki et al. (2013) suggested incorporating a hydrograph separation scheme of Hino and Hasebe (1984) with the storage estimation scheme by Kirchner (2009). The methodological effectiveness would be supported by the field-based study by Kosugi et al. (2011) who demonstrated how different aquifers within a watershed contribute to the form of hydrograph at the watershed outlet. Okazaki et al. (2013) succeeded in estimating a better hourly storage change than the original method by Kirchner (2009), however they did not describe their scientific contribution nor the further potential for applications well. In addition, the relationship between the methodologies of Okazaki et al. (2013) and storage-based models such as non-linear reservoir model (i.e. Nagai et al., 1982) is unexplained, which makes it difficult to understand the nature of the methodology and their scientific contributions.
To resolve the above unexplained problems in Okazaki et al. (2013), we traced their approach and discussed the relationship with a non-linear reservoir model (i.e. Nagai et al., 1982) that has been applied in many watersheds. In addition, we attempted to further translate the hydrological contributions of this approach for exploring potential applications of this approach to dominant runoff process modeling, watershed classifications, and practical watershed management.
Japan. The altitude of the station is 217 m surrounded by the Ou Mountains ranging from 600 to 1,900 m on the West and the Abukuma mountains ranging from 400 to 1,200 m on the East of the Abukuma main channel flowing to the North.
We downloaded hourly flow data at the Akutsu monitoring station and hourly precipitation data at the neighboring Koriyama precipitation monitoring station both for the period of 2005-2007 from the Water Information System maintained by the Ministry of Land, Infrastructure, Transport and Tourism (MLIT), Japan.

Hydrograph separation method
We applied the hydrograph separation method named "filter-separation auto-regressive method" by Hino and Hasebe (1984) to the hourly flow data after converting the units from m 3 /s to mm/h by dividing by the watershed area. Details of the method are described in Hino and Hasebe (1984), and we briefly explain it in the following.
The separation method calculates a slow flow component Q s by a convolution integral in equation (1). (1) In this equation, A, ω(τ), Q(t), t, τ, and n are respectively parameter, numerical filter, flow data, time, time starting from the beginning of filter, and the length of filter defined to be 5.0T c where T c is the time constant that is the inverse of β in equation (2).
The above equation (2) models only a recession period of a hydrograph. We evaluated the parameters α and β by fitting an exponential recession model in equation (2) to a target recession period, where t r become zero at the peak of a recession period. Hino and Hasebe (1984) introduced the numerical filter as in equation (3).
In equation (3), δ (> 2.0) is a damping factor, c 0 = (δ/T c ) 2 , and c 1 = δ 2 /T c . Hino and Hasebe (1984) used 2.5 for the parameter δ but we set it at 2.1 consistently. This difference in δ would not be effective, as we confirmed that the results of filtering were not sensitive to the magnitude of δ in our preliminary study. We applied the above hydrograph separation method repeatedly for obtaining multiple discharge sub-components that are characterized by the different magnitudes of T c (or β). We used only 2007 data, as the data in 2005-2006 were only used for removing the effects of initial conditions of the hydrograph separations. There must be the effect of seasonality of evapotranspiration (ET) in the result of hydrograph separation, but we did not consider the effect in our hydrograph separations leaving the effect of seasonal ET in the uncertainties in the time constant T c .

Estimation of storage change
We estimated watershed-scale storage change using the method suggested by Kirchner (2009). The following is a summary of the method. Kirchner (2009) firstly assumes that discharge Q is a function of storage S as per equation (4) where P is precipitation. This equation means that the S-Q relationship at a watershed scale can be approximated by the relationship between dQ/dt and Q in rainless nighttime recession periods without estimating ET. Kirchner (2009) further demonstrated that we can quantify watershed-scale storage by introducing the traditional recession model in Brutsaert and Nieber (1977) as an idealized approximation, to obtain the equation (6). (6) We can easily obtain the magnitude of the parameters a and b by finding the best-fit power function of equation (5) to the data in the rainless nighttime period. Kirchner (2009) integrated equation (6) as in equation (7), where S 0 is an integral constant that is equivalent to the storage level at which discharge becomes zero.
Although we applied an exponential function for modeling discharge recession period in hydrograph separation, with dQ/dt becoming an exponential function in equation (2), we introduced a power function between dQ/dt and Q. This may not allow for consistent discussion. Instead, it allows direct comparison with non-linear reservoir models (i.e. Nagai et al., 1982) and we decided to use the power function in equation (5) for the discharge recession model in the storage-discharge formulizing method. By comparing our model with non-linear reservoir models (i.e. Nagai et al., 1982) with model parameters k and p as in equation (8), We applied this method for all the discharge subcomponents calculated by the above hydrograph separation method of Hino and Hasebe (1984), where the relationship between dQ/dt and Q was plotted with the data in rainless nighttime recession periods for all the components. Nighttime in the present study was defined to be 19:00 to 6:00 (of the next day) at the local time.
In the present study, the time-variant overall watershedscale storage is calculated as, where k indicates the number of discharge sub-components. Figure 1a shows the 2007 hydrograph of the Akutsu monitoring station illustrated on semi-logarithmic axes with its characteristics exponential recession curves. It shows there are 4 different characteristic exponential recession curves that indicate this hydrograph consists of 4 discharge sub-components. We explored the number of the subcomponents by searching 3 to 5 characteristic exponential recession parts from a hydrograph on the assumption that the number of sub-components becomes 3 to 5 in Japanese watersheds as is often seen in hydrological modeling in Japan. The total number of sub-components was decided by applying 3 to 5 different linear segments on the recession curves plotted on semi-logarithmic axes. We summarized the magnitudes of the time constants T c (= 1/β) obtained from the recession curves for discharge sub-components in Table I. Figure 1b shows the result of hydrograph separations. We assumed that there existed 4 different discharge components and we applied the filter separation method 3 times. The figure shows discharge sub-components exhibit different roles in constituting the total discharge.

Hydrograph separation
Relationship between discharge and discharge recession rate Figure 2 shows the relationship between discharge Q and discharge recession rate dQ/dt on the recession curves in rainless nighttime periods for the 4 different discharge sub-components. The regression curves are obtained as the best fit power functions for discharge sub-components, which helps for obtaining the proportionality constants a and the exponents b of the power functions to be used in the equations (5)-(7). The magnitudes of a and b for each panel are also summarized in Table I.
Note that the diagrams in Figure 2 are all highly scattered and hence the obtained power law regression curves, which is the bases for the storage-discharge formulizations, cannot be considered representative for all the data in the diagrams. Hence this is one of the limitations of the current status of our method and we are currently exploring how to reduce the level of scatter.
The magnitudes of the exponent b of the regression curves in Figure 2 indicate the recession characteristics of discharge sub-components. The data for discharge subcomponents in panels (c) and (d) make us feel confident that the power functions better represent the characteristic scatters, whereas the power functions in panels (a) and (b) do not represent the scatters as closely, indicating potential directions for further improvement of this methodology.  Relationship between storage and discharge Figure 3a shows the relationships among watershed-scale storage and discharge for each discharge sub-component calculated by equations (7) and (9) with parameters a and b in Figure 2 and Table I. We can see that the gradients of discharge sub-components are higher in quicker flows and lower in slower flows, indicating faster flows are sensitive to changes of source storage amounts whereas slower flows are insensitive to them. This result supports the suitability of parameter magnitudes in the vertical series Tank model (Sugawara, 1995).
Although Okazaki et al. (2013) did not adequately consider storage-discharge relationships, they did provide valuable information on the watershed-scale storages. We plotted instantiations storage-discharge relationships for the year 2007 in Figure 3a, hence showing the annual variation of storages that generates different discharge subcomponents. The annual maximum storage was about 75 mm in the case of the quickest flow Q 4 . Those for Q 3 , Q 2 and Q 1 are respectively about 160 mm, 50 mm and 205 mm. If we consider that sediment disasters often occur if accumulated rainfall exceeds around 100 mm in Japanese mountainous watersheds, only the storage amounts that generate Q 4 and Q 3 would potentially explain the critical watershed scale storage amount that can cause sediment disasters.  Table I. The maximum and minimum storages became respectively 78.1 and 1.15 mm, which are different from our total storages and seem unrealistic as confirmed by . We can see that the annual storage change in the original method resembles that of precipitation and has no seasonal discharge variation in the original method, whereas our method exhibits seasonal discharge variation mainly caused by the storage change of the slowest discharge subcomponent. We think that our method is more natural and realistic because watershed-scale storage should be seasonal as rainfall has seasonality in Japan.

Adequacy of estimated storage amount
Here we want to recall that the storage-discharge relationships in Figure 3a are not representative enough to cover all the scatter diagrams in Figure 2. Currently we are investigating the reasons for the highly scattered diagrams and we cannot discuss the effects of uncertainties in the parameters a and b in Figure 2. Instead, we explored the effect of uncertainty in the magnitudes of recession time constant T c on the storage change estimations. The magnitudes of T c were calculated from visually fitted line segments on discharge data plotted on semi-logarithmic axes, and hence T c has uncertainties that can change the scatters in Figure 2 and estimates in Figure 3. Changing the level of uncertainties of T c from −10% to +10%, we recalculated the annual storage change where α in equation (2) was kept the same as those in Table I to focus on the effect of uncertainties in T c only. The results showed that the effects of uncertainties in T c appeared as storage level differences less than 10 mm in the cases of Q 2 and Q 3 . In the case of Q 1 , the effect of uncertainties in T c appeared as storage level differences less than about 120 mm. This difference would be smaller in actual modeling as α in equation (2) would also change when T c changes, but this result suggests estimation of T c for the slowest discharge sub-component is a key for reducing the uncertainties originating from the hydrograph separation method. More detailed results are summarized in Supplement Figure S1.

DISCUSSIONS AND CONCLUSIONS
The present study successfully demonstrates a methodology for estimating watershed-scale storage by combining a hydrograph separation method by Hino and Hasebe (1984) with the method by Kirchner (2009). Based on our results, we discuss below (1) the relationship with non-linear reservoir model, (2) adequacy of estimated watershed-scale storage change, and (3) possible future directions.

Relationship with non-linear reservoir model
The magnitudes of parameter b summarized in Table I allow comparison with the existing non-linear reservoir models. Parameter b ranged from 0.712 to 1.27 in our results, which corresponds to the range of the parameter p in the non-linear reservoir model from 0.73 to 1.29. Nagai et al. (1982) estimated parameter p to be 0.3 for laminar-type flows, 0.6 for Manning equation type flows, and 1.0 for Darcy type flows. Therefore, our estimates of b would be closer to Manning equation type flows or Darcy type flows, which seems reasonable. Yet, we should note that there is a structural difference between the non-linear reservoir model and our model. The non-linear reservoir model compared in the present study is a single storage model for a watershed, whereas ours is a multiple storage model and our parameters are evaluated for each discharge subcomponent. Hence, quantitative parameter comparisons would not be suitable between the two models. Rather, the authors consider that all the parameters b are close to 1.0 and hence discharge sub-components are all approximated as liner-reservoir models at least in our study area as we can see in Figure 3a. This has also been pointed out by Hino and Hasebe (1984).
For a quantitative comparison between our model and non-linear reservoir models, we need to compare watershedscale storage-discharge relationships between the two models as in Figure 4 that has not been explored so far. In the figure, we drew the watershed-scale storage and discharge relationship as a cumulative form assuming that watershed scale storage is to be filled up from the slowest discharge sub-component. Hence from the maximum point of the storage-discharge relationship for Q 1 , we added the storage-discharge relationships for Q 2 on top of the  (Nagai et al., 1982) where parameter k was kept constant at 100 and parameter p was changed as 0.3 for laminar type flows (SFM p = 0.3), 0.6 for Manning type flows (SFM p = 0.6), and 1.0 for Darcy type flows (SFM p = 1.0) maximum of the storage-discharge point for Q 1 . This figure shows that our storage-discharge relationship for Q 1 seems similar to the Darcy type flow (p = 1.0) and those for Q 2 , Q 3 , and Q 4 seem closer to laminar flow (p = 0.3). In addition, we can see that the Manning type flow (p = 0.6), among the 3 kinds of non-linear reservoir models (p = 0.3, 0.6, and 1.0), should be the most suitable relationship to explain our overall watershed-scale storage-discharge relationship. We can also see that, as a whole system, our storage-discharge relationship seems nonlinear as it is comprised of 4 different linear storage-discharge sub-systems.

Adequacy of estimated watershed-scale storage change
Currently it is impossible to evaluate watershed-scale storage amount because of the lack of our ability to specify the lower bottom boundary of the aquifer. Yet, we can discuss the adequacy of the annual watershed-scale storage change. Assuming that storage is 0 mm when discharge is 0 mm/h, we estimated maximum storage depth to be 327 mm in 2007. Note that the meaning of "storage" corresponds to the easily movable water storage that eventually contributes to river discharge. The cumulative precipitation of the heavy rainfall event of Tochigi-Fukushima of 1998 in our study area ranged from about 600 to over 1200 mm. This event caused landslide disasters in Shirakawa City, hence storage depth would have been nearly full. Therefore we think that the estimated maximum storage of 327 mm in 2007 is in a reasonable range as an estimate for a non-severe-flood year. Yet, it must be noted that we need further analysis of many different records of flood years to achieve more precise results.

Possible future directions
We achieved a realistic estimation of watershed-scale storage and storage change in a watershed with complex hydrological processes by incorporating a hydrograph separation scheme by Hino and Hasebe (1984) with storagedischarge formulization by Kirchner (2009). The next step we need to tackle is verifying the adequacy of watershedscale storage estimated by our methodology. Applying this methodology to heavy rainfall events that generated flood inundations or sediment disasters would be helpful for evaluating the maximum storage amount as mentioned above, which is demonstrated in Yokoo et al. (2013). Another possible application would be of drought situations. If we can use Figure 3 with confidence, we would be able to discuss the discharge sustainability under drought conditions based on the residual storage in each aquifer that generates total discharge. Such application studies would help us to identify directions for further improvements.
Another possible way forward would be using this methodology for dominant hydrological process modeling of a watershed considering the issue of prediction of ungauged basins (Sivapalan et al., 2003). Our methodology allows hydrograph separations and quantification of multiple storage levels, hence we can use this information as a kind of "soft data" (Seibert and McDonnell, 2002) which constrains model structure and parameter values. As a preliminary study, we have already succeeded in identifying structure and parameters in the case of the Tank model (Sugawara, 1995), which dramatically reduces optimizing parameters from more than 10 parameters (Yokoo et al., 2001;Yokoo and Kazama, 2012) to only 3 parameters. We intend to report the results in a separate paper.
The authors believe this methodology has potential as a powerful tool for watershed classification. We can obtain the number of watershed-scale storages and their capacities which could provide a measure suitable for classifying a range of the world's watersheds.