SCS-CN parameter determination using rainfall-runoff data in heterogeneous watersheds – the two-CN system approach

The Soil Conservation Service Curve Number (SCS-CN) approach is widely used as a simple method for predicting direct runoff volume for a given rainfall event. The CN parameter values corresponding to various soil, land cover, and land management conditions can be selected from tables, but it is preferable to estimate the CN value from measured rainfall-runoff data if available. However, previous researchers indicated that the CN values calculated from measured rainfall-runoff data vary systematically with the rainfall depth. Hence, they suggested the determination of a single asymptotic CN value observed for very high rainfall depths to characterize the watersheds’ runoff response. In this paper, the hypothesis that the observed correlation between the calculated CN value and the rainfall depth in a watershed reflects the effect of soils and land cover spatial variability on its hydrologic response is being tested. Based on this hypothesis, the simplified concept of a two-CN heterogeneous system is introduced to model the observed CNrainfall variation by reducing the CN spatial variability into two classes. The behaviour of the CN-rainfall function produced by the simplified two-CN system is approached theoretically, it is analysed systematically, and it is found to be similar to the variation observed in natural watersheds. Synthetic data tests, natural watersheds examples, and detailed study of two natural experimental watersheds with known spatial heterogeneity characteristics were used to evaluate the method. The results indicate that the determination of CN values from rainfall runoff data using the proposed twoCN system approach provides reasonable accuracy and it over performs the previous methods based on the determination of a single asymptotic CN value. Although the suggested method increases the number of unknown parameters to three (instead of one), a clear physical reasoning for them is presented.


Introduction
Simple methods for predicting runoff from watersheds are particularly important in hydrologic engineering and hydrological modelling and they are used in many hydrologic applications, such as flood design and water balance calculation models (Abon et al., 2011;Steenhuis et al., 1995;van Dijk, 2010).The Soil Conservation Service Curve Number (SCS-CN) method was originally developed by the SCS (US Department of Agriculture), to predict direct runoff volumes for given rainfall events and it is documented in the National Engineering Handbook, Sect.4: Hydrology (NEH-4) (SCS, 1956(SCS, , 1964(SCS, , 1971(SCS, , 1985(SCS, , 1993(SCS, , 2004)).It soon became one of the most popular techniques among the engineers and the practitioners, because it is a simple but well-established method, it features easy to obtain and well-documented environmental inputs, and it accounts for many of the factors affecting runoff generation, incorporating them in a single CN parameter.In contrast, the main weaknesses reported in the literature are that the SCS-CN method does not consider the impact of rainfall intensity, it does not address the effects of spatial scale, it is highly sensitive to changes in values of its single parameter, CN, and it is ambiguous considering the effect of antecedent moisture conditions (Hawkins, 1993;McCuen, 2002;Michel et al., 2005;Ponce and Hawkins, 1996).
The SCS-CN method was soon adopted for various regions, land uses and climate conditions (Elhakeem and Papanicolaou, 2009;King and Balogh, 2008;Mishra and Singh, 1999;Romero et al., 2007).It was also evolved well beyond its original scope and it became an integral part of continuous simulation models (e.g.Adornado and Yoshida, 2010;Holman et al., 2003;Mishra and Singh, 2004;Moretti and Montanari, 2008;Soulis and Dercas, 2007).Many studies aiming at finding a theoretical basis for the method, facilitating its use in regions and for climate conditions not previously evaluated, and supporting its further improvement, were carried out as well (Hjelmfelt, 1991;Tramblay et al., 2010;Yu, 1998).
However, in spite of its widespread use, there is not an agreed methodology to estimate the CN parameter values from measured rainfall runoff data.Such a method would be important for two main purposes: (a) it would allow the determination of the CN parameter values from measured rainfall runoff data of local or nearby similar watersheds when suitable data were available and (b) it would facilitate studies aiming at the extension of the SCS-CN method documentation for different, soil, land use, and climate conditions.Though, the main difficulty is that the CN values calculated from measured rainfall runoff data actually vary significantly from storm to storm on any watershed.This effect posed in doubt the adequacy of curve number model itself to predict runoff.Antecedent Moisture Condition (AMC) was initially assumed to be the primary cause of storm to storm variation.However, this effect is of questionable origin and it is not recommended for use anymore (Hjelmfelt et al., 2001;McCuen, 2002;Ponce and Hawkins, 1996).In the latest version of the NEH-4 the reference to AMC was revised as follows.Variability is incorporated by considering the CN as a random variable and the AMC-I and AMC-III categories as bounds of the distribution.The expressions of AMC-I and AMC-III were considered as measures of dispersion around the constant tendency (AMC-II) (Hjelmfelt et al., 2001).Ponce and Hawkins (1996) reported as possible sources of this variability the effect of the temporal and spatial variability of storm and watershed properties, the quality of the measured data, and the effect of antecedent rainfall and associated soil moisture.Soulis et al. (2009) and Steenhuis et al. (1995) also noted that the variation of CN value, according to AMC category alone, cannot justify the observed CN values variability in every case.Hawkins (1993) in his study on the asymptotic determination of runoff curve numbers from measured runoff analysing a significantly large number of watersheds, where CNs are calculated from real rainfall-runoff data, concluded that a secondary systematic correlation almost always emerges in watersheds between the calculated CN value and the rainfall depth.In most of the watersheds, these calculated CNs approach a constant value with increasing rainfall depth that is assumed to characterize the watershed.The three different behaviours that have been observed are described as follows: the most common scenario is that at small rainfall depths correspond larger values of calculated CNs, which decline progressively with increasing storm size, approaching a stable near constant asymptotic CN value with increasingly larger storms.This behaviour appears most frequently and it is characterized as "standard".An example of this pattern is given in Fig. 1.Hawkins (1993) suggests the identification of a single asymptotic CN value observed for very large  Hawkins (1993) for the "standard" (Coweeta watershed #2, North Carolina) and the "complacent" (West Donaldson Creek, Oregon) behaviour watersheds.storm sizes to characterize such watersheds.In less common cases of watersheds the observed CN declines steadily with increasing rainfall with no appreciable tendency to approach a constant value ("complacent" behaviour, Fig. 1).According to Hawkins (1993), an asymptotic CN cannot be safely determined from data for this behaviour.In the last case, concerning also a small number of watersheds, the calculated CNs have an apparently constant value for all rainfall depths except very low rainfall depths where CN increases suddenly ("violent" behaviour).
Additional examples of watersheds featuring similar behaviours are presented by Hjelmfeld et al. (2001).Bonta (1997) proposed an improvement to the Hawkins (1993) method for the asymptotic determination of CNs from measured data in "violent" and "standard" watersheds using derived distributions.
All previously developed methodologies for estimating CNs from measured data focus mainly on the determination of a single asymptotic CN value characterizing the watershed hydrologic response for high rainfall depths.The observed deviations from the asymptotic behaviour for lower rainfall depths are not essentially taken into consideration and are rather attributed to various sources of temporal variability.For this reason, the resulting CN values fail to describe the watershed response in small and medium rainfall events, limiting the applicability of the method to its original scope, namely the estimation of peak runoff values.Furthermore, the above methods fail to determine a final CN value in "complacent" watersheds.The CN varies as a function of the soil infiltration capacity and the land cover of the watershed, which are two essentially time invariant factors.
Various sources of temporal variability, such as the effect of spatio-temporal rainfall intensity variability, the effect of antecedent rainfall, etc., make CN to be considered as a random variable with bounds of distribution AMC-I and AMC-III.The SCS-CN method was originally developed as a lumped model and up to this date it is still primarily used as a lumped model.In natural watersheds, however, spatial variability (at lower or higher level) with regard to the soil-cover complex is inevitable (such spatial heterogeneity in the watershed could be considered temporally invariant).
In this paper, a novel hypothesis is proposed suggesting that the intrinsic correlation between calculated CN value and rainfall depth observed in watersheds corresponding to the "standard" and "complacent" cases is essentially the natural consequence of the presence of soils and land cover spatial variability along the watersheds.It is shown that the presence of spatial variability (at low or high level) in the watersheds produces a progressive decrease in the calculated CNs as the storm size decreases and for excessively large storm sizes the CN tends to stabilize in an asymptotic CN value.The proposed hypothesis is approached theoretically, it is analysed systematically using synthetic data, it is studied in two natural experimental watersheds with known spatial heterogeneity characteristics and it is evaluated using natural watersheds examples.The results of the analysis provide evidence that the spatial variability of the watershed can influence the CN determination procedure from measured rainfall-runoff data and that the estimation of more than one CN values is needed in order to describe the spatial variability of the watershed and to facilitate the determination procedure.Based on the above hypothesis, the simplified concept of an equivalent two-CN heterogeneous system is introduced to model the CN vs. rainfall depth variation.This new evolution takes into consideration the soil-cover complex spatial variation in the estimation of CN values from measured rainfall-runoff data, in order to extend the applicability of the SCS-CN method for a wider range of rainfall depths and to provide improved simulations in heterogeneous watersheds.

SCS-CN method
The SCS-CN method is based on the following basic form calculating runoff from rainfall depth, where P is the total rainfall, I a is the initial abstraction, Q is the direct runoff and S is the potential maximum retention.
Based on a second assumption, that the amount of initial abstraction is a fraction of the potential maximum retention Eq. ( 1) becomes The potential retention S is expressed in terms of the dimensionless curve number (CN) through the relationship taking values from 0, when S → ∞, to 100, when S = 0.This definition was originally applied to the English metric system (with S in inches).In the SI units (with S in mm) the following definition should be used: The determination of all the NEH-4 SCS-CN values commonly used in hydrologic practice, assume the initial abstraction rate to be set to the constant value, λ = 0.2, in order that S (or its transformation CN) remains the only free unknown parameter of the method.Recently, Woodward et al. (2003) analysing event rainfall-runoff data from several hundred plots recommended using λ = 0.05.The CN values corresponding to the various soil types, land cover and land management conditions can be selected from the NEH-4 tables.However, it is preferable to estimate the CN value from recorded rainfall-runoff data from local or nearby similar watersheds.When rainfall-runoff data are available for a watershed, P and Q pairs are used directly to determine the potential retention S characterizing the watershed (Chen, 1982) Combining Eq. (4b) with Eq. ( 5), CN value can be directly calculated from rainfall-runoff data CN = 25 400 2.2 Runoff prediction errors related to the use of single composite CN values Grove et al. (1998) in their study investigated the effect of using single composite CN values (i.e. the area-weighted average of the CN values in the watershed) instead of weighted runoff estimates, indicating that significant errors in runoff estimates can occur when composited rather than distributed CNs are used.Lantz and Hawkins (2001) also discussed the possible errors caused by the use of a single composite CN value.
The main reason for the errors produced using the composite CN value instead of weighted-Q is the non-linear form of the SCS-CN formula.As an example, the case of a virtual watershed divided into two equal sub-areas characterized by different CN values is illustrated in Fig. 2. In this figure the relative percentage error of the runoff predictions using a single composite CN value is plotted against the range of CN variation, for various total rainfall depths and for various average CN values.The above figure clearly illustrates that the percentage error increases as the range of CN variation increases and decreases as the average CN value and the rainfall depth increase.It is also clearly shown that for low rainfall depths significant errors are observed, even for small CN variation ranges.These results are in agreement with the results of Grove et al. (1998).

The two-CN heterogeneous system
In order to investigate the consequence of spatial variability on the CN vs. P relationship in a watershed, in a first stage of the analysis it is assumed the simplified scheme, according to which the entire area of the watershed under consideration is composed from relatively homogeneous sub-areas.Each sub-area is assigned a CN value obtained from a specific set of two CN values CN a and CN b with CN a > CN b .If a denotes the area fraction of the watershed with CN = CN a , then (1-a) is the area fraction of the watershed with CN = CN b .It seems obvious that CN must be taken constant for a relatively homogeneous soil-cover complex.Various temporal effects such as the effect of the spatiotemporal variability of given storm, the effect of storm intensity, the effect of antecedent rainfall and others are considered as random effects on the CN calculation.
Traditionally the runoff equation for a heterogeneous watershed is described by using a single composite value of the different CN-areas, this being an area-weighted CN value.However, runoff is more accurately estimated using individually calculated weighted runoff for the array of different subareas as it was shown in the previous section.Therefore, the runoff, Q responded to the causative rainfall event, P generated by the two-CN system is described by the following equation, for λS a ≤ P < λS b (7b) where S a and S b are the potential maximum retention values corresponding to the two homogeneous sub-areas characterized by the CN α and CN b values respectively, and λ is a constant value (usually λ = 0.2 or λ = 0.05).S a and S b are calculated from the corresponding CN values using Eq. ( 4b).Following, it will be pointed out that such a two-CN heterogeneous system is characterized by a secondary relationship that always emerges between calculated CN and rainfall depth, P .The particular behaviour of this relationship will be analysed in detail as well.
It is considered that for various rainfall events of depth P , realized on the two-CN heterogeneous system, the corresponding "actual" observed runoff, Q, is obtained by Eq. (7a, b, c).Then the CN for this system can be calculated by Eq. ( 6) containing only P and Q; thus any "realized" P -Q data pair can be used to calculate what should be the CN for that particular rainfall-runoff event in the heterogeneous system.

Large-P behaviour -Asymptotic CN
Equation (7c) can be standardized by using the reduced variables (P /S a ), and (P /S b ), (S a < S b ).The resulting relationship becomes: while using the auxiliary variables X 1 = P /S a + (1 − λ) and For asymptotic large values of P and consequently asymptotic large values of X 1 and X 2 , the corresponding value of Q ∞ approaches asymptotically the value or equivalently By following a similar procedure assuming a perfectly uniform watershed characterized by a single CN-value (or its simple transformed S), the value of Q ∞ for large values of P approaches asymptotically By putting S ∞ = aS a +(1−a)S b in Eq. ( 11) the two-CN heterogeneous system behave asymptotically for large P values as a single CN value system with equivalent potential retention S ∞ and equivalent CN value Only for large values of P the heterogeneous system can be characterized by a single asymptotic CN value that could be obtained using the specific "composite" CN value (Eq.13).However, even in this case this asymptotic value does not characterize a single specific soil but it is the superposition of different complexes.Systematic analysis indicates that the value of CN ∞ given by Eq. ( 13) is sufficiently close to the usual composite CN value Further analysis based on systematic generation of Q-P synthetic data for various combinations of a, CN a and CN b input parameters characterizing the two-CN system indicates that CN approaches the asymptotic value given by Eq. ( 13) for unrealistic, extremely large values of P , P > 3000 mm.Alternatively the CN approaches the composite asymptotic value given by Eq. ( 14) for more realistic large values of P .Note that the composite value given by Eq. ( 14) is traditionally used to characterize an heterogeneous system by a single-CN value.

Low-P behaviour -Envelope curve
For a two-CN system, as P decreases the calculated values of CN increase, as illustrated in Fig. 3.For some threshold value of P , the CN value becomes maximum equal to the larger CNcategory, CN a , whereas for any smaller P < P o value the CN is not defined since it will give no runoff Q = 0. Indeed for P -Q pairs generated by Eq. (7a, b, c), when P decreases approaching asymptotically the value of P o , then Q → 0 therefore the asymptotic threshold value of S,S o , calculated by Eq. ( 5 and 3)

Illustration of the two-CN heterogeneous system behaviour
In order to illustrate the behaviour of the secondary relationship between the calculated CN and the rainfall depth, P in the above described two-CN heterogeneous system, "actual" observed runoff values, Q, were obtained by Eq. (7a, b, c) for various rainfall depths P , by varying systematically the a, CN a , and CN b parameters.Then the corresponding CN values for this system were calculated by Eq. ( 6) and a series of CN-P curves were produced.It must be noticed that hereafter, the standard case of λ=0.2 is examined.However, the following analysis is also valid for other λ values, as well.In Fig. 3 the calculated CN values for the various values of a, CN a , and CN b parameters are plotted against the rainfall depth P .In this figure, a significant variation of the estimated CN values for various rainfall depths can be observed.The variation increases as the difference between CN a and that the factors associated with significant variation of the estimated CN values for various rainfall depths, are also associated with significant errors when runoff estimations are made using composited rather than distributed CNs, as it was shown in Sect.2.2.This observation provides a strong indication that the observed correlation between the calculated CN values and the rainfall depth should be associated with the presence of soil-cover complex spatial variability in the watershed.
In Fig. 3 can be also observed that the shapes of the CN-P curves produced by the two-CN heterogeneous system are quite similar with the shapes of the "standard" and "complacent" watersheds correlation curves presented by Hawkins (1993).When Q-P data are available, the two-CN system can be viewed as a fitting model to the transformed CN-P data with free parameters a, CN a , and CN b (the equations of the two-CN system fitting model that can be used in a non-linear least squared procedure, are given in the Appendix A).Thus, in order to highlight further the similarity observed in Fig. 3, the two-CN hypothetical watershed curves were fitted to the CN-P curves presented by Hawkins (1993) as examples of the "standard" (Coweeta watershed #2, North Carolina) and of the "complacent" (West Donaldson Creek, Oregon) behaviour, by adjusting the values of the a, CN a , and CN b parameters.As it can be clearly seen in Fig. 1, the CN-P curves are fitted very well by the two-CN system model in both cases.These results provide further evidence that the spatial variability of the watershed can influence the CN determination procedure.In this case the estimation of more than one CN values is needed in order to describe the spatial variability of the watershed and to facilitate the determination procedure.

Generalization
Although the previous analysis is initially restricted for two-CN idealized watershed examples, generally, in natural watersheds could appear more than two CN value categories.However, every added CN category requires the determination of two more parameters (the corresponding CN value and the area it covers), giving rise to the overparameterization problem.Therefore, in a second stage it is investigated if a heterogeneous watershed characterized by three different CN values can be approached with sufficient accuracy using two CN value categories.For this purpose, synthetic runoff data for 21 hypothetical watersheds that are characterized by three CN value categories have been created.The selected examples cover a wide variety of possible cases including watersheds with various ranges of CN variation and watersheds dominated by the lower, the medium or the higher CN value (Table 1).The synthetic runoff data were calculated as the weighted average of the runoff values resulted by the SCS-CN method for the three CN values characterizing each hypothetical watershed, for rainfall depths ranging from 0 to 300 mm.Then, the corresponding a, CN a , and CN b parameters were determined by fitting the two-CN system model to the synthetic CN-P data.As it can be observed in Fig. 4 and at the results presented in Table 1, the synthetic CN-P curves are fitted very well by the two-CN system model in all the examples examined.
In Fig. 5, the synthetic runoff data for six characteristic examples of hypothetical watersheds comprising three CN value categories are plotted in comparison to the runoff predictions of the SCS-CN method using the single composite CN value, the single asymptotic CN value according to Hawkins (1993), the best fitted single CN value, and the CN values determined with the two-CN system model.In this figure it can be observed that the SCS-CN method using a single CN value category can provide adequate results only in the case that one CN category dominates runoff production in the watershed (e.g. in case 3).In all other cases the use of two CN categories provides much better results.

Case studies
The validity of the above analysis in natural watersheds is investigated in two representative examples, the Little River N Experimental Watershed and the Lykorrema Experimental Watershed.These watersheds were selected because they have been presented in the literature as examples of the "standard" and the "complacent" behaviour respectively, and for both of them, detailed geographical data were available (Hjelmfelt et al., 2001;Soulis et al., 2009).

Little river subwatershed N
The Little River Experimental Watershed (LREW) (Fig. 6a), is one of twelve national benchmark watersheds participating in the Conservation Effects Assessment Project-Watershed Assessment Studies (CEAP-WAS) (Bosch et al., 2007a).It is located near Tifton, Georgia, in the western headwaters area of the Suwannee River Basin, centred at approximately 31.61 • N and 83.66 • W. The Suwannee River Basin is completely contained in the Gulf-Atlantic Coastal Plain physiographic region, which is characterized by low topographic relief (Sheridan, 1997).Climate in this region is characterized as humid subtropical with an average annual precipitation of about 1167 mm.Hydrology, climate and geographical data at LREW have been monitored by the ARS Southeast Watershed Research Laboratory (SEWRL) since the 1960s (Bosch and Sheridan, 2007;Bosch et al., 2007a, b;Sullivan and Batten, 2007;Sullivan et al., 2007).
The 15.7 km 2 Little River subwatershed N (LRN) (Fig. 6a) was presented by Hjelmfelt et al. (2001) as a characteristic example of a "standard" watershed.The main soil series in LRN are Tifton loamy sand (48 %), Alapaha loamy sand (16 %), and Kinston and Osier fine sandy loam (6 %).The agricultural lands are mostly covered by Tifton series soils having moderate infiltration rates (hydrologic soil group B), while the areas around the stream and wetland areas are covered by Alapaha and Kinston-Osier soils (hydrologic soil group D).As it is reported by Lowrance et al. (1984), row crops, pasture, and riparian forests cover approximately 41, 13, and 30 % of LRN, respectively, while the remaining 16 % includes roads, residences, fallow land, and other land cover types.

1
Figure 5. Synthetic runoff data in comparison to the runoff predictions of the SCS-CN method 2 using the single composite CN value, the single asymptotic CN value according to Hawkins 3 (1993), the best fitted CN value, and the proposed two-CN system model, for six 4 characteristic cases as described in Table 1. 5 6 Fig. 5. Synthetic runoff data in comparison to the runoff predictions of the SCS-CN method using the single composite CN value, the single asymptotic CN value according to Hawkins (1993), the best fitted CN value, and the proposed two-CN system model, for six characteristic cases as described in Table 1.

Lykorrema, Penteli
The small scale experimental watershed of Lykorrema stream (15.2 km 2 ), situated in the east side of Penteli Mountain, Attica, Greece, centred at approximately 38.02 • N and 23.55 • E (Fig. 6b).The watershed is divided in two sub-watersheds.The Upper Lykorrema watershed (7.84 km 2 ) and the Lower Lykorrema watershed (7.36 km 2 ).The Upper and Lower Lykorrema experimental watersheds are operated from the Agricultural University of Athens, Greece and the National Technical University of Athens, Greece, respectively.
The region is characterized by a Mediterranean semi-arid climate with mild, wet winters and hot, dry summers.The yearly average precipitation value is 595 mm.The watershed presents a relatively sharp relief, with elevations ranging between 146 m and 950 m.The watershed is dominated by sandy loam soils with high infiltration rates (hydrologic soil group A, 64 %) and a smaller part is covered by sandy clay loam soils presenting relatively high infiltration rates (hydrologic soil group B, 29 %).The dominant land cover type in the watershed is pasture with a few scattered tufts of trees (93 %).The remaining 7 % includes roads, residences, bare rock and other land cover types.Detailed description of the hydrology, climate and physiography of Lykorrema experimental watershed and of the available geographical and hydro-meteorological databases are provided by Baltas et al. (2007), Soulis (2009), andSoulis et al. (2009).

Identification of spatial distribution of CN along watersheds from measured data using the two-CN system
In a first attempt a simplified identification procedure is proposed for spatially distribute along the watershed the two-CN categories using the measured P -Q data.The simplified procedure includes the following steps: 1.The measured P and Q values are sorted separately and then realigned on a rank order basis to form P -Q pairs of equal return period following the frequency matching technique (Hawkins, 1993;Hjelmfelt et al., 1980Hjelmfelt et al., , 2001)).Then the measured P -Q data are transformed in the equivalent P -CN data using Eq. ( 6).
2. The two-CN system model (Eq.A1, A2, A3) is fitted to the transformed CN-P measured data curve yielding a first set of best estimates of parameters a (•) , CN 3. The watershed is divided in a set of n relatively uniform subareas with constant soil-cover complex.The subareas are clearly spatially identified along the watershed.For each subarea characterized by a specific soilcover complex an initial approximate CN (table) value is attributed based on the NEH-4 tables.The areas of all subareas characterized by each specific CN (table) value are also determined.The m different CN (table)   In order to more closely describe the real conditions of natural watersheds it could be proposed using as free parameters three or even four CN categories to be spatially distributed along the watershed, however such a procedure has an additional risk to appear non-convergence and non-unique solution problems when the inverse solution procedure is applied.Hjelmfelt et al. (2001), using the measured P -Q data obtained the transformed CN-P measured data curve for the LRN watershed, in a similar way to the first step of the proposed methodology (Fig. 7).Applying the second step of the proposed methodology, the two-CN system model (Eq.A1, A2, A3) was fitted to the above mentioned CN-P measured data curve presented by Hjelmfelt et al. (2001) (Fig. 7) yielding the best estimates of the three fitting parameters:  , 2007).Each subarea characterized different CN (table) (as selected from the NEH-4 tables) was spatially identified along the watershed.Fig- 8a presents the categories spatial distribution along the watershed.Then the cumulative fraction area for each CN (table) category was determined.The cumulative area fractions distribution curve for the various approximate CN values is presented in Fig. 9.The single composite CN value was also determined equal to CN (table ) = 71.From the cumulative area fraction distribution curve (Fig. 9) the value of A = 0.137 was selected as the closest value to the value of a (•)  = 0.151 obtained using the P -Q measured data, as it is described in the fourth step of the proposed methodology.Then, the two-CN system model For comparison reasons, the two composite CN values corresponding to the area fractions of the watershed equal to a and 1-a were also calculated according to the tables and the methodology provided in NEH-4, and based on the available soil and land cover data.The resulted CN values were equal to 83 and 69 respectively.These values are comparable to the best estimates of CN a , and CN b parameters' values obtained from the measured P -Q data.The LRN watershed is clearly a heterogeneous watershed with CN varying between 100 and 55 according to the tables and the methodology provided in NEH-4.The above results provide strong indica- tions that the observed correlation between the CN values and the rainfall depths presented in Fig. 7 is essentially related to the spatial variability of the watershed.Additionally, it can be noticed that the estimation of two CN values can sufficiently describe the spatial variability of the watershed.

Lykorrema, Penteli
Following the first step of the proposed methodology, the measured Q-P data presented by Soulis et al. (2009), were sorted separately and then realigned on a rank order basis to form P -Q pairs of equal return period and then were transformed in the equivalent P -CN data curve using Eq. ( 6) (Fig. 10).At the next step, the two-CN system model (Eq.A1, A2, A3) was fitted to the produced CN-P curve (Fig. 10) yielding the best estimates of the three fitting parameters: a (•)  = 0.068, CN Then, in the same way as in the previous case study, the approximate values of curve numbers and their spatial distribution along the watershed were initially estimated by selecting them according to the tables and the methodology provided in NEH-4, based on the available soil and land cover data (Soulis, 2009;Soulis et al., 2009).Each subarea characterized by different CN (table) (as selected from the NEH-4 tables) was spatially identified along the watershed.Figure 11a presents the CN (table) categories spatial distribution along the watershed.Then the cumulative fraction area for each CN (table) category was determined.The cumulative area fractions distribution curve for the various approximate CN values is presented in Fig. 12.The single composite CN values were also determined equal to CN (table ) = 51 and CN (table ) = 55 for the Upper Lykorrema watershed and for the entire watershed, respectively.
From the cumulative area fraction distribution curve (Fig. 12) the values of A = 0.052 and A = 0.075 are selected as the closest values to the corresponding a (•) values for the Upper and the Entire Lykorrema watershed respectively, as it  is described in the fourth step of the proposed methodology.Then, the two-CN system model (Eq.A1, A2, A3) was once again fitted to the transformed CN-P measured data leading to the parameters CN  The Lykorrema watershed is also a heterogeneous watershed with CN varying between 100 and 45 according to the tables and the methodology provided in NEH-4.Furthermore, it can be observed that the area fractions of the watershed corresponding to the higher best estimate CN value (CN a ) are comparable to the area fractions of the watersheds covered with impervious or nearly impervious surfaces (e.g.roads, buildings, bare rock and stream beds), which are equal to 0.051 and 0.075 for the Upper and the Entire Lykorrema watershed respectively.
In an analogous way as in the LRN case study, the obtained results provide strong indications that the observed correlation between the CN values and the rainfall depths presented in Fig. 10 is essentially related to the spatial variability of the watersheds and that the estimation of two CN values can sufficiently describe the spatial variability in both cases.

Discussion
In this work it is assumed that the specific behaviour in watersheds, according to which CN systematically varies with rainfall size (Hawkins, 1979(Hawkins, , 1993)), reflects the effect of the inevitable presence of spatial variability of the soil -cover complexes of watersheds.Since this characteristic of the watershed can be considered invariant in time, therefore in all statistical studies concerning the variation of CN in a watershed, the produced effect of heterogeneity (e.g. the CN-P relationship) should be included as a deterministic part of the analysis.Other, temporally variant, causes of variability (e.g.rainfall intensity and duration, soil moisture conditions, cover density, stage of growth, and temperature) can explain the remaining scatter around the main rainfall-CN correlation curve.The concept of a simplified idealized heterogeneous system composed by two different CN values is introduced.The behaviour of the CN-P function produced by such a system was analysed systematically and it was found similar to the CN-P variation observed in natural watersheds (Fig. 1, 7, 10).Measured P -Q data can be used to identify the two different CN values and the corresponding area fractions of the simplified two-CN system.Then the initial threshold value CN o and the asymptotic large S value of CN ∞ are also obtained and the characteristics of the CN(P ) as well as Q(P ) functions are determined.
The proposed method is advantageous over previous methods suggesting the determination of a single asymptotic CN ∞ value to characterize the watershed runoff behaviour as it permits the accurate prediction of runoff for a wider range of rainfall depths (including low and medium rainfall depths) and not for excessively large storms only (it must be noticed that the asymptotic CN ∞ value is essentially observed for excessively large P > 3000 mm).Therefore, the proposed method can be also used in continuous hydrological models.
To illustrate if the proposed method of CN determination in heterogeneous watersheds provides improved runoff predictions over a wider range of rainfall depths than the traditional method that is based on the determination of a single asymptotic CN value, in Fig. 13, the measured runoff is plotted against the rainfall depth for two "standard" and two "complacent" watersheds' examples presented in the literature.At the same figure the runoff predictions of the SCS-CN method using the CN values obtained by the proposed CN determination methodology assuming a two-CN system as well as the runoff predictions of the SCS-CN method based on the determination of a single asymptotic CN value proposed by Hawkins (1993), are also plotted.
In Fig. 13a can be observed that the proposed methodology over performs the previous original CN determina-tion method even if the "Coweeta" watershed was selected as a characteristic example of the "standard" behaviour in the study of Hawkins (1993) concerning the asymptotic CN determination method.Furthermore, significant errors are observed for low and medium runoff predictions (for P < 100 mm) when the traditional asymptotic method is used.Similar observations can be made in Fig. 13b for the LRN watershed, which was also presented as a characteristic example of the "standard" behaviour by Hjelmfelt et al. (2001) even if the difference in this case is small.
The advantages of the proposed method are more evident in Fig. 13c and d, where two characteristic examples of "complacent" behaviour watersheds presented by Hawkins (1993) and Soulis et al. (2009), respectively, are demonstrated.As it can be clearly seen, satisfactory runoff predictions can be obtained using the CN values determined by the proposed methodology.In contrast, the CN values determined with the asymptotic method completely fail to predict runoff.It must be noticed that according to Hawkins (1993) and Hjelmfelt et al. (2001), an asymptotic CN cannot be determined from data for "complacent" watersheds.For this reason, the runoff predictions obtained based on the best fitted single CN values were also plotted in Fig. 13c and d.It can be seen once again that the runoff predictions obtained are very poor in both cases as well.These results are in agreement with the results of the detailed analysis based on synthetic data (Fig. 5) presented in the Sect.2.4.
In previous analysis it is demonstrated that the presence of heterogeneity produces CN-P correlations that stabilize to a steady state regime (asymptotic value) for large values of P .Therefore the "complacent" behaviour could be considered as a specific case, in which the available range of rainfall measurements dataset is restricted in such a way that the steady state regime is not yet established and thus an asymptotic CN value cannot be determined from this dataset.
In Figs.7b and 11b the spatial distribution of the estimated CN values in the two case studies is presented.In these figures, the association of the a, CN a , and CN b parameters to the actual characteristics of the watersheds is highlighted.The ability of the proposed methodology to provide information on the spatial distribution of the estimated CN values is also demonstrated.

Conclusions
Considering the theoretical analysis, the systematic analysis using synthetic data and the detailed case studies it can be concluded that the observed correlation between the calculated CN value and the rainfall depth in a watershed can be attributed to the soils and land cover spatial variability of the watershed and that the proposed two-CN system can sufficiently describe the CN-rainfall variation observed in natural watersheds.The results of the synthetic data analysis (Fig. 5) and the results of the real watersheds examples (Fig. 13  indicate that the SCS-CN method using the CN values obtained by the proposed CN determination methodology provides superior runoff predictions in most cases and extends the applicability of the original SCS-CN method for a wider range of rainfall depths in heterogeneous watersheds.Furthermore, the proposed methodology allows the CN determination in "complacent" watersheds.Although the sug-gested method increases the number of unknown parameters to three, a clear physical reasoning for them is presented.A simplified procedure to identify the spatial distribution of the two different CN values along the watersheds (Fig. 8b, 11b) is also presented.Taking into consideration this additional capability, i.e. to provide information on CN values spatial distribution and thus spatially 1014 K. X. Soulis and J. D. Valiantzas: The two-CN system approach distributed runoff estimations, the proposed method can be used in other environmental applications e.g.water quality studies or estimation of erosion hazard.
The next step of this approach could be the validation of the proposed methodology to additional experimental watersheds with known characteristics.This is needed for a more definitive validation, and might lead to some adaptations of the proposed conceptual model for explaining the intrinsic correlation of CN-P data.However, despite these reservations, it is quite interesting that the observed CN-P correlation in watersheds can be the effect of an intrinsic characteristic of the natural watersheds, which is the spatial heterogeneity.This observation may facilitate future studies aiming at the extension of the SCS-CN method documentation for different regions and different soil, land use, and climate conditions.

1Figure 2 . 4 Fig. 2 .
Figure 2. Relative percentage error against the range of CN variation, for various total rainfall 2 depths and for various average CN values.3 4

Figure 3 .Fig. 3 .
Figure 3. Calculated CN values against rainfall depth for various values of the a, CNa, and 2 CNb parameters.3 4 Fig. 3. Calculated CN values against rainfall depth for various values of the a, CN a , and CN b parameters.
) is S o = P o /λ.Since S a is also given by S a = P o /λ, therefore the threshold value of CN, CN o = CN a .The values of threshold maximum curve number, CN o as function of P o is threshold CN o (P o ) curve is an envelope curve that could be interpreted as the intrinsic CN(P ) variation for a two-CN system with asymptotic characteristics CN a → 100, CN b → 0, and a → 0. It is the curve defining the position of max CN o = CN a value at the threshold P = P o =λS a (seeFigs.1

Figure 4 .Fig. 4 .
Figure 4. Two-CN system model curves fitted to the synthetic rainfall-CN data created for the 221 examples of hypothetical watersheds that are characterized by three CN value categories as 3 described in Table1.4 5

1Figure 7 .Fig. 7 .
Figure 7. Two-CN system model fitted to the data presented by Hjelmfelt et al. (2001) for the 2 LRN watershed.3 4 •) a = 86, and CN (•) b = 63.At the next step, approximate values of numbers and their spatial distribution along the watershed were initially estimated them according to the tables and the methodology provided in NEH-4, based on the soil and land cover data in the LREW geographical database et al.

Fig. 9 .
Figure 9. LRN watershed cumulative area fraction distribution curve. 2 3 for the Upper and the Entire Lykorrema watershed respectively.
Figure 10.Two-CN system model fitted to the rainfall-CN data presented by Soulis et al. 2 (2009) for the (a) Upper and (b) Entire Lykorrema watersheds.3 4

Fig. 10 .
Fig. 10.Two-CN system model fitted to the rainfall-CN data presented by Soulis et al. (2009) for the (a) Upper and (b) Entire Lykorrema watersheds.

=
100 and CN (distr) b = 40 for the Upper and the Entire Lykorrema watershed respectively (step 5).The resulted spatial distribution of the estimated CN (distr) a and CN (distr) b parameters is presented in Fig. 11b.

Figure 13 .
Figure 13.Measured runoff against the rainfall depth in comparison to the runoff predictions 2 of the various CN value determination methods for two "Standard" (a, b) and two 3 "Complacent" (c, d) watersheds' examples.4

Fig. 13 .
Fig. 13.Measured runoff against the rainfall depth in comparison to the runoff predictions of the various CN value determination methods for two "standard" (a, b) and two "complacent" (c, d) watersheds' examples.

Table 1 .
Characteristics of 21 examples of hypothetical watersheds that are characterized by three CN value categories and best fitted values of the a, CN a , and CN b parameters.
b parameters value increases and decreases as the rainfall depth and the weighted CN value increase.It is clearly shown as well that for very high weighted CN values, the estimated CN value is almost invariable.It can be observed ,A 2 , ..., A m cumulative fractions area.4.The A( i=1,m ) values are compared to the best estimate fraction parameter a(•)and the A i value closer to the a (•) (e.g.A j ) is selected.5.The two-CN system model is once again fitted to the CN-P measured data curve by fixing the parameter a = A j and treating CN a and CN b as free parameters leading to CN