A Polynomial Method Approximating S-Curve with Limited Availability of Reliable Rainfall Data

A methodology named the step response separation (SRS) method for deriving S-curves solely from the data for basin runoff and the associated instantaneous unit hydrograph (IUH) is presented. The SRS method extends the root selection (RS) method to generate a clearly separated S-curve from runoff incorporated in mathematical procedure utilizing the step response function. Significant improvements in performance are observed in separating the S-curve with rainfall. A procedure to evaluate the hydrologic stability provides ways to minimize the oscillation of the S-curve associated with the determination of infiltration and baseflow. The applicability of the SRS method to runoff reproduction is examined by comparison with observed basin runoff based on the RS method. The SRS method applied to storm events for the Nenagh basin resulted in acceptable S-curves and showed its general applicability to optimization for rainfall-runoff modeling.


Introduction
This paper is the sequel of a previously suggested methodology for smoothing the oscillatory S-curve [1], in which a model was presented for determining the oscillationreduced S-curve and associated instantaneous unit hydrograph (IUH) model using the Savitzky-Golay smoothing and differentiation filter.
The concept of the S-curve has been widely applied in rainfall-runoff analysis since its introduction in the unit hydrograph (UH) method [2,3]. The S-curve represents a direct runoff due to the effective rainfall applied over an infinite time, and its intensity is one-unit depth per unit of time. The role of the S-curve is to alter a UH of arbitrary duration into a UH of desired duration. The conventional method to generate an S-curve requires surface runoff data and corresponding effective rainfall. However, estimating accurately the effective rainfall for S-curve generation is difficult mainly because of problems related to the estimation of a real rainfall and to the determination of the reduction in total rainfall due to infiltration, which is affected by the uncertain soil moisture condition of basins and groundwater behavior [4,5]. Mathematically, these difficulties result in a hydrologically unstable system. Following, a problem arises in determining the S-curve of a linear rainfallrunoff system, since an oscillatory function is likely to be obtained. This problem can be circumvented by developing a method deriving an S-curve without using rainfall data.
The main objectives of this study are to estimate an S-curve without observed rainfall data and to investigate the hydrologic stability of rainfall-runoff data by using the derived S-curve. The approach adopted in this study relies on the root selection (RS) method [6]. This study names the suggested procedure as the step response separation (SRS) method. The RS method only requires runoff data from a single storm of the basin of interest to give a UH estimate [3]. In the RS method, the UH is determined based on the shape traced by the complex roots of the runoff polynomial on an Argand diagram, which yields some considerable applications. For example, the method was adopted in developing a design hydrograph in an arid basin [7], in identifying the quick and the slow components of runoff [4], and in synthesizing the Clark IUH [8]. However, very little research has been conducted on the RS method despite its advantages, mainly because of practical reasons since very high-order polynomials reflecting direct runoffs should be manipulated to derive UH ordinates [9]. Another reason is that the method of using an Argand diagram is subjective and cumbersome. Therefore, even when there is direct runoff data, synthetic UHs are often used in practice. However, in this process, important information about watersheds and rainfall-runoff may be lost. This important information includes the runoff characteristics and hydrological stability of the basin. Therefore, deriving UH through polynomial analysis is hydrologically meaningful.
In the present study, this paper develops the procedure for estimating an S-curve solely from basin runoff data by improving the defects of the existing method and examines the hydrologic stability of the resulting rainfall-runoff system. To that goal, this study uses cumulative runoff in the RS method. The applicability of the proposed methodology is evaluated using rainfall-runoff data for the Nenagh River at Clarianna, Ireland [10,11]. The performance of the derived S-curve and corresponding IUH is evaluated by the standard errors and Nash-Sutcliffe criterion. The proposed method is found to be more objective and distinct to generate hydrologic responses and be capable of checking the hydrologic stability of rainfall-runoff system.

Procedure of the Proposed Method
The impulse response function of a linear system indicates the system response to an instantaneous impulse of unit volume of input applied at time t = 0. In hydrology, it is assumed that basins behave as linear system, and the impulse response function represents an IUH for the system-theoretical modeling of rainfall-runoff processes. Denoting the direct runoff by q(t), the effective rainfall by x(t), and the IUH by h(t) for stationary and casual systems, q(t) can be decomposed into IUH through the convolution integral: The unit step response function defines the response to a constant unit input for an infinite time. In hydrology, the unit step response function corresponds to the S-curve caused by a continuous effective rainfall of unit intensity. Thus, the S-curve can be obtained by substituting x(t) = µ(t) in Equation (1), which leads to: where s(t) is an S-curve, and µ(t) is the unit step function defined as follows: In Equation (2), s(t) is the response to effective rainfall applied over an infinite time, and its intensity is one-unit depth per unit of time (e.g., 1 cm/h), which is equivalent to the unit step response defined as the output of a system due to the unit step function µ(t).
This study suggests a convolution approach determining s(t) without using rainfall data. For this purpose, µ(t) is applied in Equation (1), which leads to where r(t) denotes a cumulative function of the runoff (CFR). Using the properties of convolution (commutative law or associative law), Equation (4) is rearranged as: Thus, the S-curve s(t) can be determined by solving an inverse problem if r(t) is given. When the effective rainfall x(t) is given as a hyetograph, the corresponding CFR can be represented as the discrete convolution of the rainfall hyetograph and S-curve. In the transformation technique for discrete time function, the Z-transform [12] is preferably used. Based on the Z-transformation theory, Equation (5) can be expressed in term of polynomials in z −1 , where R z −1 , S z −1 , and X z −1 are the Z-transforms of r(t), s(t), and x(t), respectively. When CFR data {r(iT), i = 0 · · · P} are collected repeatedly at fixed time interval T, R z −1 in Equation (6) can be formulated as: where r(iT)z −1 represents the Z-transform of r(iT) calculated at t = iT. Equation (7) is a P-order polynomial in z −1 , which implies the existence of P roots. These roots consist of pairs of conjugate complex numbers so that they can be represented on a plane called the Argand diagram [12]. As with the existing RS method, the method of this study also ensures that the roots of Equation (7) depict a circular pattern on the Argand diagram with some points showing abnormal feature. The roots of the S-curve polynomial S z −1 , or the step response function, are selected based on the circle pattern, while the roots of the rainfall polynomial X z −1 are regarded as remaining roots. Using similar notation, S z −1 and X z −1 can be represented by the polynomials shown in Equations (8) and (9): where N is the number of roots of the S-curve polynomial, where M is the number of remaining roots. According to the linear convolution theory, the following relationship between the polynomials is satisfied: This study names the above-mentioned procedure as the SRS method.

Derivation of CFR Polynomial
To obtain the CFR polynomial R z −1 of Equation (7), the runoff polynomial Q z −1 , that is the Z-transform of q(t), is also calculated according to Equation (4). When runoff data {q(iT), i = 0 · · · P} are collected repeatedly at fixed time interval T, Q z −1 can be formulated as: where q(iT)z −1 represents the Z-transform of q(iT) sampled at t = iT. The Z-transform of the unit step function µ(t) is defined as follows [13]: Thus, Equation (4) with respect to CFR becomes: whereŜ z −1 is an estimate of S z −1 . However, when using Equation (13), one may face the problem of factoring the polynomial. As already mentioned, S z −1 indicates the S-curve hydrograph, which implies thatŜ z −1 should be a real number coefficient polynomial in z −1 . To satisfy this condition, 1 − z −1 should be a factor of the polynomial Q z −1 . However, in general, it is not the case, because the polynomial division of Equation (13) yields a remainder. So, this research uses a power expansion technique to approximate the rational function R z −1 . In fact, this expansion technique is mathematically plain because the divisor 1 − z −1 constructs R z −1 by adding together the coefficients of Q z −1 lagged by one time step. Hence, the following equation for CFR can be formulated: Therefore, Equation (14) can replace Equation (7). Solving for R z −1 = 0, Equation (14) yields P roots. These P roots are drawn on the Argand diagram, and the pattern of these roots is investigated.

Inverse Transformation of S-Curve into Time Domain
As previously explained, similar to the RS method, the method of this study determines the S-curve polynomial S z −1 by selecting the complex roots based on the a priori expectation that the roots will draw a circle in the Argand diagram. The P roots obtained by solving Equation (14) are drawn on the Argand diagram. Among them, N roots pertain to the S-curves polynomialŜ z −1 and M roots pertain to X z −1 . Hence, the S-curve can be written in the form of the polynomial: where the coefficients s 1 , s 2 , . . . , and s N represent the ordinates of the S-curve. Note that these coefficients are not the exact values but the relative magnitude of each coordinate of the S-curve. Therefore, these values should be adjusted using the equilibrium of the S-curve. Based on S-curve theory, the coefficient s N of the highest term of Equation (15) should attain the equilibrium discharge q eq occurring at the base time of UH, and the value is determined from: where A is the watershed area (km 2 ); d is the unit rainfall depth (1 cm); and, T is the sampling interval. Adjusted with q eq by a linear programming technique, the ordinates of the S-curve can be defined as follows:

Model Performance Indices
Two different criteria are considered for evaluating the model performance of the SRS method. To do this, this study uses the derived IUHs for all storms by means of numerical differentiation of S-curves and their reproduced runoff.
As a model performance criterion, the standard errors (STDER) [14] are evaluated. STDER is defined by the following equation: where h n is the IUH ordinate at the nth time interval; h re f ,n is the ordinate of the referenced IUH at the same time; and w n is the weight given by: where h avg is the average value of the referenced IUH ordinates. A low value of STDER implies a good model, and the zero value of STDER means that the computed IUH perfectly matches with the referenced IUH. In addition, the performance of the SRS method is evaluated by the Nash-Sutcliffe criterion based on the reproduced runoff [15] given in Equation (20).
where E is the Nash-Sutcliffe model efficiency coefficient: J is the number of observations; q j is the observed runoff at time j;q j is the predicted runoff at time j determined from the proposed IUH; and q is the mean of the observed runoff. The closer E is to 100, the better the performance of the model [16]. The peak relative error, QB, is the difference in magnitude between the observed and estimated peaks and is used for model performance analysis [17].
where q p is the peak value of observed runoff; andq p is the peak of estimated runoff. The closer QB is to 0, the better the performance of the model.

Results
The proposed methodology using CFR and the Argand diagram is applied to the estimation of the S-curveŜ z −1 . The considered basin is the Nenagh River at Clarianna, Ireland, with an area of 295 km 2 [6]. The catchment of Clarianna has four rain gauges reading daily records and one recording continuously. The present study refers to the data for five of 22 storms on the Nenagh River basin, which were examined by using the least squares analysis [10]. The considered five storm events present strong difficulty in identifying the roots of CFR and are referred to as events 5, 13, 14, 15, and 20. Note that storm event 13 will be used for identifying the sensitivity of the proposed method. For the selected events, the peaks of the direct runoffs recorded a range from 13.19 to 23.52 m 3 s −1 when the measurements are made at 3 h (= T) interval. The equilibrium discharge q eq is 273 m 3 s −1 . The major features of the selected storm events are listed in Table 1. The proposed method determines the S-curve S z −1 by selecting the complex roots based on the a priori expectation that the roots will describe a circle in the Argand diagram for the roots of the CFR polynomial. The Argand diagrams of storm events 5, 14, 15, and 20 are plotted in Figure 1. The root separation procedure by the SRS method developed by modifying the existing RS method [6] is as follows: (1) select the S-curve roots based on the circle pattern; (2) construct the S-curve using the so-separated roots; and (3) combine the remaining roots to form an estimate of effective rainfall. The complex roots for the direct runoff polynomial for each event plotted in Figure 1 are relevant to the S-curve ( ) and rainfall ( ). In Figure 1, the roots of ( ) represented by dots are spread at a practically equal angle from each other and describe almost perfect circles with a slightly increasing radius from left to right. Approximating the polynomial for ( ) is achieved by combining the roots on these circles. To check the validity of the approximation, ( ) is illustrated in Figure 2.  The root separation procedure by the SRS method developed by modifying the existing RS method [6] is as follows: (1) select the S-curve roots based on the circle pattern; (2) construct the S-curve using the so-separated roots; and (3) combine the remaining roots to form an estimate of effective rainfall. The complex roots for the direct runoff polynomial for each event plotted in Figure 1 are relevant to the S-curveŜ z −1 and rainfall X z −1 . In Figure 1, the roots ofŜ z −1 represented by dots are spread at a practically equal angle from each other and describe almost perfect circles with a slightly increasing radius from left to right. Approximating the polynomial for S z −1 is achieved by combining the roots on these circles. To check the validity of the approximation,Ŝ z −1 is illustrated in Figure 2. The root separation procedure by the SRS method developed by modifying the existing RS method [6] is as follows: (1) select the S-curve roots based on the circle pattern; (2) construct the S-curve using the so-separated roots; and (3) combine the remaining roots to form an estimate of effective rainfall. The complex roots for the direct runoff polynomial for each event plotted in Figure 1 are relevant to the S-curve ( ) and rainfall ( ). In Figure 1, the roots of ( ) represented by dots are spread at a practically equal angle from each other and describe almost perfect circles with a slightly increasing radius from left to right. Approximating the polynomial for ( ) is achieved by combining the roots on these circles. To check the validity of the approximation, ( ) is illustrated in Figure 2.  In Figure 2, the estimated S-curves for the four storms were derived by the SRS method and are compared with those obtained by the conventional method in which the 3 h S-curves are computed by successively adding the coordinates of 3 h UHs lagged by the duration (=3 h) until 24 time steps. The 3 h UH of each storm is calculated by means of the least square method. From the figure, it is seen that the SRS method generates smooth S-curves when using the roots ofŜ z −1 represented in the Argand diagram in Figure 1. Moreover, the coordinates of S-curves using the SRS method agree well with those of the conventional method in all storms except for storm event 20. From the comparison of the coefficients of theŜ z −1 polynomial with the other S-curves, this study can infer that a nonlinear dynamic system of rainfall-runoff structure to obtain accurate surface runoff coordinates and relevant rainfall values gave rise to this inconsistent behavior.
Since numerical comparison is difficult, Figure 3 compares the IUH ordinates for the four selected storms [1,11]. In Figure 2, the estimated S-curves for the four storms were derived by the SRS method and are compared with those obtained by the conventional method in which the 3 h S-curves are computed by successively adding the coordinates of 3 h UHs lagged by the duration (=3 h) until 24 time steps. The 3 h UH of each storm is calculated by means of the least square method. From the figure, it is seen that the SRS method generates smooth S-curves when using the roots of ( ) represented in the Argand diagram in Figure 1. Moreover, the coordinates of S-curves using the SRS method agree well with those of the conventional method in all storms except for storm event 20. From the comparison of the coefficients of the ( ) polynomial with the other S-curves, this study can infer that a nonlinear dynamic system of rainfall-runoff structure to obtain accurate surface runoff coordinates and relevant rainfall values gave rise to this inconsistent behavior.
Since numerical comparison is difficult, Figure 3 compares the IUH ordinates for the four selected storms [1,11]. In theory, the ordinate of an IUH at time t is the slope of the S-curve of intensity 1 cm/h at the corresponding time. If the time interval of ordinates of the S-curve is infinitesimal, the slope is equivalent to the differentiation. In this study, the differentiation of ( ) is performed by the digital filtering method as well as the method of a prior research [1], which uses the Savitzky-Golay smoothing and differentiation filter. These IUHs based on the SRS method were compared with Nash IUHs reported in the literature [11]. The Nash IUH is represented in Equation (22).
where ℎ( ) is the IUH ordinates at time t, n, and K are parameters that define the shape and scale of the Nash IUH. The Nash IUH requires a complete set of rainfall-runoff data to obtain the model parameters.
In Figure 3, the values of the time-to-peak of IUH based on the SRS method appear to be very similar (within single time step of calculation) to those of the Nash IUHs in all selected cases. The rising and falling limbs of the proposed IUH are underestimated near the peak of the Nash IUH for each storm event except storm event 20. To check the performance index of STDER, the Nash IUH is adopted as the reference IUH. The STDER of In theory, the ordinate of an IUH at time t is the slope of the S-curve of intensity 1 cm/h at the corresponding time. If the time interval of ordinates of the S-curve is infinitesimal, the slope is equivalent to the differentiation. In this study, the differentiation ofŜ z −1 is performed by the digital filtering method as well as the method of a prior research [1], which uses the Savitzky-Golay smoothing and differentiation filter. These IUHs based on the SRS method were compared with Nash IUHs reported in the literature [11]. The Nash IUH is represented in Equation (22).
where h(t) is the IUH ordinates at time t, n, and K are parameters that define the shape and scale of the Nash IUH. The Nash IUH requires a complete set of rainfall-runoff data to obtain the model parameters.
In Figure 3, the values of the time-to-peak of IUH based on the SRS method appear to be very similar (within single time step of calculation) to those of the Nash IUHs in all selected cases. The rising and falling limbs of the proposed IUH are underestimated near the peak of the Nash IUH for each storm event except storm event 20. To check the performance index of STDER, the Nash IUH is adopted as the reference IUH. The STDER of each storm is listed in Table 2 along with the adopted parameter values to set up the Nash IUH [11].  Table 2 indicates that when the SRS is applied to each storm, the estimated IUHs are generally in good agreement. Interestingly, storm event 20 shows the smallest STDER even if it provided rather poor accuracy in estimating the peak value.
A plot showing the comparison between the observed and reproduced discharges for the four storms is presented in Figure 4. each storm is listed in Table 2 along with the adopted parameter values to set up the Nash IUH [11].  Table 2 indicates that when the SRS is applied to each storm, the estimated IUHs are generally in good agreement. Interestingly, storm event 20 shows the smallest STDER even if it provided rather poor accuracy in estimating the peak value.
A plot showing the comparison between the observed and reproduced discharges for the four storms is presented in Figure 4. When the SRS method is applied to measured rainfall-runoff data, the estimated UHs should be consistent with the UHs resulting from the conventional method to provide reproduced flows close to the measured runoff. An overview of runoff reproduction by SRS is presented to appreciate the advantage of the SRS method. For given storm events, the calculated results listed in Table 2 are compared through the Nash-Sutcliffe model efficiency coefficient formulated in Equation (20) and the peak relative error formulated in Equation (21). It can be clearly observed from Figure 4 as well as Table 1   When the SRS method is applied to measured rainfall-runoff data, the estimated UHs should be consistent with the UHs resulting from the conventional method to provide reproduced flows close to the measured runoff. An overview of runoff reproduction by SRS is presented to appreciate the advantage of the SRS method. For given storm events, the calculated results listed in Table 2 are compared through the Nash-Sutcliffe model efficiency coefficient formulated in Equation (20) and the peak relative error formulated in Equation (21). It can be clearly observed from Figure 4 as well as Table 1 that the UH based on the SRS method reproduced runoff responses very close to the observed runoff hydrograph, while the reproduced runoff responses are not consistent with the observed ones for storm events 5 and 20 by showing E = 77.6% and 40.2%, respectively. The reproduced runoff generally underestimates within a narrow range, which is between 0.007 and 0.420. However, the simulated runoff using the Nash model provided in prior research [11] shows virtually consistent results with the runoff responses using the SRS method. The reproduced runoffs for four storms were computed for the Nash model, and the results were compared with those of the SRS method in Figure 4 and listed in Table 2. The figure shows that the reproduced ordinates by the SRS method and Nash model agree well with each other. Interestingly, although it is not a remarkable difference in values, the SRS method shows a better or comparable E value than the Nash model, even though the SRS method is a restricted method not using rainfall information. However, the SRS method is far less sensitive to fit the peak ordinate of measured runoff than the Nash model when comparing QBs in Table 2. As long as accurately measured effective rainfall and direct runoff are used, the reproduced runoff will be consistent with the result of existing UH model, as seen in the case of storm event 14. The discrepancy of storm event 15 is due to the oscillation of the S-curve. This will be discussed in detail in the discussion section. The problem of event 20 can be resolved if a delay is given to S-curve. In Figure 1, the point z −1 = 0 represents unit delay, and this delay is equally applied to each event. Therefore, this discrepancy problem of storm event 20 appears to be caused by the nonlinearity of this event. Thus, the inconsistency problem of storm events 5 and 20 is not inherent to the proposed method itself, but it is associated with the fundamental problem related to the nonlinearity of the system, the base-flow separation, and infiltration for these storms. It is an important aspect of the SRS method which can give information on the usefulness level of the rainfall-runoff data we are dealing with.

Comparison of Argand Diagrams
The conventional RS method [6] determines the roots of the UH polynomial H z −1 , in which H z −1 is the Z-transform of h(t) in Equation (11). In the present method, the selection of the UH roots is performed based on the a priori expectation that the roots will draw a skewed circle in the Argand diagram. The Argand diagrams of storm events 5, 14, 15, and 20 are plotted in Figure 5. ones for storm events 5 and 20 by showing E = 77.6% and 40.2%, respectively. The reproduced runoff generally underestimates within a narrow range, which is between 0.007 and 0.420. However, the simulated runoff using the Nash model provided in prior research [11] shows virtually consistent results with the runoff responses using the SRS method. The reproduced runoffs for four storms were computed for the Nash model, and the results were compared with those of the SRS method in Figure 4 and listed in Table 2 The figure shows that the reproduced ordinates by the SRS method and Nash model agree well with each other. Interestingly, although it is not a remarkable difference in values the SRS method shows a better or comparable E value than the Nash model, even though the SRS method is a restricted method not using rainfall information. However, the SRS method is far less sensitive to fit the peak ordinate of measured runoff than the Nash model when comparing QBs in Table 2. As long as accurately measured effective rainfall and direct runoff are used, the reproduced runoff will be consistent with the result of existing UH model, as seen in the case of storm event 14. The discrepancy of storm event 15 is due to the oscillation of the S-curve. This will be discussed in detail in the discussion section. The problem of event 20 can be resolved if a delay is given to S-curve. In Figure  1, the point z −1 = 0 represents unit delay, and this delay is equally applied to each event Therefore, this discrepancy problem of storm event 20 appears to be caused by the nonlinearity of this event. Thus, the inconsistency problem of storm events 5 and 20 is not inherent to the proposed method itself, but it is associated with the fundamental problem related to the nonlinearity of the system, the base-flow separation, and infiltration for these storms. It is an important aspect of the SRS method which can give information on the usefulness level of the rainfall-runoff data we are dealing with.

Comparison of Argand Diagrams
The conventional RS method [6] determines the roots of the UH polynomial ( ) in which ( ) is the Z-transform of ℎ( ) in Equation (11). In the present method, the selection of the UH roots is performed based on the a priori expectation that the roots will draw a skewed circle in the Argand diagram. The Argand diagrams of storm events 5, 14 15, and 20 are plotted in Figure 5.  All the storm events in Figure 5 show a number of roots of the polynomial H z −1 lying outside the expected ideal skewed circle. The skewed circles are not easily identifiable for storm events 15 and 20 in Figure 5c,d, respectively. Furthermore, the abnormal roots, representing the roots of effective rainfall polynomial X z −1 can be over-counted according to the researcher's subjective decision. These over-counted roots sometimes produce a higher degree polynomial X z −1 than the actual one, and they may lead to unrealistic effective rainfall distribution. Furthermore, the relevant UH polynomial is never stable unless smoothing methods are additionally introduced into the RS procedure. In contrast, the proposed method differs from the conventional method in that SRS makes it possible to separate clearly X z −1 from runoff, as shown in Figure 1. In fact, the SRS method automatically separates X z −1 from R z −1 prior to the analysis. Naturally, the coordinates of roots for X z −1 obtained by SRS are completely included in the roots of the RS method. This explains that SRS does not affect X z −1 of the RS method. Therefore, it is reasonable to suppose that SRS is more objective by reducing some ambiguity in separating transfer functions from runoff data.
Unusual rainfall pattern, errors in the baseflow separation method, and measurement error may result in a pattern of chaotically scattered roots. As a result, the S-curve oscillates rather than approaching an equilibrium value and can produce largely erroneous ordinates at the tail. As mentioned earlier, the circle pattern of roots is not more easily identifiable for storm event 15, as seen in Figure 5c. This leads to the oscillations of the S-curve. In this study, an attempt was made to stabilize the oscillations occurring at the ends of the S-curve using a simple linear programming method, but it did not give better results than those shown in Figure 2c. Therefore, the S-curve needs to be smoothened graphically before taking derivatives to give better results than the IUH and runoff predictions shown in Figures 3c and 4c, respectively.

Examples for Different River Basins
To evaluate the generality of the SRS method, additional analyses for the storms of other basins were carried out. For this, Argand diagrams for direct runoff of the three different river basins using the SRS method were compared to the results obtained by the RS method. The basins considered here are dealt in the previous study [1] by the authors, and they show a wide range of watershed areas. The first result is from the storm event occurring in the Almond and Almondell catchment (229 km 2 ), U.K. [18]; the second result is obtained using a storm event in the Shaol Creek watershed (11.25 km 2 ), USA [5]; the third result is obtained based on the runoff data of the North Potomac River watershed (2266 km 2 ), USA [19]. The complex roots for each basin resulting from both SRS and RS methods are plotted on Argand diagrams in Figure 6.
From Figure 6, it appears clearly that the SRS method leads to an apparently better separation of hydrologic response than the RS method in each considered basin, as observed previously in the Nenagh basin. As a result, the roots corresponding to the S-curve and rainfall can be readily identified in each basin when using SRS. Another advantage is that the specific hydrologic response feature of each basin can be readily distinguished from each other when the SRS method is used. The left panels of Figure 6 show that the position of the circular pattern of roots differs slightly in each case. This implies that there exists a specific root pattern of hydrologic response in each basin. This can be identified in the Nenagh basin cases by showing the roots pattern of S-curves to be virtually identical to each other, as illustrated in Figure 1. In fact, a similar result has already been revealed in a previous study [6] using RS. However, SRS provides better results than using the RS method.

Consideration of Stability Using Argand Diagram
In the frequency domain, the system response functions, such as S-curve or UH, are given by the quotient obtained when dividing the input (rainfall) by the output (CFR or runoff). Therefore, the coefficients of these functions are positive and real. The system response (S-curve) can be given by rearranging Equation (6), which leads to For a stable hydrological system, the zero-crossings of the polynomial denominator must be equal to those of the polynomial numerator; otherwise, the system cannot be stable because of the occurrence of oscillations in the system responses. The zero-crossings of the numerator and denominator polynomials in Equation (23) are conjugate complex numbers. In the SRS or RS method, the polynomial denominator of Equation (23) is a mathematical rainfall polynomial ( ) of which the zero-crossings are those of the polynomial numerator. In a previous study [20], this ( ) was used as a replacement of the observed rainfall polynomial ( ) to get complete stability of the system. However, it is hardly hydrologically justifiable unless the roots of ( ) are close to those of ( ). However, instead of using this replacement method, using the least square method with an average rainfall of ( ) and ( ) seems to be more rational if hydrologic restrictions such as infiltration assumptions or continuity condition

Consideration of Stability Using Argand Diagram
In the frequency domain, the system response functions, such as S-curve or UH, are given by the quotient obtained when dividing the input (rainfall) by the output (CFR or runoff). Therefore, the coefficients of these functions are positive and real. The system response (S-curve) can be given by rearranging Equation (6), which leads to For a stable hydrological system, the zero-crossings of the polynomial denominator must be equal to those of the polynomial numerator; otherwise, the system cannot be stable because of the occurrence of oscillations in the system responses. The zero-crossings of the numerator and denominator polynomials in Equation (23) are conjugate complex numbers. In the SRS or RS method, the polynomial denominator of Equation (23) is a mathematical rainfall polynomial X math z −1 of which the zero-crossings are those of the polynomial numerator. In a previous study [20], this X math z −1 was used as a replacement of the observed rainfall polynomial X obs z −1 to get complete stability of the system. However, it is hardly hydrologically justifiable unless the roots of X math z −1 are close to those of X obs z −1 . However, instead of using this replacement method, using the least square method with an average rainfall of X math z −1 and z −1 seems to be more rational if hydrologic restrictions such as infiltration assumptions or continuity condition are a priori considered. In this case, stability can be achieved by using S-curve smoothing [1,20].
In order to investigate the hydrologic restrictions, the effect of a small difference in the runoff value on the hydrologic stability is examined. To do this, this study selected an event where the real part of the rainfall roots is negative. For example, for storm event 13 in Table 1, increasing the initial non-zero runoff value (=ordinate at t = 3 h) by 1.1 m 3 s −1 from 1.9 to 3.0 m 3 s −1 produces the circle pattern of the runoff polynomial shown in Figure 7. The increased runoff value is at most 4% of the peak flow rate (25.32 m 3 s −1 ), but it can be seen from Figure 7 that a hydrologically sensitive change occurs.
Water 2021, 13,3447 12 of 14 are a priori considered. In this case, stability can be achieved by using S-curve smoothing [1,20]. In order to investigate the hydrologic restrictions, the effect of a small difference in the runoff value on the hydrologic stability is examined. To do this, this study selected an event where the real part of the rainfall roots is negative. For example, for storm event 13 in Table 1, increasing the initial non-zero runoff value (=ordinate at = 3 h) by 1.1 m 3 s −1 from 1.9 to 3.0 m 3 s −1 produces the circle pattern of the runoff polynomial shown in Figure  7. The increased runoff value is at most 4% of the peak flow rate (25.32 m 3 s −1 ), but it can be seen from Figure 7 that a hydrologically sensitive change occurs. By performing this alteration, the root patterns of the S-curve polynomial ( ) become closer to a smoother circle for the initial runoff of 3.0 m 3 s −1 , as shown in Figure  7a,b. Even so, this small variation of initial runoff value does not affect significantly the root pattern of ( ). This is thought to be because a large number of S-curve ordinates share the alteration of runoff. Accordingly, both IUHs are readily identical, as shown in Figure 8a. However, the result for rainfall polynomial ( ) is quite different from that of the S-curve. The coefficient of the first-order term in polynomial ( ) changed from 0.3322 to −0.0640. As a result, the rainfall intensity after altering runoff has a negative value, as can be seen in Figure 8b. This change in sign is due to the shift of the positions of the roots of ( ) on the positive real axis resulting from the alteration of the initial runoff value. By performing this alteration, the root patterns of the S-curve polynomial H z −1 become closer to a smoother circle for the initial runoff of 3.0 m 3 s −1 , as shown in Figure 7a,b. Even so, this small variation of initial runoff value does not affect significantly the root pattern of S z −1 . This is thought to be because a large number of S-curve ordinates share the alteration of runoff. Accordingly, both IUHs are readily identical, as shown in Figure 8a. are a priori considered. In this case, stability can be achieved by using S-curve smoothing [1,20]. In order to investigate the hydrologic restrictions, the effect of a small difference in the runoff value on the hydrologic stability is examined. To do this, this study selected an event where the real part of the rainfall roots is negative. For example, for storm event 13 in Table 1, increasing the initial non-zero runoff value (=ordinate at = 3 h) by 1.1 m 3 s −1 from 1.9 to 3.0 m 3 s −1 produces the circle pattern of the runoff polynomial shown in Figure  7. The increased runoff value is at most 4% of the peak flow rate (25.32 m 3 s −1 ), but it can be seen from Figure 7 that a hydrologically sensitive change occurs. By performing this alteration, the root patterns of the S-curve polynomial ( ) become closer to a smoother circle for the initial runoff of 3.0 m 3 s −1 , as shown in Figure  7a,b. Even so, this small variation of initial runoff value does not affect significantly the root pattern of ( ). This is thought to be because a large number of S-curve ordinates share the alteration of runoff. Accordingly, both IUHs are readily identical, as shown in Figure 8a. However, the result for rainfall polynomial ( ) is quite different from that of the S-curve. The coefficient of the first-order term in polynomial ( ) changed from 0.3322 to −0.0640. As a result, the rainfall intensity after altering runoff has a negative value, as can be seen in Figure 8b. This change in sign is due to the shift of the positions of the roots of ( ) on the positive real axis resulting from the alteration of the initial runoff value.  However, the result for rainfall polynomial X z −1 is quite different from that of the S-curve. The coefficient of the first-order term in polynomial X z −1 changed from 0.3322 to −0.0640. As a result, the rainfall intensity after altering runoff has a negative value, as can be seen in Figure 8b. This change in sign is due to the shift of the positions of the roots of X z −1 on the positive real axis resulting from the alteration of the initial runoff value. This may violate the hydrologic stability condition. However, if this negative value is considered to have occurred due to a very high rate of surface evaporation in the middle of a complex rainfall event, the hydrologic stability conditions can be satisfied [21].
However, relating directly a high rate of evaporation occurring during a storm to the alteration of initial runoff seems quite delicate. Rather, this case seems to be more related to the method of calculating the direct runoff from a measured runoff. Such a negative coefficient in polynomial X z −1 can happen because a small amount of variation for runoff is acceptable when considering the common base flow separation process. Thus, it is recommended to derive the response functions with caution by referring to the SRS method proposed in this study.

Conclusions
This paper suggested a reliable method that can be used to estimate the S-curve in an objective manner from storms with unreliable rainfall data. This method depends only on the separation of an S-curve from the runoff data of a given basin similar to the root selection method. As an analysis input, cumulative runoff is used to promote less subjectivity in separating the roots of the S-curve and rainfall from the input approximated by power series expansion. This method was designed to outperform the root selection method, so that it was possible to extract the S-curve clearly by avoiding the unclear roots, and to solve the hydrologic stability problem. The applicability of the proposed approach was demonstrated for the Nenagh River basin. Based on the results, the following conclusions can be drawn: (1) The proposed method provided a clear and effective tool identifying the S-curve solely from the effective runoff of given storm data. (2) The roots of the S-curve polynomial described a clearer circular pattern that allowed root separation for the polynomial without any additional efforts. (3) The proposed method did not affect the root pattern of rainfall provided by the root selection method (4) The estimates of IUHs were in general agreement with those of previous studies. (5) The ordinates of the reproduced direct runoff were in good agreement with those of previous studies. (6) Hydrologic stability evaluation could be performed by comparing the roots of observed rainfall and abnormal roots on the Argand diagram. (7) The sensitivity of roots for rainfall was relatively larger than that of the S-curve polynomial with respect to a small variation of initial runoff ordinate.