Retention and Curve Number Variability in a Small Agricultural Catchment: The Probabilistic Approach

Department of Land and Water Resources Management, Slovak University of Technology,Radlinskeho Street 11, Bratislava 813 68, Slovak Republik; E-Mail: silvia.kohnova@stuba.sk* Author to whom correspondence should be addressed; E-Mail: rmrutkow@cyf-kr.edu.pl;Tel.: +48-12-6624544.Received:31 October 2013 / in revised form: 16 April 2014 / Accepted: 21 April 2014 /Published: 29 April 2014Abstract: The variability of the curve number (CN) and the retention parameter (S) ofthe Soil Conservation Service (SCS)-CN method in a small agricultural, lowland watershed(23.4 km

Keywords: curve number; retention parameter; theoretical distribution; antecedent runoff conditions

Introduction
For the estimation of design floods and volumes, often, event-based hydrological models are applied. As inputs, they often use the components of rainfall, infiltration, soil properties, land use, etc. The variability of these components can cause the instability of model parameters. This relates especially to small catchments, which are sensitive to environmental changes [1], observed as, e.g., land use and land cover changes, as well as changes in short-duration intense rainfalls, which considerably influence flood frequency curves for small catchments to a greater degree than for large catchments [2,3].
The rainfall-runoff curve number method was introduced in the 1950s by the United States Department of Agriculture-Soil Conservation Service (USDA-SCS, now NRCS (Natural Resources Conservation Service)) to estimate the runoff depth from small catchments. It is included also in many complex models, e.g., for the estimation of the rainstorm-generated sediment yield [4] or SWAT [5], which predicts the impact of land use on water, sediment and chemicals.
The method allows the relationship among the rainfall depth, P , the retention parameter, S, the curve number CN , the initial abstraction, I a , the initial abstraction ratio λ = I a /S and the direct runoff, H, to be collected in deterministic equations. In the model, the retention parameter, S (maximum potential catchment storage), expresses the response of the watershed to the rainfall event. The basic equations for the SCS-curve number (CN) method are given by [6]: where H, P and S are in millimeters. The parameters, CN and S, characterize the reaction of the catchment to the rainfall event. The method assumes them to be constant in each watershed (or in a part of it), being functions of geomorphologic and physiographic properties and land use. The tabulated value, CN tab , is determined on the basis of the National Engineering Handbook, Section 4 (NEH4) tables [7], due to soil type and land use. The primary source reference describing the CN method was the NEH4, later followed with a number of updates. The current, comprehensive overview of investigations, their applications, calibration and practice are the ASCEitems [6,7]. Many studies were done, e.g., by Hawkins [8][9][10], who investigated CN with declining rainfall depth, dealt with the CN infiltration rate equation and the infiltration rate data and formulated CN for antecedent moisture conditions. Several modifications of the method were conducted, e.g., in [11] the inclusion of the static portion of infiltration and the antecedent moisture was proposed and, in [12], the storm intensity incorporated. Several authors investigated the method of SCS-CN and tested it for various hydrological conditions of events. The testing of the method showed its deficiencies, pointed out by several authors. Therefore, in 1990, a working group was created, whose task it was to explore and improve the SCS-CN method. Their contributions were studies of seasonal and regional variations of the runoff curve number and a review of the ratio, I a /S.
In practice, a variability of S and CN is observed, as a result of seasonal and storm morphology changes. To describe this variability and assess errors in the estimation of CN , many statistical studies were performed; Hjelmfelt [13] introduced error bands in the studies on the variable nature of CN ; McCuen [14] dealt with confidence intervals of 100-CN , studied a peak discharge model with the SCS method and revealed its great sensitivity to CN changes. Later, Schneider and McCuen [15] formulated guidelines for CN calibration. Many statistical methods in practical applications of the model were introduced, e.g., by Tedela et al. [16] and Bondelid et al. [17], who revealed the decrease of the sensitivity of the SCS-CN model to CN variation with increasing rainfall depth, and by Epps et al. [18], who analyzed the influence of seasonal trends in water table elevation in coastal watersheds.
The asymptotic approach recognizes the drift of the CN s with increasing rainfall depth [19]. This study is based on the assumption that the return period of runoff depth equals the return period of rainfall depth, i.e., that the frequencies of rainfalls and runoffs are associated. Three typical responses of watersheds were specified by Hawkins: standard, complacent and violent. The standard behavior justifies the usage of the CN method and is often treated as a calibration [6] method.
The Zagożdżonka watershed, presented in this paper, is, since 1962, an experimental watershed of the Department of Water Engineering of the Warsaw University of Life Sciences (Szkoła Główna Gospodarstwa Wiejskiego), Poland. The studies of CN in this watershed were conducted mainly by Banasik, Ignar and Woodward, e.g., [20,21], and resulted in a determination of the tabulated value. It was here that the parameter CN was also investigated directly from rainfall-runoff events. The studies revealed the seasonal variations and the standard behavior of CN in the asymptotic approach.
The aim of the paper was to assess the variability of CN and contribute to a better understanding and estimation of CN for its practical application in engineering hydrology. In particular, a discordance of CN related to moderate and heavy rainfalls was analyzed. The variability of CN is a known fact; therefore, the authors' intentions were rather: (i) to assess the accuracy of CN tab and the antecedent runoff conditions (ARCs) of Hawkins, under very heavy or moderate rainfall; and (ii) to assess whether the application of the ARCs of Hawkins is justified, even for a well-calibrated model. To examine and describe the variability of CN , the probabilistic approach had to be taken and statistical methods used. This enabled the estimation of ranges of values within which the true CN was expected to lie with high probability. Generally, these ranges are crucial for design engineers, because ranges of runoff volumes, especially flood volumes, can be further calculated through the application of the SCS-CN model. The research presented in this paper related to investigations carried out by Tedela et al. [22]. They revealed, using statistical methods, a serious inaccuracy in the estimation runoff in most of the examined catchments when tabulated numbers were applied and claimed that many watersheds require variable curve numbers, associated with different magnitudes of rainfalls.
The paper is structured as follows: in the second part, the watershed is described; in the third section, the data and methods used are depicted; in the fourth section, the results are presented and discussed, and the last part contains the conclusions of the study.

Watershed Characteristics
The Zagożdżonka watershed is located in central Poland, with an area of 23.4 km 2 to the Czarna gauging station (Figure 1). Its topography is typically lowland, with the absolute relief ranging from 147.5 m to 184.5 m. The mean slopes of the main streams range from 2.5 m to 3.5 m per 1,000 m. A significant part of the area to the Czarna gauging station, approximately 16.8%, constitute local depressions, which do not contribute to direct runoff; therefore, only the hydrologically active part of the watershed was considered in the following study. Seventy percent of this part covers arable land with potatoes, grain and alfalfa, 9.4% pastures and 20% forests. Loamy sandy soils are abundant here, with soil groups B and D, with respect to the hydrologic soil groups introduced in [7]. In this system, recommended by the SCS-CN method [6], soils are placed into four groups, A, B, C, D, with A being the least runoff-prone and D the most runoff-prone. The tabulated CN value for the Zagożdżonka watershed, calculated as the weighted average from soil groups and land use, is CN tab = 74.6, and the retention parameter is S tab = 86.5 mm [21,23]. The mean annual runoff coefficient is low, approximately 0.17 at the Płachty gauging station, 2.3 km downstream from Czarna. This is because of the sandy soils, the lowland character of the watershed, as well as agricultural land use, with some parts being forests and wetlands. In the period 1962-2011, the total annual precipitation amount varied from 414 mm to 941 mm, with a mean value of 612 mm. The winter season (November-April) was characterized by the average precipitation amount of 223 mm, which was lower than 389 mm for the summer season. The driest month is January (32 mm, which is 5.2% of the annual value), and the wettest is July (82 mm, 12.9%). The long-term runoff varied from 52 mm to 209 mm, with a mean value of 107 mm. The month with the largest/lowest runoff is March, with 15.5 mm and July, with 5.8 mm, respectively. Analogously the month with the largest/lowest discharge at Płachty is March (0.477 m 3 /s) and July (0.166 m 3 /s), respectively, and the mean discharge is 0.277 m 3 /s. The SCS-CN method could be adapted here, due to the small area and agriculture-dominated land use [16]. The hydrological regime of the watershed was exhaustively presented in [21,24,25]. It is worth pointing out that in the studies of Banasik and Woodward [21], the standard behavior of CN in the asymptotic approach was revealed.

Material and Methods
The study was carried out using empirical CN and S, obtained directly from observed rainfall-runoff (P, H) pairs. The heavy rain totals were recorded at the precipitation gauging stations located in the watershed. The amounts were checked, and only clearly defined, single rainfall events from the period 1980-2008 were chosen. These events had to be visibly distinguishable in the records of runoff data and independent of other events. Subsequently, for each event, the base flow and direct runoff were separated by drawing a straight line from the point of rise to the beginning point of a fitted part of a groundwater recession curve matched to the recession segment; the second point was specified on the semilogarithmic plot of discharge [26]. Initially, a data set of 62 events, described in the work of Banasik and Woodward [21], was used for further analysis. Considering the fact that in agricultural lowland catchments, even large rainfall events do not cause surface/direct runoff, especially in the vegetation period, when the rainfalls usually appear [21,27], the set of 62 rainfall-runoff events recorded during the period 1980-2008 is regarded as a large one.
The calibration of the initial abstraction coefficient, λ, used the asymptotic value of CN , due to the asymptotic fitting method, and the medians of the CN samples, calculated from Equations (1) and (2). Both the asymptotic value and median were calculated for λ varying from 0.05 to 0.25. Among them, λ = 0.2 resulted in a CN nearest to CN tab ; therefore, λ = 0.2 was fixed in this study. In this case, 0.2S theor = 17; hence, events with P > 17 mm were selected for the main (P, H) series, while events with lower rainfall were omitted. This curtailed significantly the number of events to 34, but preserved the initial requirement in the SCS-CN method to estimate the direct runoff for stormy rainfalls [28]. On the other hand, the consequence of a reduction of the sample was a lessening of the power of the statistical methods, which were intended to apply. This problem is more widely discussed in Section 3.1.
Eventually, the main CN and S samples were calculated from Equations (1) and (2) for λ = 0.2: for P > 17 0 otherwise This data set of 34 elements selected for further analysis had a large range of rainfall events, i.e., from 18 mm to 96 mm. As we intended to assess the values of CN and S associated with different magnitudes of rainfalls, we divided this set of 34 elements with P 50% = 37 mm. The rainfall events of a depth not larger than 37 mm were considered as moderate rainfalls and larger than 37 mm as heavy rainfalls. We used another set of data from the period 1963-2011, namely the annual maxima of daily rainfall depth for the estimation of P 50% . Eventually, three cases, A, B, C, of CN s and Ss were considered in the study, of P > 17, 17 < P ≤ 37 and P > 37. The sets of CN s and Ss contained 34, 22 and 12 elements, respectively. Cases B and C expressed the response of the catchment to moderate and to heavy rainfalls, respectively.

Identification of the Theoretical Distribution
The identification of the theoretical distribution function was executed with statistical tests, to estimate the quantiles of CN directly from the pdf and to assess other characteristics of CN . Both distributions of CN and S can be identified in two ways: (i) for each variable separately or; (ii) when one distribution is known, for example the distribution of S, the distribution of CN can be obtained from the relationship between CN and S. To explain the second method more precisely, we denote, by g S , g CN , the densities of S and CN and; additionally, define y = p(x) = 25400 x+254 , as in Equation (5). Then, q, the inverse function of p, is expressed as q(y) = 254( 100 y − 1) (Equation (2)). The relationship between the densities is given by g CN (y) = g S (q(y)) · |q (y)| = y 2 25400 · g S ( 25400 y − 254) [29]. If the distribution of S is from a known family (e.g., lognormal, GEV, Weibull, etc.); then, this formula is impractical, being computationally highly complex and not leading to a distribution from a known family. This is because the distributions of S and CN were identified separately in this paper.
Empirical CN s were characterized by a negative skewness; therefore, the variables, 100-CN (instead of CN ), were considered to provide a positive skewness. This also allowed a choice of the parent distribution function of 100-CN between a large number of existing theoretical, positively skewed distribution functions. This is the simplest transformation, which leads to a positively skewed variable, with a lower limit equal to zero (because of the upper limit of 100 for CN ). This idea was introduced by McCuen [14]. Empirical Ss were positively skewed, and the distribution was directly identified.
Subsequently, the families F of competing theoretical distributions of S and 100-CN were built. The histograms and Cunnane plotting position plots had a supporting role in this step. The parameters of the distributions were estimated using the maximum likelihood estimation (MLE) method [30,31]. To explain the method, assume that {x 1 , ..., x n } are the observed values and that θ is the vector of unknown parameters of a theoretical density, g. The likelihood function is defined as: The maximum of L is achieved at someθ. This value is the estimate of θ in the MLE method. Often, ln L is considered, instead of L, to provide an effective method for estimation.
The candidates were the distributions, which passed the Kolmogorov-Smirnov (KS) [32], Cramer-von Mises (CM) [33] and Anderson-Darling (AD) [33,34] tests at a significance level equal to 0.05. Subsequently, the distributions, which achieved p−values of not less than 0.8 in all three tests were selected for the families F . This requirement was considered as the entrance condition.
The above three tests compare G(x, θ) and G n (x), the hypothetical and empirical cumulative distribution functions (cdf). The formulas for the test statistics are the following [35]: The AD test gives more weight to the tails than the KS and CM tests do, which was essential in further interval estimation. Therefore, the AD test was joined to the study, although it was not used in the further selection criteria. If the distribution was supposed to be normal, then the sample was additionally required to pass the Shapiro-Wilk (SW) test [36].
The identification of the distribution closest to the one supposed to have generated the data was performed using two criteria. The criterion W is based on the KS and CM statistics and on r, the correlation coefficient between the theoretical and empirical quantiles. It takes the form: This value is recommended for flood discharge estimation [37]. The second measure is the Akaike information criterion [38], which is also applied to flood frequency analysis [39,40]: where g(x i ) is the theoretical probability density function evaluated at x i and p is the number of parameters of the theoretical distribution. The lowest W and AIC indicate the best fit.
We faced the problem that the above tests allowed for accepting various distributions. Furthermore, measures W and AIC, even though they selected the best fit, did not show apparent differences between distributions. This is caused by the fact that statistical tests and W and AIC are not able to discriminate between the exact shapes of the distributions on the basis of small samples (especially in Cases B and C), because: (i) the estimates of parameters for each distribution are biased; (ii) several distributions may agree in the range of small samples; and (iii) the power of the goodness of fit tests, understood as the effectiveness in the discrimination between distributions, is low. Therefore, the tests, rather, indicated an acceptable group of distributions, because the selection decision of the one best distribution was risky in this situation. However, this decision was not needed, because, in fact, we did not require the exact information about the distribution type, but rather, sought confidence intervals, even if they were approximately similar for different distributions.

Confidence Intervals and ARC
The confidence interval (a, b) (CI) covers the unknown value with a large probability, say 95% or 99%, and is constructed directly from the theoretical cdf. The Hjelmfelt proposition [13] for the determination of CN (I) and CN (III) is based on the statistical approach [7]. Hjelmfelt recognized a lognormal distribution of S, used the logarithmic transformation of S, estimated the 10% and 90% quantiles of the fitted normal distribution and then used the reverse logarithmic transformation. After this transformation, he had two estimates of the quantiles of S. These quantiles, equal to S(III) and S(I), respectively, represented wet or dry conditions in the catchment. Subsequently, CN (I) and CN (III) were computed from Equation (5). They always provide CN values that are located in the 95% and 99% CIs (of CN ). Hence, they estimate error bands of CN [6]. Observe that the ARCs of Hjelmfelt vary with the sample, while the ARCs of Hawkins do not.
In the following study, the assumption of lognormality of S was verified; the 10% and 90% quantiles were estimated directly from the theoretical cdf of S and transformed to CN (III) and CN (I) from Equation (5). Then, they were compared to the 90% and 10% quantiles of CN calculated directly from the distribution of 100-CN . They showed a good level of agreement. Hence, the Hjelmfelt method led to: where x 10% and x 90% are the 10% and 90% theoretical quantiles of CN , which are the 90% and 10% quantiles of 100-CN . Subsequently, a comparison between CN (I) and CN (III) of Hawkins and of Hjelmfelt was completed.

Results and Discussion
The results of the analysis of Samples A, B and C were collected in Table 1. A good agreement (fewer than 1.5 CN ) can be observed in Case A between CN tab = 74.6 and the mean and median. Recall that, as was mentioned in Section 1, the asymptotic fitting method also resulted in a value similar to CN tab , namely 73.4. These three methods, the median, the arithmetic mean and the asymptotic fitting, are treated as methods of calibration of CN [6,7,41]. Therefore, CN tab was assumed to be unquestionable in this paper as a very good calibrated parameter. These results differed from those obtained in [22], where empirical CN s were studied for 10 forested-mountainous watersheds in the USA. For several watersheds, a serious disagreement between the tabulated and the asymptotic or median/average CN s was observed, reaching over a dozen or several dozen CN s. Soil misclassification and differences in watersheds elevation were supposed to cause these discrepancies in [22].

Assessment of the Accuracy of CN tab and the ARCs of Hawkins under Very Heavy and Moderate Rainfall
Division into two subsamples revealed a discrepancy: the mean and median were larger in B and smaller in C when compared to CN tab , with moderate difference, i.e., greater than five and less than 10. The medians and means in B and C differed, however, by more than 10 CN . For the parameter, S, the relationship was reversed. The samples were characterized by a moderate variability and by a large range. In particular, the range of CN and S in Subsamples B and C was wide. The skewness in all three cases was negative for CN , being not substantially different from zero in Case B and high in C, which was caused by a low CN during heavy rainfall. The parameter, S, was positively skewed.
The families, F S , F 100−CN , of competing distributions of the variables, S and 100-CN , were defined as: where the following abbreviations are used: GEV, generalized extreme value; N, normal; GLO, generalized logistic; WE, Weibull; LN, lognormal; G, gamma. Both measures W and AIC are presented in Table 2. Obviously, the distributions differ by parameters, e.g., GEV in Case A has other parameters than in B and C. Criterion W indicated GLO as the first parent distribution of S, GEV as the second and LN as the third in A, B and C. Measure AIC, in turn, recognized the lognormal distribution as the best fit, and GEV and GLO as nearly equivalent. The distribution G gave a comparable fit to LN in Cases A and C. The fit with WE was just somewhat worse than that with LN and G. The fit with WE and G in Case B was not considered, as the condition with p-value was not fulfilled. The results partially agree with Hjelmfelt's [13], who identified LN as the parent distribution function of S.
Criterion W aimed at GEV for the first parent distribution of 100-CN in Cases A, B, C. It located N as the second and GLO as the third in A, GLO as the second and N as the third in C and N and GLO as equivalent in B. In all three cases, the AIC measure recognized N as the best fit, then GEV and GLO. The distributions LN, WE and G gave a slightly lower fit and were considered only in C.
The differences in measures were small between distributions; therefore, the first three, GEV, GLO, LN for S and GEV, N, GLO for 100-CN , were allowed for further study. For safety, to overcome the problem of non-uniqueness of the method discussed in Section 3.1, the CIs were constructed for each of three distributions, GEV, N and GLO for CN , by transforming them from 100-CN .
The GEV distribution is featured by a thick, moderate or thin tail. It relates to the shape parameter greater than zero (Frechet distribution), when the pdf decreases slower than exponential, equal to zero (Gumbel distribution) or lower than zero when the pdf has an upper limit (reversed Weibull distribution). In Cases A, B and C, the results for 100-CN indicated a reversed Weibull. For parameter S, the reversed Weibull distribution was obtained for B. Frechet was obtained for both A and C, which was caused by the large values in the samples. These very large values of S related to heavy rainfall.
A similar analysis, based on events with the largest discharges, was performed by McCuen [14], who considered empirical CN in five watersheds in the Piedmont Region, USA. In his study, the gamma distribution gave the best fit to the 100-CN data.
Subsequently, the 95% and 99% confidence limits of CN were computed, with GEV, N and GLO for 100-CN as the probable parent distribution. Next, to apply the Hjelmfelt method, the 10% and 90% quantiles of lognormally distributed S were calculated, and then, using Equation (2), the ARC for CN were determined. The 10% and 90% quantiles of CN were computed also directly from the parent distributions of 100-CN . A very good agreement was observed between the quantiles obtained using these two methods, with differences less than 0.5. This procedure was repeated for GLO as the first theoretical distribution function of S, due to measure W. The results were analogous, which implied that the Hjelmfelt method of quantiles of S could be replaced by the method of quantiles of CN , as they provided very similar CN (I) and CN (III) values. Eventually, the 10% and 90% GEV quantiles, denoted by CN (I) Hjel and CN (III) Hjel , were applied in subsequent analysis.
Equations (12) and (13) The above values were located approximately symmetrically around CN tab . Table 3 presents the characteristics resulting from distribution fitting. They are depicted in Figure 2 to enable the comparison between them and the ARCs of Hawkins. Table 3 and Figure 2 showed a large disagreement between medians of CN for moderate and heavy rainfalls, reaching 15. This related also to differences in the ARCs of Hjelmfelt. The theoretical range of CN was larger in Case B than in C. Below, A, B and C are discussed.
Case A (P > 17): The median of the distribution approximately matched CN tab , which confirmed the right determination of the tabulated value. The 95% and 99% confidence limits were similar for both GEV and N. The lower confidence limits obtained from GLO were lower than those from GEV and N, which was caused by a long tail of the GLO pdf. Both CIs covered CN tab , CN (I) Haw and CN (III) Haw for any selected distribution, and additionally, the differences between CN (I) Haw and CN (I) Hjel and between CN (III) Haw and CN (III) Hjel were only minor. The results obtained in this case suggest that the Hawkins formulas perfectly reflect the error bands.
Case B (17 < P ≤ 37): CN tab was located in the lower part of both CIs, being less than the median of each theoretical distribution. The disagreement between CN tab and the median was 7.3 for GEV and N and 7.2 for GLO. For all distributions, the CN (I) Haw was out of the 95% CI and was lower than the 99% lower confidence limit for both GEV and N. Hence, the CN (I) Haw was inappropriate, being too low. The CN (III) Haw was located in the upper part of each CI and was lower by four than CN (III) Hjel . Generally, the CN tab and CN (I) Haw were strongly underestimated, and CN (III) Haw was slightly underestimated for moderate rainfalls.
Case C (P > 37): CN tab was located in the upper part of both CIs, being larger than the theoretical median. This was the source of a discrepancy between CN tab and the median, which reached 8.1 for GEV, 8.5 for N and 7.3 for GLO. As CN tab nearly reached CN (III) Hjel , then, for about 90% of heavy rainfalls, CN was suspected to be lower than CN tab . CN (III) Haw was out of the 99% CI, and the difference between CN (III) Haw and CN (III) Hjel was large: 11.9 for GEV, 11.2 for N and 12.1 for GLO. CN (I) Haw , located in the lower part of the CIs, was greater by three than CN (I) Hjel . Hence, CN tab and CN (III) Haw strongly overestimated, and CN (I) Haw slightly overestimated CN for heavy rainfalls.

Assessment of the Applicability of the ARCs of Hawkins
A general observation is that in Case A, both of the ARCs of Hawkins were covered by all three CIs, while in Case B, the ARC for dry conditions and in C the ARC for wet conditions were out of the CIs. Large differences between the ARCs of Hawkins and Hjelmfelt indicate the inapplicability of CN (I) Haw for moderate and CN (III) Haw for heavy rainfall. The usefulness of CN tab for heavy rainfall is also doubtful.
The study completed by McCuen [14], mentioned above, also led to the construction of confidence intervals. The ARCs of Hawkins provided ranges that entirely covered 95% CIs, which was unlike our results, where CIs in Case A had a broader range than that given by the ARCs, and only one end point of the CI was covered by the range given by the ARCs in cases B and C.
To support our conclusions, an analogous study was conducted for the Stupavský catchment, located in the Inner Western Carpathians in Slovakia. The stream length is 16.7 km and the catchment area 33 km 2 for the gauging station. The catchment is forested, with a forest cover ratio equal to 0.9. The tabulated CN and the ARCs of Hawkins are CN tab = 67, CN (I) Haw = 49.4, CN (III) Haw = 82.6. Only one-day precipitation events were selected for the study. The two-year daily maximum rainfall depth was estimated as 46.6 from the observation period 1983-2006. Twenty eight rainfall-runoff data samples were selected for the main sample. After division, the samples with P ≤ 46.6 (B, 16 elements) and with P > 46.6 (C, 12 elements) were created. Afterwards, the GEV distribution was fitted and the confidence intervals constructed. The value CN (I) Haw = 49.4 was not covered by the 99% confidence interval (58.8; 77.2) for P ≤ 46.6. For P > 46.6, both CN tab = 67 and CN (III) Haw = 82.6 were not covered by the 99% confidence interval (30.0; 66.2). Therefore, the use of CN (I) Haw for P ≤ 46.6 and CN (III) Haw for P > 46.6 is suspected to be inappropriate.
A comparison of the ARCs of Hawkins and Hjelmfelt is the motivation for the formulation of the statement that even the perfect calibration of the SCS-CN model, based on soil and land use and on empirical CN for various rainfall conditions, may not provide the best possible estimation of the curve number. Different responses of the catchment to moderate and heavy rainfall do not justify a thoughtless use of the ARCs of Hawkins, but require a good understanding of the properties of CN .

Conclusions
The conclusions, presented below, are catchment specific. To observe the CN behavior in other catchments, similar investigations should be carried out using empirical CN values. Such further research would allow for estimating design floods in other small ungauged catchments, in order to inform flood protection measures.
The main conclusion follows from the comparisons of CN tab and the median and of CN (I), CN (III) of Hawkins and of Hjelmfelt. When the sample was considered as a whole (Case A), CN was correctly calibrated, and a very good agreement between CN tab and the median and between CN (I), CN (III) of Hawkins and of Hjelmfelt held. This agreement was reliable, as it was observable for all three distribution functions. However, when the sample was divided in B∪C, a serious inconsistency arose, especially between the CN (III) of Hawkins and Hjelmfelt in Case C. Additionally, CN (III) Haw was located outside the 99% CIs of three distributions; hence, it was close to impossible.
Similar considerations applied to Case B with the CN (I) of Hawkins and Hjelmfelt, although this disagreement was lower.
Generally, ARCs should be applied with caution. Even if they agree with the 10% and 90% quantiles for all collected rainfall-runoff data, which is what suggests the justification for their use, a serious inconsistency is observed in subsamples. It is worth pointing out that, although the inconsistency of CN for moderate and heavy rainfalls is generally a known fact, the above probabilistic approach allowed us to make constructive corollaries of the practical usage of ARCs: • instead of CN (I) Haw = 59.6 and for a rainfall depth not much higher than 17 mm, consider a value between 70 or 71, which is near CN (I) Hjel for moderate rainfall, • instead of CN (III) Haw = 87.3 and for a very high rainfall depth, consider a value between 75 or 76, which is near CN (III) Hjel for stormy rainfall.
Moreover, the Hjelmfelt method may be simplified. The direct calculation of quantiles of CN is advisable instead of using the transformed quantiles of S. Additionally, as the normal distribution gave a reasonable fit, the direct use of its cdf substantially simplifies the estimation of the quantiles. The normality of CN should, however, still be verified when new data is obtained. The method of Hjelmfelt may be applied even when (P, H) data are not accessible. Where data is not available, an assessment of the hydrologic similarity of the watershed to that of another watershed for which the distribution of S and CN is known can be performed to determine its CN values.