Improvement of SCS-CN Initial Abstraction Coe ﬃ cient in the Czech Republic: A Study of Five Catchments

: The SCS-CN method is a globally known procedure used primarily for direct-runo ﬀ estimates. It also is integrated in many modelling applications. However, the method was developed in speciﬁc geographical conditions, often making its universal applicability problematic. This study aims to determine appropriate values of initial abstraction coe ﬃ cients λ and curve numbers ( CNs) , based on measured data in ﬁve experimental catchments in the Czech Republic, well representing the physiographic conditions in Central Europe, to improve direct-runo ﬀ estimates. Captured rainfall-runo ﬀ events were split into calibration and validation datasets. The calibration dataset was analysed by applying three approaches: (1) Modifying λ , both discrete and interpolated, using the tabulated CN values; (2) event analysis based on accumulated rainfall depth at the moment runo ﬀ starts to form; and (3) model ﬁtting, an iterative procedure, to search for a pair of λ , S ( CN , respectively). To assess individual rainfall characteristics’ possible inﬂuence, a principal component analysis and cluster analysis were conducted. The results indicate that the CN method in its traditional arrangement is not very applicable in the ﬁve experimental catchments and demands corresponding modiﬁcations to determine λ and CN (or S , respectively). Both λ and CN should be viewed as ﬂexible, catchment-dependent (regional) parameters, rather than ﬁxed values. The acquired ﬁndings show the need for a systematic yet site-speciﬁc revision of the traditional CN method, which may help to improve the accuracy of CN -based rainfall-runo ﬀ modelling.


Introduction
The Soil Conservation Service Curve Number (SCS-CN) method is an empirical lumped rainfall-runoff model developed in 1954 by the US Department of Agriculture's Soil Conservation Service (USDA SCS, now the Natural Resources Conservation Service NRCS). Originally, the method entailed converting basic descriptive inputs characterising catchment features into numeric values of curve numbers (CNs) [1]. For any given rainfall amount, the corresponding runoff level is estimated. The procedure is based on a combination of the water-balance equation and two fundamental assumptions, in which Q is direct runoff (mm); P is precipitation (mm); Ia is initial abstraction (mm) that has to be exceeded so that direct runoff can start to form; F is cumulative infiltration (mm), not including Ia; S is maximum potential retention of a watershed after the start of the runoff and λ is initial abstraction ratio (dimensionless). The runoff depth Q is expressed using Equations (1)-(3) as for P > Ia, Q = 0 for P ≤ Ia.
The maximum potential retention S is expressed as a function in which CN is the curve number, theoretically ranging between 0 and 100. Over the years, the method has undergone numerous modifications [2][3][4][5][6], during which time the scientific community adopted the method globally for its simplicity, applicability and relatively low demands on input data [7][8][9][10][11].
Attention has been payed to the ability of a single CN value to characterise the runoff response of a watershed correctly. While the CN may be calculated from observed P-Q data, practical experience shows that the CN values vary between events [11,24]. Therefore, various approaches were introduced aiming to estimate an appropriate single CN value or single asymptotic CN, separately [1,[25][26][27][28][29]. However, the aforementioned methods showed certain shortcomings addressed by a study that analysed the effect of soil and land cover complexes on the hydrologic responses in watersheds [30]. Specifically, the authors proposed a concept of a two-CN heterogenous system. This study was subsequently generalised [31]. The latter two methods along with other CN estimation methods were used and compared in a study [32] of a watershed affected by a major forest fire, so that the effect of both spatial heterogeneity in land cover and the fire could be investigated.
Another problematic aspect of the CN method lies in the characterisation of previous rainfall. It is expressed using three antecedent moisture condition (AMC) classes (I: Dry; II: Normal; and III: Wet), each of which has its own set of CNs. The discrete character of the CN-AMC relationship is manifested by sudden jumps in estimated runoff [33]. On the basis of previous five-day rainfall totals, CN II is converted to CN I and CN III using relationships proposed by [24][25][26][27][28][29][30][31][32][33][34][35][36][37] and directly from the NEH-4 tables [4,38,39]. Despite the existence of different conversion methods, the discontinuity problem remained unresolved, and the concept of soil moisture accounting (SMA) became a subject of research. Williams and LaSeur [40] presented the first SMA concept. Their procedure was based on the idea that the CN varies continuously, along with soil moisture. Hawkins [41] then developed a continuous soil moisture accounting concept, but this modification did not take into account the SCS-CN procedure's structural shortcomings in terms of SMA, a problem addressed by Mishra et al. [42], who introduced an improved SMA-based SCS-CN model. Subsequently, Mishra and Singh [43] developed a versatile SCS-CN model (VSCS-CN) based on the assumption that higher antecedent soil moisture before a rainfall event will provide greater runoff and lower the soil's capacity to store and retain rainfall water.

Study Area
The study was carried out in five experimental watersheds in the Czech Republic (Figure 1), differing mainly in size and in land use. This made it possible to assess the CN method's performance and interpret the findings of this research with regard to each watershed's specific geographical conditions. All the studied watersheds' basic characteristics are provided in Table 1. widen downstream, and the surface becomes slightly curved to flat. Two main soil types are present in the catchment. Cambisols are found mainly in the upper parts, while luvisols are prevalent in the lower parts. As for land use, a specific spatial distribution exists. Arable land, usually in the form of large plots, is found both in the upper and lower parts of the catchment on the flat relief. Forests and meadows are located mostly in the transition between the upper and lower parts, in the valleys and along streams. This specific environment creates favourable conditions for generating flash floods from intense rainfall. Two major flash-flood events were recorded in summer 2009 and 2010 [72,73]. On 2 July 2009, the region experienced extreme rainfall (local precipitation reached up to 50 mm) after a very rainy period, causing a severe flood event with peak discharge of 33.1 m 3 /s and corresponding to a 50-year flood event. On 1 June 2010, another rainfall event was recorded (with a maximum raingauge measurement of 55 mm), but due to the precipitation's differing spatial distributions and lower previous rainfall, peak discharge reached only 14.4 m 3 /s, corresponding to the return period of 5-10 years. Nevertheless, both events resulted in heavy damage to municipalities along the stream.

Suchý Creek Catchment (NE)
The study area is located approximately 30 km north of Brno in the Drahany Highlands. The catchment is drained by a minor stream, Suchý Creek, and forests cover the western part of the watershed. More than half the area is arable land. Due to large fields, soil erosion poses an imminent threat to a little dam that was constructed right under the catchment's outlet to protect the areas downstream from flash floods [74].

Kopaninský Creek Catchment (KO)
This catchment is located in the Křemešnice Highlands in the southern part of the Sázava River basin, a right-side tributary of the Moldau (Vltava) River. It is a source of drinking water for the City of Prague, located approximately 80 km to the northwest. From a geological perspective, this is a very old area formed by igneous and metamorphic rocks during the Paleozoic Era, including granite and schist. Rolling hills and low mountains dominate this region's surface, with the altitude mostly ranging  Table 1. Basic characteristics of the experimental catchments and percentage of land use categories (AL-arable land; BG-bare ground; FO-forest; GR-grassland and meadows; GA-gardens and orchards; PA-paved and built-up areas; SH-shrubbery; WA-water bodies).

Husí Creek Catchment (FU)
Husí Creek is a left-side tributary of the Oder (Odra) River, located in northeastern Czech Republic, approximately 30 km southwest of the city of Ostrava. The upper parts belong to the easternmost edge of the Bohemian Massif in the Nízký Jeseník Highlands, while the lower parts lie in the Outer Subcarpathia, a depression at the outer base of the Carpathian arc. The catchment's position predetermines its specific relief; it is characterised by a flat or slightly wavy surface in the upper parts that gradually passes into narrowing valleys with greater longitudinal slope. The valleys widen downstream, and the surface becomes slightly curved to flat. Two main soil types are present in the catchment. Cambisols are found mainly in the upper parts, while luvisols are prevalent in the lower parts. As for land use, a specific spatial distribution exists. Arable land, usually in the form of large plots, is found both in the upper and lower parts of the catchment on the flat relief. Forests and meadows are located mostly in the transition between the upper and lower parts, in the valleys and along streams. This specific environment creates favourable conditions for generating flash floods from intense rainfall. Two major flash-flood events were recorded in summer 2009 and 2010 [71,72]. On 2 July 2009, the region experienced extreme rainfall (local precipitation reached up to 50 mm) after a very rainy period, causing a severe flood event with peak discharge of 33.1 m 3 /s and corresponding to a 50-year flood event. On 1 June 2010, another rainfall event was recorded (with a maximum rain-gauge measurement of 55 mm), but due to the precipitation's differing spatial distributions and lower previous rainfall, peak discharge reached only 14.4 m 3 /s, corresponding to the return period of 5-10 years. Nevertheless, both events resulted in heavy damage to municipalities along the stream.

Suchý Creek Catchment (NE)
The study area is located approximately 30 km north of Brno in the Drahany Highlands. The catchment is drained by a minor stream, Suchý Creek, and forests cover the western part of the watershed. More than half the area is arable land. Due to large fields, soil erosion poses an imminent threat to a little dam that was constructed right under the catchment's outlet to protect the areas downstream from flash floods [73].

Kopaninský Creek Catchment (KO)
This catchment is located in the Křemešnice Highlands in the southern part of the Sázava River basin, a right-side tributary of the Moldau (Vltava) River. It is a source of drinking water for the City of Prague, located approximately 80 km to the northwest. From a geological perspective, this is a very old area formed by igneous and metamorphic rocks during the Paleozoic Era, including granite and schist. Rolling hills and low mountains dominate this region's surface, with the altitude mostly ranging between 500 and 650 m. As for the soils, cambisols are found most frequently here, along with pseudogleys and gleys, in valleys and along streams. The landscape is covered predominantly by arable land, forests and, to a lesser extent, grasslands and meadows. The settlement is scattered in rather small villages.
The KO catchment is drained by the Kopaninský Creek. The main catchment outlet is KO-3, another two sub-catchments analysed in this study are KO-2 and KO-1. The KO-2 catchment is located in the uppermost part of the KO-3 catchment and has the greatest portion of forested areas, at more than 44%. The KO-1 catchment, located in the lower part, is covered almost entirely by arable land (more than 97%) and is partly tile-drained.
All five studied watersheds' basic geographical characteristics are listed in Table 1.

Data Collection
Automatic rain gauges were used to measure precipitation heights. The monitoring network included up to 17 rain gauges (automatic 0.2 mm tipping buckets 200 cm 2 and 0.1 tipping buckets 500 cm 2 ) in catchment FU and its nearby surroundings. The Thiessen polygon interpolation procedure was used to calculate average precipitation levels [74]. A single rain gauge (automatic 0.2 mm tipping bucket 200 cm 2 ) was used to measure rainfall in catchments NE and KO. Using the precipitation data, individual rainfall events were identified, for which five-and 10-day antecedent precipitation totals (P 5d , P 10d ) were calculated.
The water level data were measured using calibrated weirs with ultrasound probes (Fiedler company)-a Thompson weir was used with in catchments NE, KO-1 and KO-2, a compound Cipolletti weir was used in catchment KO-3. In catchment FU, water levels were measured using hydrostatic level sensors (Fiedler company). Discharge values were calculated by means of rating curves set for all monitored outlets on the basis of hydrometric measurements, levelling and using the Chézy formula [75].
An overview of measuring devices is provided in Table 2.

Estimate of CN Values
The hydrologic soil groups (HSGs) were derived based on Valuated Soil-Ecological Units (VSEUs), a whole-country database (1:5000) determined by the Research Institute for Soil and Water Conservation (RISWC). Land-cover units were identified based on aerial images and a reconnaissance field survey. The corresponding percentages are provided in Table 1, and the spatial distribution is depicted in Figure 2. Moreover, cropping and tillage history also was considered in the NE catchment. Subsequently, the final identification of CNs for each unit was done according to Janeček et al.'s [56] methodology, which frequently is applied in the Czech Republic, in which the tabulated CN values correspond to different HSGs, land-use units and crop types. Moreover, moisture conditions are reflected through antecedent moisture conditions (AMCs). AMCs are characterised by three degrees (I: Dry; II: Normal; and III: Wet), delimited by a five-day precipitation total of P 5d of 36 mm and 53 mm prior to a rainfall (P 5d < 36 mm-dry conditions; P 5d ≥ 53 mm-wet conditions). Maps of CN II values are provided in Figure 2 for each catchment.  The aforementioned parameters were applied in the analyses described in detail in the following sections. To verify the adjusted CN method's performance, events' datasets were split into training and validating datasets (see the numbers provided in Table 2), so that they were similar in terms of level and variability of P, P5d and P10d.

Principal Component Analysis, Cluster Analysis
In order to reveal possible concealed relationships between parameters of rainfall-runoff events and to assess objectively the role that the individual parameters may play in the rainfall-runoff process,

Events, Determination of λ
To avoid biased λ values, rainfall events recorded between April and September with P ≥ 5 mm were chosen for the analysis. In order to exclude rainfall events that caused compound hydrographs, rainfall-runoff events were selected individually, so that the influence of runoff from previous rainfall events could be obviated. For each event, the following parameters were determined: precipitation level P (mm), 2.
hydrograph was drawn and base flow was separated from the hydrograph based on the hydrograph's inflection points, which indicated the start of the rising limb and the end of the recession limb for each runoff event [69,70], 6.
initial abstraction Ia (mm), which was calculated as the rainfall level at the moment when direct runoff begins. Such a procedure has been used for instance in research by Shi et al. [11]. However, in large watersheds and/or in case of uneven spatial distribution of precipitation such a calculation of Ia might be problematic as the runoff needs some time to reach the watershed's outlet [24].
The aforementioned parameters were applied in the analyses described in detail in the following sections. To verify the adjusted CN method's performance, events' datasets were split into training and validating datasets (see the numbers provided in Table 2), so that they were similar in terms of level and variability of P, P 5d and P 10d .

Principal Component Analysis, Cluster Analysis
In order to reveal possible concealed relationships between parameters of rainfall-runoff events and to assess objectively the role that the individual parameters may play in the rainfall-runoff process, principal component analysis (PCA) and cluster analysis (Ward's method) were conducted for the training datasets from all experimental catchments. The cluster analysis allowed to explore groups (clusters), which the events tend to form. The optimal number of clusters was set using Pham et al.'s [76] k-selection method.

Modified λ with Tabulated CNs
The procedure begins with the determination of tabulated CN values [65] for individual units in a catchment following the procedure described in Section 3.2. For each unit, from the CN values, maximal potential retention S values are derived using Equation (5). Average S for the entire watershed area is calculated using weighted arithmetic mean. In the following step, the corresponding initial abstraction coefficient is calculated using Equation (4), which is converted into the quadratic formula producing a pair of initial abstraction coefficient values for each event. From these, the final selection is carried out following the assumption that the lowest nonnegative value must be chosen. Numerically, the solution of the quadratic equation is correct. However, initial abstraction coefficient λ must be higher than or equal to zero. In case of two non-negative roots, choosing the higher one causes initial abstraction Ia being too high, resulting in zero direct runoff generated as the condition from Equation (4) is not met. The determination of S values and subsequent calculation (Equation [6]) were conducted twice: (1) Taking into consideration the AMC classes (i.e., using CN I , CN II and CN III ); and (2) neglecting AMC classes and using only CN II , assuming that the AMC determination might not be a deciding factor in calculating direct runoff.
Discrete λ This procedure's objective was to identify a discrete representative λ for a watershed calculated as mean and median of values obtained using Equation (6). Mean and median λ were calculated for the entire training dataset of events and for clusters delimited by the cluster analysis.
Interpolated λ This procedure's objective was to express λ as a linear function of P [66], instead of using a single value. The regression equations were derived using λ:P plots. Again, the calculations were carried out with and without determination of AMC.

λ Modifications Not Dependent on Tabulated CNs Event Analysis
This procedure was an inverted process. The maximum potential retention, S, and λ were calculated on the basis of rainfall and runoff parameters. The calculation included the following steps: 1.
The total runoff Q and the initial abstraction Ia were calculated using the procedures described above (Section 3.3).

2.
The maximum potential retention, S, was an unknown parameter that was calculated using equation obtained by modifying Equation (4).

3.
The initial abstraction coefficient λ values were determined by dividing Ia by S for each rainfall-runoff event. Mean and median λ were calculated for each experimental watershed. 4.
Regression of S, according to P 10d , was calculated and used for the validation together with mean and median λ.

Model Fitting
Model fitting is an iterative least-squares procedure for determining a pair of λ and S such that the sum of the squared differences between the runoff observed and modelled as a minimum. The computation of the model-fitting procedure was carried out using a loop-based script in Python 3.7.

Comparison of Estimated Runoff
For each experimental catchment, a training dataset of events was chosen to identify the key relationships characterising the CN method's proper setting. The remaining events were used for the parameter validation.
The quality of the runoff levels estimated by the modified CN method was evaluated by means of Nash-Sutcliffe model efficiency coefficient (NSE) (9), and bias (e) (10).
NSE is given as in which Q pi is the i-th modelled runoff depth, Q i is the i-th observed runoff depth and Q a is the mean observed runoff depth. Bias (e) is given as: NSE can range from −∞ to 1. An NSE equal to 1 expresses a perfect match of the modelled and observed values. An NSE ≥ 0.5 indicates a very good prediction ability for the model. A lower NSE indicates a worse performance by the simulations. The e is the average difference between the modelled and predicted runoff, thereby expressing the model's under-or over-estimating tendencies.

HEC-HMS Simulations
To test the applicability of the modified SCS-CN parameters for simulation of surface runoff, HEC-HMS 4.3 software was used. The simulations' aim was to verify the model's ability to deliver a runoff volume/level that did not differ from the observed one by more than 10%, with an NSE ≥ 0.500 (i.e., corresponding to the observed hydrograph's time course). Such simulations are viewed as sufficiently accurate. The computations were carried out only in the NE catchment, which is relatively small (A = 2.8 km 2 ), thereby allowing for better simulation of the surface runoff and identifying the surface-runoff simulations' response to changes in individual parameters provided in Table 3. Table 3. Overview of methods employed in surface runoff simulations using HEC-HMS. The catchment was split into three sub-basins. Calibration was carried out for two types of runoff events: (1) Those with a gradual onset, and (2) those with both a rapid onset and rapid recession limb.

Loss
The initial abstraction, Ia, was calculated using the curve numbers estimated for the sub-basins through the method described in Section 3.3, the λ values from the original SCS-CN method (λ = 0.2) and the median λ values (λ = 0.0142) from the analysis of which the results are provided in Section 4.2 The percentage of impervious surfaces (buildings and sealed roads) was calculated using ESRI ArcGIS 10.7 software. Time of concentration, T c , was calculated using Dingman's [77] formula, in which L is the longest flow distance path in a sub-basin (m), A is the sub-basin's area (km 2 ) and Y is mean slope (%). Dingman's formula (11) commonly is used in calculations by the Czech Hydrometeorological Institute (CHMI).
The initial discharge at each sub-basin's outlet was calculated using the main outlet's value according to the percentage of the sub-basin areas in the entire NE catchment. The remaining parameters (i.e., storage coefficient, recession constant and ratio) were set so that the calibration delivered the best possible result expressed by the bias of the runoff volume/level and NSE. Calibration coefficients were calculated for Ia, CN and T c as rates of the calibrated values to the original ones.
Finally, two sets of these values were available: (1) For the runoff events with slow onset, and (2) for runoff events with rapid onset.
In the subsequent step, two verification events were chosen-one with a slower rising limb and another with rapid onset. Having known the P 5d (AMC, respectively), the original CN, Ia and T c values were used for both λ = 0.2 and λ = 0.0142, corresponding to each verification event's year. These values were multiplied by the calibration coefficients acquired in the calibration simulations. The storage coefficient, the recession constant and the ratio were calculated in the calibrations.

Principal Component Analysis
To reveal various rainfall and runoff parameters' influence on runoff response to rainfall at event level, principal component analysis (PCA) and cluster analysis were conducted for the studied watersheds.
The PCA analyses' outputs are presented in Figure 3. Table 4 provides an overview of the variability percentages that the first two principal components (PC1 and PC2) explain, including all input parameters' loading values in relation to PC1 and PC2. The variability percentages explained by PC1 and PC2 differ between the studied catchments. In the FU catchment, PC1 showed a value of over 71% and PC2 less than 14%. In the other catchments, PC1 and PC2 percentage differences were not as big. The least difference was recorded in the KO-1 catchment (PC1 49%, PC2 40%). The vectors' orientation in Figure 3 and the loading values in Table 4 indicate that P 10d and P 5d present the most influential parameters in the direction of PC1, with maximal rain intensity and P being the most influential parameters in the direction of PC2. The KO-1 catchment has an opposite arrangement, as maximal rain intensity and P have the highest loading values in PC1 and P 10d and P 5d in PC2. This might be caused by the high percentage of arable land with presence of tile drainage system [78]. In such environment, P 10d and P 5d do not play the most influential role. Instead, direct runoff formation is dominated by the rain intensity along with P. Higher rain intensity produces faster and higher surface runoff and vice versa. This applies in general, but is probably accentuated in this specific catchment, making it different from the other ones.
To sum up, PCA indicates that P 10d and P 5d are the most influential parameters. The effects from maximal rain intensity and P, the other two most influential parameters, can be viewed as opposites, as they act along PC2.
Water 2020, 12, x FOR PEER REVIEW 11 of 28 To sum up, PCA indicates that P10d and P5d are the most influential parameters. The effects from maximal rain intensity and P, the other two most influential parameters, can be viewed as opposites, as they act along PC2.     2005), the optimal number of clusters was determined to be two for the experimental watersheds. Subsequently, non-hierarchical k-means clustering was conducted. Mean and median values of variables in the delimited clusters are listed in Table 5. A difference in P 5d and P 10d between the delimited clusters in the studied watersheds is apparent. Such a finding corresponds to the output of the cluster analysis. Given that P 10d was identified as the most important factor, P 10d threshold values were calculated for the catchments to be able to distinguish the cluster to which an event belongs. To determine thresholds, the mean of the upper quartile P 10d from the cluster with smaller P 10d and the lower quartile P 10d from the other cluster (i.e., a cluster with higher P 10d ) was calculated. The same procedure was applied to calculate the corresponding P 5d thresholds, as the P 5d normally is used in the traditional SCS-CN procedure to determine AMC classes. The P 5d and P 10d thresholds were 26 mm and 47 mm in the FU catchment, 18 mm and 39 mm in the NE catchment, 35 mm and 42 mm in the KO-1 catchment, 41 mm and 50 mm in the KO-2 catchment and 20 mm and 28 mm in the KO-3 catchment. For each catchment, these P 5d and P 10d thresholds were used for separation of rainfall-runoff events into two groups in analyses across the paper.

Tabulated CN-Based λ Modification
The first way to improve the CN method is through incorporation of tabulated CNs [65] along with modified λ values. The following calculations were carried out with and without taking into consideration the separation of rainfall events according to the P 5d and P 10d thresholds determined through the cluster analysis. Moreover, AMC classes' effect was examined by performing the calculations with and without the AMC classification. Modified λ was calculated as a discrete value and as functionally dependent on P.

Discrete λ
For the comparison, the original CN method with λ = 0.2 was applied, and the results showed that in this form, the model mostly was unable to deliver good results, except for the KO-1 catchment (for λ = 0.2 NSE = 0.861). Partial improvements were recorded in simulations with the AMC classification omitted. Only in the KO-1 catchment, the accuracy reduced significantly (NSE dropped from 0.861 to −16.769, see Table 6 and Figure S1a,b). It is important to emphasize that the traditional CN method is limited in terms of "zero runoff" cases since the condition from Equation (4) is often not met due to high λ, and also due to the effect of spatial variability in catchments [30,31], which, however, was not addressed in the present study. The relatively good results for the traditional CN method needs to be confronted with high numbers of events for which the model was unable to deliver non-zero direct runoff. This is apparent in Figure S1a,b. Specifically, the number of non-zero cases ranged from 0 to 13 for simulations with AMC classified and from 7 to 27 for simulations with AMC implicitly set as normal (AMC = II).
To improve the CN method's performance, mean and median λ values were determined in the training datasets as a whole and in groups separated according to the P 5d and P 10d thresholds. Generally, the analysis revealed that in most instances λ was lower than 0.2. In the unseparated datasets, mean λ ranged from 0.0365 to 0.1623 and median λ from 0.0142 to 0.0913. Except for one catchment (KO-1), the reduction of λ led to improvements in estimated runoff, since NSE exceeded 0.500, thereby indicating sufficient accuracy. In simulations with separation of events according to P 5d and P 10d thresholds, a greater variability in λ values existed, ranging from 0.0183 to 0.4231. Despite the expectation that this greater diversification of λ values will lead to more accurate results (in relation to the non-cluster-based approach-all, i.e., without the separation according to P 5d and P 10d ), this assumption only partially was fulfilled. In the FU catchment, the cluster-based approach improved all direct runoff estimates. In the NE catchment, better results were delivered for median. In catchments KO-2 and KO-3, only simulations with median λ applied in events separated according to P 10d led to increased accuracy in estimates. No positive response was recorded in simulations using the modified λ values in the KO-1 catchment (NSE ranged from −7.993 to 0.323).
In examining the calculations while not considering AMC classes (AMC implicitly was set to II, normal), higher accuracy in calculated runoff was recorded in the original CN approach, except for the KO-1 catchment. The cluster-based simulations led to better runoff matches only in the FU and NE catchments. In the other catchments, accuracy even decreased significantly.
In summary, the results show that replacing λ originally set to 0.2 with lower values can lead to improvements in direct-runoff simulations. However, the KO-1 catchment was the only one in which no positive effect from the λ modification was recorded. High e values indicated the model's over-estimating tendency in this specific catchment, probably due to improper CN determination. This might have been caused by the high percentage (97.3%) of arable land along with the presence of the tile drainage systems [78]. In other words, the results revealed quite specific hydrologic conditions in KO-1 catchment compared to the other two KO catchments. Moreover, this well illustrates that CN determination based solely on the prescribed criteria [65] might be misleading. Instead, it should take into consideration all available data on an area of interest.
It should be noted that some studies stated that the modification of λ should be accompanied by a corresponding adjustment in CN or S, separately [2,16]. Therefore, such an approach also was used in this research, but this did not lead to better results than those presented above. Table 6. Accuracy of direct runoff simulations using mean and median λ values in the experimental catchments with (antecedent moisture condition (AMC) all) and without (AMC = II) determination of AMC classes. NSE ≥ 0.500 indicating sufficient accuracy are highlighted.

Interpolated λ
Besides discrete values, λ can be expressed as a linear function of P [66,67]. The linear λ-P regressions were calculated for each catchment with and without AMC classification. Analogously to the previous approach (Section 4.2.1), regression equations were calculated for datasets delimited by the P 5d and P 10d thresholds determined through cluster analysis. All the regression equations are provided in Figure 4 and Table 7. Water 2020, 12, x FOR PEER REVIEW 16 of 28   Generally, application of the λ-P linear regressions increased the runoff estimates' accuracy. However, for higher Q the estimates showed systematically underestimating tendency (see Figure 5). Moreover, the effect differed between the experimental catchments. In the FU catchment, the best performance (NSE = 0.670) was recorded with using regression for the entire dataset (all-i.e., without the separation according to P 5d and P 10d ). In the other four catchments, the cluster-based regressions (i.e., those considering the separation according to P 5d and P 10d ) delivered the best outputs. However, only in catchments NE and KO-3 did NSE values exceed 0.500 (specifically, 0.670 and 0.540 in the NE catchment, 0.536 and 0.620 in the KO-3 catchment). Applying the all regression in the NE catchment reduced NSE from −0.070 (original CN method) to −0.740. In the KO-3 catchment, the effect was opposite as NSE increased from −26.3 to 0.386. In the remaining catchments KO-1 and KO-2 the results were similar, but NSE values did not indicate sufficient accuracy. Better results were obtained from simulations that omitted the AMC classification in general (see Figure S2), but NSE exceeded 0.500 only in the FU catchment (all regression-NSE = 0.600) and the KO-3 catchment (cluster-based P 5d regression-NSE = 0.507).
Compared with the "discrete λ" approach, the regression-based simulation appears to be more stable in terms of overall enhancement of the CN method's performance. The results indicate that λ should be a variable and event-dependent value. However, the irregular arrangement of the results hints that reliable usage of robust regression-based approach will require a subsequent detailed research of more rainfall-runoff data from various watersheds.
Water 2020, 12, x FOR PEER REVIEW 17 of 28 Figure 5. Comparison of observed direct runoff heights (mm) with those simulated in the experimental catchments using regressions of λ according to P for events with AMC classified (AMC all). Comparison of direct runoff estimates with AMC not classified (AMC = II) is provided in Figure  S2.

Event Analysis
Unlike the previous analyses, this approach does not start with determining CNs in the watersheds. Instead, it is based on analysing the rainfall-runoff events' parameters (P, Ia, Q), from which λ and S are derived. This procedure aims to analyse the array of λ and S to reveal a natural relation between them.
PCA revealed that P 10d poses the most important variable in the arrangement of events (see Figure 3). Therefore, regressions for S according to P 10d were calculated for each catchment. Mean and median λ identified in the training dataset were used (see Table 8). Unlike the tabulated CN-based approaches, the event analysis delivered more accurate direct runoff estimates (see Table 8, Figure 6). Again, the KO-1 catchment showed lower NSE values (0.536 and 0.511) than for the original CN method applied (0.861). On the other hand, the informative value of NSE in case of the traditional CN method with λ = 0.2 was addressed already (Section 4.2.1). Despite the overall improvement of the estimates, there is an underestimating tendency for higher Q evident, similar to the one seen in the tabulated CN-based approach.
Water 2020, 12, x FOR PEER REVIEW 19 of 28 sufficiently accurate hydrographs in terms of runoff volume/level, although the culmination is inaccurate. Figure 6. Comparison of observed direct runoff heights with those simulated in the experimental catchments using regressions of S according to P10d (event analysis) with application of mean and median λ. The results presented in Table 8 show that using P 10d -S regressions, along with modified λ, improved direct-runoff estimates' accuracy in all five studied catchments. NSE increased (for λ = 0.2, NSE ranged from −2.012 to 0.129) both with mean λ (NSE ranged from 0.159 to 0.683) and median λ (NSE ranged from 0.511 to 0.758). The observed-vs-modelled runoffs are plotted in Figure 6.

Model Fitting
The iterative model-fitting procedure was applied to training datasets in all studied catchments to identify pairs of S and λ values. As in previous analyses, computations were performed with and without the clustering (separation), according to the P 5d and P 10d thresholds. The S:λ pairs are provided in Table 9. Verification has shown that using the identified S:λ pairs leads to enhancement of estimated direct runoff, compared with those calculated through the original CN method. However, performance differs among catchments (Table 10, Figure 7). For simulations without the separation (all), NSE ranged from 0.146 to 0.585. NSE exceeded 0.500 in the FU and KO-3 catchments. For the P 5d -based separation, NSE ranged from 0.362 to 0.671, and NSE was higher than 0.500 in the FU and KO-3 catchments. Finally, for the P 10d -based separation, NSE ranged from 0.391 to 0.768, with NSE exceeding 0.500 in the NE and KO-3 catchments. In addition, here, the underestimating tendency is apparent for higher Q. Table 9. Pairs of S and λ values identified by the model fitting iterative procedure in the experimental catchments. Clustered (cl.) and non-clustered (all) events were used. The cluster delimitation was done using the P 5d and P 10d thresholds.   Figure 6. Comparison of observed direct runoff heights with those simulated in the experimental catchments using regressions of S according to P10d (event analysis) with application of mean and median λ.

Figure 7.
Comparison of observed direct runoff heights with those simulated using pairs of S and λ identified by means of model fitting procedure. The validation was carried out both with (all) and without delimitation of clusters according to P5d (clustered P5d) and P10d (clustered P10d).

Figure 7.
Comparison of observed direct runoff heights with those simulated using pairs of S and λ identified by means of model fitting procedure. The validation was carried out both with (all) and without delimitation of clusters according to P 5d (clustered P 5d ) and P 10d (clustered P 10d ).

HEC-HMS Simulations
To simulate a surface-runoff process using HEC-HMS, calibration was conducted for two types of events: (1) Those with a slowly rising limb of the hydrograph, and (2) those with a rapid onset. In the calibrations, all parameters' values (see Section 3.5 were optimised, and calibration coefficients were calculated for Ia, CN and T c ( Table 11). The remaining parameters' values were taken as they had been calculated in the calibrations. Table 11. Calibration coefficients for transformation of Ia, CN and T c in the three sub-basins of the NE catchment. The coefficients are provided for runoff events with both slow and rapid onset and for the λ = 0.2 and λ = 0.0142.

Slower Onset
Rapid Onset The HEC-HMS simulations indicate that changes in all surface-runoff-driving parameters led to improved results. For the event with the rising limb of the hydrograph rising rather slowly, significant improvement was achieved even for λ = 0.2 with modified parameters. In the original setting, the bias was −52.29% and NSE was −4.655 (see Figure 8a). In the simulation that employed the modified parameters, bias was only −3.17%, and NSE reached 0.918 (see Figure 8b). For the same event, a simulation using λ = 0.0142 with the original parameters' values delivered solid output, with bias at −7.01% and NSE reaching 0.845 (see Figure 8c). Finally, computations using λ = 0.0142 and the recalculated parameters led to improvement in the simulation, as bias fell to −2.12%, and NSE increased to 0.912 (see Figure 8d).   The simulation of runoff response with a rapid onset was less accurate, as indicated in Figure 9. The computation employing λ = 0.2 with the original setting of the surface-runoff-driving parameters was insufficient (see Figure 9a). After recalculating the parameters, accuracy increased (NSE = 0.408), but absolute bias worsened. The modelled runoff depicted in Figure 9b shows an approximation to the real data in terms of temporal course. Analogous use of λ = 0.0142 delivered a more pronounced improvement in the modelled runoff, with bias falling from 55% to 3%, and NSE increasing from −0.194 to 0.695 (Figure 9c,d).
Bearing in mind that the number of simulations by the HEC-HMS performed in this study is low, one should avoid general interpretations. Nevertheless, the presented results outline how reducing Ia improves the model's ability to deliver more accurate output, which is significantly more pronounced in the case of slow onset runoff, in which both the temporal course and runoff volume/level matched well with real runoff. In the case of rapid onset events, reducing λ from 0.2 to 0.0142 played a more significant role, as the best simulation was the one with λ = 0.0142 and recalculated parameters. However, the simulated hydrograph did not match well with the real results (Figure 9d). This indicates that the HEC-HMS simulation of such types of runoff events can deliver sufficiently accurate hydrographs in terms of runoff volume/level, although the culmination is inaccurate.

Discussion and Conclusions
Experience from various watersheds worldwide has highlighted the sensitivity of runoff to the initial abstraction coefficient [2,24,64]. In this study, initial abstraction coefficients λ were estimated for five experimental catchments in the Czech Republic, which well reflect the physiographic conditions (in terms of soil and geomorphology) in the Central European watersheds, using various approaches. From a methodological perspective, they can be divided into two types: (1) Approaches based on a classical determination of tabulated CN numbers [65], and (2) those based on a CN-independent analysis of rainfall-runoff data. The adjusted λ values were determined as discrete or function-dependent (interpolated). Contribution of this research lies in the comparison of different approaches of CN estimations applied in various watersheds.

Tabulated CN, Discrete λ
When using discrete λ values, these general conclusions based on our findings can be stated: λ differed both between studied catchments and between events, and the original value of 0.2 should be reduced in studied catchments, except the KO-1 catchment. Mean λ ranged from 0.0365 to 0.8330 and median λ from 0.0142 to 0.0913 for simulations that take the AMC classification into consideration. In simulations that omitted the AMC classification (AMC implicitly was set to II), λ values were slightly higher. Mean λ ranged from 0.0818 to 0.2846, and median λ ranged from 0.0479 to 0.2001; λ exceeded 0.2 in the KO-1 catchment, which is almost entirely covered by arable and land with specific infiltration conditions due to presence of tile drainage system [66], making this catchment very different from the others. The example of the KO-1 catchment well documents that correct CN determination must take into consideration artificial elements that are not reflected in the traditional methodology, yet their influence on runoff formation is significant.
Generally, the obtained results correspond with those of other authors, such as Baltas et al. [62], who reported mean λ = 0.034 and 0.014 for two watersheds. They attributed the various λ results to differences in the presence of impervious formations and urban development. Similarly, Shi et al. [11] reported λ ranging from 0.010 to 0.152, with a median of 0.048 due to landscape and geological characteristics. Lal et al. [13] conducted research in different agricultural plots in India and found that the mean and median for λ were 0.0 and 0.034. Krajewski et al. [24] compared λ values in 2 catchments in Poland and found average λ = 0.026 and 0.047.

Tabulated CN, Interpolated λ
Some approaches exist for expressing λ as being functionally (linearly) dependent on precipitation. For example, Kohnová et al. [69] conducted research in Slovakian watersheds and expressed λ as a linear function of precipitation for each AMC class. Moreover, they also discussed the low number of II and III AMC classes. Therefore, they stated that no apparent relationship existed between CN and AMC, and that other parameters' roles should be emphasised. Application of linear λ:P relationships in this study increased runoff estimates' overall accuracy compared with the original CN method. Like the discrete λ, NSE only exceeded 0.500 in one catchment-the KO-3 catchment with AMC classified, and in the KO-1 catchment without AMC classification. Generally speaking, the ability of λ:P linear regressions used for runoff simulations appears to be rather weak, which can be attributed to unknown interplay among other characteristics taking part in the direct-runoff formation, or to inappropriate tabulated CNs. Moreover, a question arises of the need of training datasets large enough, from which representative regression equations can be derived. Despite the sufficient amount of training events in this study (from 44 to 80), there is a problem of lack of extreme/severe events, in which the runoff formation may be very specific.
Recently, a different interpolating approach has been successfully tested by Krajewski et al. [24], approximating the λ-P relationship with an asymptotic formula. They found greater variability in λ values for smaller P. After exceeding a certain threshold P, λ approaches a constant value. The threshold depends on a catchment's characteristics, namely land cover that strongly influences the direct runoff formation or initial abstraction, separately. Based on the promising results of this asymptotic approach, the method definitely deserves to be tested in wider range of geographical conditions.

Summary of Tabulated CN-Based Approaches
Modification of λ in calculations using tabulated CNs led to a partial increase in the accuracy of direct-runoff estimates, but because NSE often did not exceed 0.500, overall performance seems to be rather unambiguous. More specifically, it is apparent that the estimates for smaller Q were more accurate than for those with higher Q. This issue might be resolved based on sophisticated statistical analyses that were beyond the scope of this study.
Both discussed approaches' inaccuracy indicates that the tabulated CNs [65] used in the Czech Republic often are not representative, thereby requiring a systematic revision that should be based on analyses of a large amount of data from various watersheds in the Central European region, ensuring the CN method's robustness. Moreover, just as it seems inappropriate to use the universal initial abstraction coefficient λ value for different river basins, it certainly would be more appropriate to use the CN method, taking into account the spatio-temporal variability of runoff conditions resulting from seasonal variability in CN and λ regarding vegetation variability and distribution of soil and land use [24,[30][31][32]. Extant studies have dealt with this issue, such as one conducted by Hawkins et al. [2], who concluded that lower CNs prevailed during the summer (growing season) and that higher CNs better represented runoff conditions during the dormant season. In a recent study by Muche et al. [79], CNs were derived using the Normalised Difference Vegetation Index (NDVI) to reflect seasonal land-use changes in the hydrologic circle. They concluded that the CN-NDVI method is a promising alternative to the traditional CN method. In the present study, the CORINE Land Cover dataset was refined using aerial images together with annual data on crop distribution on arable land, where available. It is reasonable to assume that by using appropriate satellite data, the estimated CNs would better represent the spatial distribution of land-cover features, vegetation and also soil moisture. However, this went beyond the scope of this research. Moreover, it should be noted that the Central European region currently is facing the bark beetle calamity and a drought affecting large areas of forest stands (mostly with spruces). This should be taken into account during determination of CNs, although the calamity did not affect the studied catchments. In this regard, it should also be added that the issue of a single CN value characterising the entire catchment's conditions has been addressed by Soulis and Valiantzas [30,31] who implemented a two-CN system approach in heterogeneous watersheds and improved the CN method's performance. Even in a watershed affected by massive changes caused by a forest fire, this approach performed well [32]. Thus, efficient methods exist, offering opportunity to be used in the process of new CN determination (revision), not only in the Czech Republic.

Tabulated CN-Independent Approaches
To estimate optimal λ and S values, two methods without using tabulated CNs were applied: (1) Event analysis based on analysing individual events' rainfall-runoff parameters, and (2) model fitting-an iterative numerical procedure searching for a pair of λ and S values that best describe a particular catchment's runoff conditions.
The event analysis yielded the best results from all the methods compared in this research. In all studied catchments, the application of S-P 10d regressions combined with mean (ranging from 0.0064 to 0.0129) and median λ (ranging from 0.0 to 0.0029) led to significant increases in estimated direct-runoff accuracy. With one exception (KO-3 catchment with mean λ applied), NSE always exceeded 0.500.
The model-fitting procedure also produced more accurate outputs compared with the original CN method, both with and without using the clustering (separation) according to the P 5d and P 10d thresholds. NSE exceeded 0.500 in only two out of the five studied catchments. This probably is caused by the fact that using a single pair of λ and S values that universally represent runoff conditions in a catchment in a wide range of rainfall and land-cover arrangements is questionable [24,30]. Another problem resides in the lack of events with high Q, of which the behaviour is strongly variable, and thus significantly influences the calculation of the λ:S pairs.

HEC-HMS Simulations Using Adjusted λ and CN
Application of adjusted λ values in HEC-HMS shows that using lower λ leads to improvements in simulated runoff accuracy, but a difference exists between hydrographs with slow and rapid onsets: Both yielded sufficiently accurate results in terms of runoff volume. However, within the hydrograph with rapid onset, the model was unable to match the culmination discharge. This might indicate a lower performance in the HEC-HMS SCS-CN-based simulations to model events with fast runoff response, especially flash floods, in terms of peak discharge, which often is the parameter of the main interest. At the same time, it should be emphasised that general interpretations would require an analysis on a representative dataset of rainfall-runoff events. In this present study, the aim was to examine preliminarily the applicability of modified λ values within HEC-HMS software.

Summary
The present study shows that the original CN methodology without adjustments to local conditions is hardly applicable in the geographical conditions of the experimental watersheds. Both CN and λ values should be based on site-specific characteristics, i.e., on rainfall-runoff data observations, hinting that they should be viewed as flexible regional parameters, rather than fixed values. In the studied catchments' specific natural conditions, the λ appears mostly to be considerably lower than the original value of 0.2. As the CN method frequently is used in long-term and event-based runoff simulations, erosion modelling and wide range of other applications, appropriate modifications of the method, such as those presented in this study, can contribute significantly to improvements in their outputs. It is the authors' intention to extent the present research and take up the unresolved problems, such as the systematic underestimation of direct runoff for events with higher runoff observed. Moreover, testing the possibilities of revisiting CN determination, mainly with regard to watersheds' heterogenity, as well as investigation of the regionalised concepts in the natural and agricultural conditions of the Czech Republic and neighbouring countries appears to be a challenging issue deserving our attention. This study's outcomes may serve as a preliminary basis for subsequent activities addressing the aforementioned tasks.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4441/12/7/1964/s1, Figure S1a: Comparison of observed direct runoff heights (mm) with those simulated in the experimental catchments using discrete mean and median λ values with AMC classified (AMC all), Figure S1b: Comparison of observed direct runoff heights (mm) with those simulated in the experimental catchments using discrete mean and median λ values with AMC implicitly set as normal (AMC = II), Figure S2: Comparison of observed direct runoff heights (mm) with those simulated in the experimental catchments using regressions of λ according to P for events with AMC classified (AMC all) and implicitly set as normal (AMC = II).