Derivation of Soil Moisture Recovery Relation Using Soil Conservation Service (SCS) Curve Number Method

: Soil moisture retention (SMR) capacity plays a key role in estimating the direct runoff when a multi-pulse storm event occurs. It is very important to know how much SMR will be recovered during the intervals of no rain of a multi-pulse storm. This study developed a new approach for derivation of the SMR recovery curve (R-curve) at sub-daily time-scales using the Curve Number (CN) method. The methodology was applied using complex storm events in the Napa River basin, California. The R-curve is classiﬁed into three sections depending on the recovery rate of SMR during the inter-storm interval of no rain (INR), and this study deﬁnes the characteristics. The ﬁrst section of the R-curve (INR 0–21 h with 0.97 mm/h) is described as gradually recovering SMR, since water is being inﬁltrated and the upper soil layer is not fully saturated. The second section (INR 21–36 h with 2.11 mm/h) is deﬁned as steeply recovering S due to downward drainage (sub-surface/inter ﬂows) and evaporation without inﬁltration. The third section (INR 36–68 h with 0.34 mm/h) is described as gradually decreasing recovery dependent on evaporation since percolation and drainage have almost stopped. (SMR) at hourly time-scales. Its application can improve continuous flood forecasting when a complex storm event occurs in a short time frame. The R-curve method was derived from a change of SMR and the interval of no rain (INR). The data were obtained for the upper Napa River basin in California and the sixteen rainfall events from 2011 to 2012 were used in estimating the event-based SMR. The methodology that is proposed in this study was based on the CN method, and is considered appropriate to apply for the other watersheds using observed direct runoff. The salient results of the study include:


Introduction
The state of soil moisture retention (SMR) plays a key role for estimating the antecedent moisture content (AMC) at storm initiation, which has a large effect on the overall runoff response [1][2][3][4][5][6]. SMR is the amount of water a soil column can retain. Continuously estimating SMR, however, is a challenge at sub-daily time resolution due to the recovery of SMR during non-rain intervals of multi-pulse storms [1,7]. The recovery is attributed to soil water drainage by subsurface flow, interflow, and evapotranspiration.
Various terms have been described in the literature focusing on the recovery of infiltration capacity to account for the water balance in soil [1,7,8]. Bauer [7] proposed a modification of the Horton's equation by introducing a soil drainage term into the equation to overcome a lack of the recovery of infiltration capacity during the interval of no rain (INR). Aron [9] suggested a modification of Soil Conservation Service Curve Number (CN) method by including a drainage rate decay factor [6]. Pitt et al. [8] conducted laboratory experiments on infiltration rates, which supported the previous two studies. Also, there are a variety of soil moisture accounting models considering percolation,

Curve Number Method and Soil Moisture Retention
The CN method is based on a water balance equation and two fundamental assumptions [6]. In Equation (1), the water balance statement equates total rainfall (P) to the sum of the amount of direct runoff (Q), the initial abstraction (I a ), and the amount of infiltration (F). The first assumption is the ratio of the amount of actual infiltration (F) to the amount of S (Equation (2)). The second assumption relates the initial abstraction (I a ) to S (Equation (3)). Ratio assumption: I a − S assumption: Here, λ is the initial abstraction coefficient, historically 0.2 [16]. Physically, this means that for a given storm, 20% of S is the initial abstraction before runoff begins [24].
The CN method calculates direct runoff from a rainfall depth, as shown in Equations (4)- (6). Q = (P − I a ) 2 P − I a + S for P > I a , and Q = 0 for P ≤ I a (4) Equation (4) becomes S is expressed in terms of the dimensionless CN taking values from 0, when S → ∞ , to 100, when S = 0. S = 1000 CN − 10 (6a) 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: In the CN method, S represents a storage capability of water infiltrated into soil [9]. Figure 1 shows a physical (a) soil profile, (b) soil-water-air layer, and (c) a conceptual diagram of soil moisture, accounting to understand the meaning of S using volumetric terms since S depends on the volumes of the soil, water, and air layer and soil textures in the soil profile.
Water 2018, 10, x FOR PEER REVIEW 3 of 22 Ratio assumption: a QF P I S   (2) a IS  assumption: Here,  is the initial abstraction coefficient, historically 0.2 [16]. Physically, this means that for a given storm, 20% of S is the initial abstraction before runoff begins [24].
The CN method calculates direct runoff from a rainfall depth, as shown in Equations (4)-(6). 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: In the CN method, S represents a storage capability of water infiltrated into soil [9]. Figure 1 shows a physical (a) soil profile, (b) soil-water-air layer, and (c) a conceptual diagram of soil moisture, accounting to understand the meaning of S using volumetric terms since S depends on the volumes of the soil, water, and air layer and soil textures in the soil profile.  In general, a soil column consists of surface soil, topsoil, subsoil, and bedrock, as shown in Figure 1a. Depending on the depth of each soil layer and their soil type, S of each soil layer may differ. However, for the CN method, S indicates the total capacity of SMR in a soil column because the method does not consider the soil profile with various soil layers.
In the soil column, infiltration and percolation rates are different since their mechanism differ, and they depend on AMC and air entrapment [25,26]. Also, the percolation rates from surface soil to topsoil and from topsoil to subsoil are different depending on the soil texture and the depth in each layer [27,28]. Generally, the upper layer has a loose texture, while the lower layer has a relatively lower hydraulic conductivity [29][30][31]. Therefore, a retardation effect on the boundary between the upper and lower layers could occur.
Each soil layer in the soil profile is classified among solids, water, and air, as shown in Figure 1b. In this figure, soil indicates the total volume of the soil particles, including its minerals and organics. The pore space of soil contains the liquid and gas phases of soil, i.e., everything but the solid phase that contains mainly minerals of varying sizes as well as organic compounds. Water means that the total water amount currently stored in the soil pore space, and represents AMC from previous precipitation [32]. Air can be part of the void space when not saturated with water. The volume of void space is determined by a degree of soil saturation under the assumption the total pore space volume of the soil column is constant. So, the void space is equal to the current potential maximum moisture retention. The sum of air and water are the absolute potential maximum moisture retention (S abs ) [14]. The total volume of voids and the absolute potential maximum moisture retention are theoretically equal, as shown by Equation (7).
In soil moisture accounting (Figure 1c), the sum of S and AMC is S abs , and the same as V v in the soil-water-air layer. An approach to estimating the current S uses a relationship between the absolute potential maximum retention and AMC (Equation (8)). According to this equation, S depends on AMC as S abs is constant. When the soil is very dry, AMC is reduced by drainage and evapotranspiration, and S approaches S abs . S = S abs − AMC (8) Using a time variable concept, Equation (8) becomes Equaiton (9), which represents the change of S after an INR. S(t + α) is the time variable S after the interval, α. However, S abs is a constant. AMC(t) is the soil moisture content at the end of antecedent rainfall and it depends on ds/d(t + α), indicating a change of S for INR. Also, it tracks the recovery of SMR. From the equation, it is evident that S depends on the recovery amount for the INR.

Methodology to Derive R-Curve
This study developed a methodology to derive the R-curve, which involved several preprocessing steps and a computation of S. Figure 2 shows the flowchart representing details of the workflow.
(a) First, we collected data for a number of complex storm events having multiple short-term independent rainfall pulses ( Figure 2(A-1)). This study employed a concept of Inter-Event Time Duration (IETD) in order to separate rainfall pulses from a complex rainfall event [33][34][35]. In this study, 6 h was used as the minimum IETD to classify an independent rainfall event. Other analysts have used 6 h as a proper IETD [1,[36][37][38][39]. Sixteen rainfall events were identified for the study period.
(b) Next, the complex hydrograph [40] from the multi-pulse rainfall event was separated into individual direct runoff hydrographs (A-2 in pre-processing in Figure 2). In the hydrograph separation process, the total direct runoff was calculated from the separated hydrograph to estimate an event-based CN for each rainfall pulse. This study used the N-day method [41] to identify the start of the event direct runoff from the complex hydrograph. Cubic spline interpolation was used to estimate the pulse event discharges on the regression limb. runoff from the complex hydrograph. Cubic spline interpolation was used to estimate the pulse event discharges on the regression limb. · Convert initial CN map to initial S (Si) map using the eq. (8b) · Using eq. (7) and initial S, estimate Sb1. At the eq. (7), gradually change the initial S by 1.0% until the best objective function is found, which is the smallest difference of Q and P in the eq. (8b). Also assume that the initial abstraction is 20% of initial S 2 1  (c) To estimate an event-based CN, an initial 1 km by 1 km CN map based on the land-use and soil maps was used (A-3 in pre-processing in Figure 2). To determine the appropriate event-based CN for each rainfall pulse, CN values in the CN map were changed uniformly and incrementally by 1.0% until a volume of simulated runoff flows by CN method matches a volume of the observed runoff flows. Areal average CN from the best appropriate event-based CN map was then used to calculate S as the representative value for the event. In addition, it should be emphasized that the initial abstraction is important to estimate CN values as it is related to the capacity of soil moisture retention. In this study for a developing purpose, it is assumed to be 20%. However, many studies have suggested various initial abstractions depending on region [6,23,42,43].
(d) Three S values for three time steps are required to derive the R-curve. The process computation of S in Figure 2B involves the definition of three S values at the specified times (Sb1, Sa1, and Sb2). Sb1, Sa1, and Sb2 are associated with the antecedent and post rainfall pulses. The first S value, Sb1, represents antecedent SMR for the antecedent rainfall pulse. Sb1 was calculated from the event-based CN in the process (c) and Equation (8b) (Figure 2(B-1)). The second S value, Sa1, represents SMR after the event rainfall pulse. Sa1 was calculated using the water balance Equation (3) and a relation of S, AMC, and infiltrated water (Figure 2(B-2)). In the computation process, infiltrated water was calculated from the (c) To estimate an event-based CN, an initial 1 km by 1 km CN map based on the land-use and soil maps was used (A-3 in pre-processing in Figure 2). To determine the appropriate event-based CN for each rainfall pulse, CN values in the CN map were changed uniformly and incrementally by 1.0% until a volume of simulated runoff flows by CN method matches a volume of the observed runoff flows. Areal average CN from the best appropriate event-based CN map was then used to calculate S as the representative value for the event. In addition, it should be emphasized that the initial abstraction is important to estimate CN values as it is related to the capacity of soil moisture retention. In this study for a developing purpose, it is assumed to be 20%. However, many studies have suggested various initial abstractions depending on region [6,23,42,43].
(d) Three S values for three time steps are required to derive the R-curve. The process computation of S in Figure 2B involves the definition of three S values at the specified times (S b1 , S a1 , and S b2 ). S b1 , S a1 , and S b2 are associated with the antecedent and post rainfall pulses. The first S value, S b1 , represents antecedent SMR for the antecedent rainfall pulse. S b1 was calculated from the event-based CN in the process (c) and Equation (8b) (Figure 2(B-1)). The second S value, S a1 , represents SMR after the event rainfall pulse. S a1 was calculated using the water balance Equation (3) and a relation of S, AMC, and infiltrated water (Figure 2(B-2)). In the computation process, infiltrated water was calculated from the water balance equation using direct runoff from a separated hydrograph, total rainfall, and 20% initial abstraction coefficient. The third S value, S b2 , represents SMR after the INR; it is the antecedent Water 2018, 10, 833 6 of 21 SMR for the post rainfall pulse. S b2 was estimated using the same process for estimating S b1 . Depending on the rainfall pulse, S b2 could be S b1 of the post rainfall pulse.
(e) The INR between two rainfall pulses was used as the time interval to recover S. In this study, the interval was assumed as the time period enabling S to recover through the downward movement of soil moisture and evapotranspiration. The INR was calculated while considering the time duration from the end of antecedent rainfall pulse to the beginning of the post rainfall pulse ( Figure 2).
(f) The R-curve was derived from the relationship between the change of S values and the INR between rainfall pulses. The change of S values was calculated as the difference between S a1 and S b2 . It indicates the change of S during the INR enabling S to recover.
Based on the methodology explained above, the relationship between the change of S and the duration as R-curve is derived for the sixteen rainfall events.

Watershed
The Napa River basin in Napa County, CA, was used as an application watershed. The gage at St. Helena (USGS 11456000), Napa County, CA, was the primary site to collect hydrographs as it is located in the upper portion of the basin and has minimal water management influences. Figure 3 shows the location of the Napa watershed and the land-use and the soil maps. The drainage area is 204.8 km 2 (79.3 mi 2 ). Annual average precipitation is 685 mm, of which 80% occurs mainly during the rainy season from November to March. water balance equation using direct runoff from a separated hydrograph, total rainfall, and 20% initial abstraction coefficient. The third S value, Sb2, represents SMR after the INR; it is the antecedent SMR for the post rainfall pulse. Sb2 was estimated using the same process for estimating Sb1. Depending on the rainfall pulse, Sb2 could be Sb1 of the post rainfall pulse.
(e) The INR between two rainfall pulses was used as the time interval to recover S. In this study, the interval was assumed as the time period enabling S to recover through the downward movement of soil moisture and evapotranspiration. The INR was calculated while considering the time duration from the end of antecedent rainfall pulse to the beginning of the post rainfall pulse ( Figure 2).
(f) The R-curve was derived from the relationship between the change of S values and the INR between rainfall pulses. The change of S values was calculated as the difference between Sa1 and Sb2. It indicates the change of S during the INR enabling S to recover.
Based on the methodology explained above, the relationship between the change of S and the duration as R-curve is derived for the sixteen rainfall events.

Watershed
The Napa River basin in Napa County, CA, was used as an application watershed. The gage at St. Helena (USGS 11456000), Napa County, CA, was the primary site to collect hydrographs as it is located in the upper portion of the basin and has minimal water management influences. Figure 3 shows the location of the Napa watershed and the land-use and the soil maps. The drainage area is 204.8 km 2 (79.3 mi 2 ). Annual average precipitation is 685 mm, of which 80% occurs mainly during the rainy season from November to March.  Table 1 lists the fractions (%) of land-use type and the hydrologic soil group in the basin ( Figure 3). The sum of forest and grassland/herbaceous areas exceeds 65%, followed by cultivated crops (9.9%), pasture/hay (8.4%), and shrub/scrub (7.1%). Based on the hydrologic soil group, the drainage capability is relatively low as the well-drained group A area is only 1.5%, and there are large areas of poorly drained group C (28.6%) and D (29.5%) soils.  Table 1 lists the fractions (%) of land-use type and the hydrologic soil group in the basin ( Figure 3). The sum of forest and grassland/herbaceous areas exceeds 65%, followed by cultivated crops (9.9%), pasture/hay (8.4%), and shrub/scrub (7.1%). Based on the hydrologic soil group, the drainage capability is relatively low as the well-drained group A area is only 1.5%, and there are large areas of poorly drained group C (28.6%) and D (29.5%) soils. The NRCS cover type influences the CN mapping since Min and Max values depend on the cover type. Table 1 shows Min and Max values that are associated with the NRCS cover types [16]. For example, the CN can range from 70 to 82 in group C of the forest category. Thus, CN maps resultant from the land-use and soil maps could differ depending on the analyst's choices. Also, it is important for estimating the event-based CN when the map is used as an initial condition.

Complex Storm Event
Complex storm events having many rainfall pulses in succession were used to derive the R-curve, with the S value of each pulse being linked to the previous and the post soil moisture contents. The rainfall pulses are established as independent of each other through hydrograph separation and a 6 h period as the minimum IETD. Figure 4 shows the time series of rainfall and streamflow for the four complex storm events that were disaggregated into 16 independent rainfall pulses. Table 2 shows the characteristics of the independent events (complex storm #01-three independent rainfalls; #02-seven independent rainfalls; #03-three independent rainfalls, #04-three independent rainfalls). In Table 2, a complex storm event consists of independent rainfalls with their total rainfall, the duration with no rain, and runoff coefficient, which is a dimensionless coefficient related to representation of a proportion of the amount of direct runoff to the amount of rainfall.

Radar-Based Rainfall Data
Multi-Radar Multi-Sensor (MRMS) was used as rainfall data in this study. MRMS is a National Weather Service (NWS) operational system having automated algorithms that integrate data streams from multiple radars, surface and upper air observations, lightning detection systems, and satellite and forecast models [44]. MRMS uses the Next-Generation Radar (NEXRAD) and rain gauge network data to produce a suite of Quantitative Precipitation Estimation (QPE) products, including radar-only, bias adjusted radar using rain gauges, and gauge-only. It generates gridded rainfall fields at 1 km resolution every hour. The MRMS archive can be downloaded for a user-selectable area at the public website (https://www.nssl.noaa.gov/projects/mrms/).
It is possible to estimate the CN value of each grid of the rainfall field since MRMS is a radar-based gridded rainfall product. Moreover, based on Willie et al. [45], the MRMS QPE was generated using the Vertical-Profile-of-Reflectivity (VPR) and the Mean Field Bias (MFB) correction method [46,47]. Figure 5 shows (a) a sample MRMS hourly rainfall map in Northern California and (b) the same sample MRMS hourly rainfall map for the application basin.

Complex Hydrograph Separation
Hydrograph separation of multiple peak flows arising from a multi-pulse rainfall event requires the calculation of direct runoff for each pulse. It should be emphasized that the complex hydrograph separation is important as the complex hydrograph separation methods are used to determine the total runoff volume. It, also, might have an effect on the R-curve derivation depending on the methods. To separate the multiple hydrograph into its components (i.e., base and direct flows) we used two methods, N-day and spline interpolation. The main reasons for selecting the two methods are as follows: (1) N-day method is convenient to estimate the end of a hydrograph using the drainage area and (2) the spline interpolation could reflect the trend of recession limb of hydrographs to a separated hydrograph. However, it is able to apply various hydrograph separation methods for deriving the R-curve [48]. When considering various methods, it is required to analyze the sensitivity of the hydrograph separation methods to examine their effect on the R-curve.
The N-day method estimates the duration between a time-to-peak and the runoff flow recession at the end of a hydrograph [41]; this is an adaptation of an equation like Equation (10)

Complex Hydrograph Separation
Hydrograph separation of multiple peak flows arising from a multi-pulse rainfall event requires the calculation of direct runoff for each pulse. It should be emphasized that the complex hydrograph separation is important as the complex hydrograph separation methods are used to determine the total Water 2018, 10, 833 9 of 21 runoff volume. It, also, might have an effect on the R-curve derivation depending on the methods. To separate the multiple hydrograph into its components (i.e., base and direct flows) we used two methods, N-day and spline interpolation. The main reasons for selecting the two methods are as follows: (1) N-day method is convenient to estimate the end of a hydrograph using the drainage area and (2) the spline interpolation could reflect the trend of recession limb of hydrographs to a separated hydrograph. However, it is able to apply various hydrograph separation methods for deriving the R-curve [48]. When considering various methods, it is required to analyze the sensitivity of the hydrograph separation methods to examine their effect on the R-curve.
The N-day method estimates the duration between a time-to-peak and the runoff flow recession at the end of a hydrograph [41]; this is an adaptation of an equation like Equation (10). For the Napa River near Helena site N-days is approximately 2.4 days (58 h), when substituting the drainage area 205.4 km 2 .
N days = (0.386 × A) 0.2 (10) Here, N is the number of days runoff ceases after a storm and A is the drainage area of the application basin (km 2 ). Figure 6 shows the process to separate a multiple peak hydrograph into its pulse components, described as follows: (1) Estimate the duration from a time to peak to the end time of runoff using the N-day method, as shown in Figure 6A. (2) Use the spline interpolation method to estimate the pulse discharges on the regression limb, as shown in Figure 6B. When applying the interpolation method, this study used the discharges for the recession limb and the end time of runoff flow from the N-day method. Finally, the separated hydrograph is shown in Figure 6C. A similar process is done for the remaining peaks.

Estimation of Soil Moisture Retention
A gridded CN for each storm event was estimated using an initial gridded CN map, as described at Section 3. The initial gridded CN map was estimated from Table 2 with both minimum and maximum CN values to determine the difference of CNs depending on NRCS cover types. Figure 7 shows the result of the initial gridded CN map and its corresponding S values. As shown in Figure 7, the areal mean CN using the minimum values is 67; the maximum mean CN value is 81. In Figure 7b, S and CN are inversely proportional, so the maximum S is from the minimum CN value. The maximum S is twice the minimum S, because the CN-S relationship is non-linear (Equation (6)). Areal mean S using minimum CN is 134.3 mm and the areal mean S using maximum CN is 61.4 mm.

Estimation of Soil Moisture Retention
A gridded CN for each storm event was estimated using an initial gridded CN map, as described at Section 3. The initial gridded CN map was estimated from Table 2 with both minimum and maximum CN values to determine the difference of CNs depending on NRCS cover types. Figure 7 shows the result of the initial gridded CN map and its corresponding S values. As shown in Figure 7, the areal mean CN using the minimum values is 67; the maximum mean CN value is 81.

Estimation of Soil Moisture Retention
A gridded CN for each storm event was estimated using an initial gridded CN map, as described at Section 3. The initial gridded CN map was estimated from Table 2 with both minimum and maximum CN values to determine the difference of CNs depending on NRCS cover types. Figure 7 shows the result of the initial gridded CN map and its corresponding S values. As shown in Figure 7, the areal mean CN using the minimum values is 67; the maximum mean CN value is 81. In Figure 7b, S and CN are inversely proportional, so the maximum S is from the minimum CN value. The maximum S is twice the minimum S, because the CN-S relationship is non-linear (Equation (6)). Areal mean S using minimum CN is 134.3 mm and the areal mean S using maximum CN is 61.4 mm. These S values indicate the amount of SMC the soil column could hold at any time.
The minimum value in Table 1 was used to map the initial CN since it is better than the maximum value to estimate event-based CN. The non-proper values resulted in a large difference of the observed direct runoff and the calculated direct runoff from the estimated CN by the maximum value. Results from the minimum value have relatively small difference.
The resulting S values (Sb1, Sa1, Sb2, and a change of S) are tabulated in Table 3. This table also lists the event-based CNs and INR for deriving the R-curve for each event. In Table 3, INR is calculated from the time interval between antecedent rainfall pulse and post rainfall pulse with no rain, the recovery of S is a recovered capacity of SMR for the duration. For the independent rainfall events No. 2 and 3 in Table 3, the recovery of SMR, 8.0 mm, occurred over a duration of 29 h. The first S value, Sb1, is estimated at the time just before the onset of the first rainfall. The second S value, Sa1, is representative of the time just after the first rainfall period ends. The third S value, Sb2, represents the first response of the second hydrograph arising from the second rainfall period. In Figure 7b, S and CN are inversely proportional, so the maximum S is from the minimum CN value. The maximum S is twice the minimum S, because the CN-S relationship is non-linear (Equation (6)). Areal mean S using minimum CN is 134.3 mm and the areal mean S using maximum CN is 61.4 mm. These S values indicate the amount of SMC the soil column could hold at any time.
The minimum value in Table 1 was used to map the initial CN since it is better than the maximum value to estimate event-based CN. The non-proper values resulted in a large difference of the observed direct runoff and the calculated direct runoff from the estimated CN by the maximum value. Results from the minimum value have relatively small difference.
The resulting S values (S b1 , S a1 , S b2 , and a change of S) are tabulated in Table 3. This table also lists the event-based CNs and INR for deriving the R-curve for each event. In Table 3, INR is calculated from the time interval between antecedent rainfall pulse and post rainfall pulse with no rain, the recovery of S is a recovered capacity of SMR for the duration. For the independent rainfall events No. 2 and 3 in Table 3, the recovery of SMR, 8.0 mm, occurred over a duration of 29 h. The first S value, S b1 , is estimated at the time just before the onset of the first rainfall. The second S value, S a1 , is representative of the time just after the first rainfall period ends. The third S value, S b2 , represents the first response of the second hydrograph arising from the second rainfall period. Table 3. Results of three S estimated and the recovery amount of S for the duration with no rain using the proposed methodology.

Derivation of the Recovery Curve
The R-curve is derived from the relationship between the change of S values (the recovery in Table 3) and INR, as shown in Figure 8a. The R-curve is divided into three sections, according to the recovery rate (mm/h) of SMR. Root Mean Square Error (RMSE) is used as a criterion to determine the best linear regression line to estimate the recovery rate in each section. For the three sections, (1) the first section is in the range 10 to 21 h and exhibits a gradually recovering S, (2) the second section is in the range 21 to 36 h and exhibits a steeply recovering S, and (3) the third section is in the range 36 to 66 h and it exhibits a gradually decreasing recovery rate. However, ΔS is negative because a1 S at the end of first storm event is larger than b1 S after an INR.
From the Napa River events, the negative S values continue to 27 h, which means that soil moisture contents are being increased without additional rainfall for 27 h. In other words, actual SMR is decreased instead of increasing. In general, SMR is decreased in the short-term by drainage arising from the movement of water flowing to lower soil layers and this process is commonly finished within 2-3 days. This physical process is difficult to describe using the CN method because it regards a soil column as a single layer in accounting for the total value of SMR. The results that are described above can be explained conceptually when dividing a soil column However, ∆S is negative because S a1 at the end of first storm event is larger than S b1 after an INR. From the Napa River events, the negative S values continue to 27 h, which means that soil moisture contents are being increased without additional rainfall for 27 h. In other words, actual SMR is decreased instead of increasing. In general, SMR is decreased in the short-term by drainage arising from the movement of water flowing to lower soil layers and this process is commonly finished within 2-3 days. This physical process is difficult to describe using the CN method because it regards a soil column as a single layer in accounting for the total value of SMR.
The results that are described above can be explained conceptually when dividing a soil column into two or more layers reflecting the soil profile and taking consideration of physical phenomenon occurring within each layer. Figure 9 shows a conceptual diagram of a soil column consisting of two soil layers with S and SMC, and showing infiltration, percolation, and drainage by arrows. In this study, the upper layer was defined as a soil layer that is saturated by infiltration, and the lower layer is defined as a soil layer saturated by downward percolation. Wetting front Figure 9. Conceptual diagram of a soil column considering the soil layers in the first section: This study assumed that the soil column consists of two layers, an upper layer (saturated by infiltration) and a lower layer (saturated by percolation). S indicates the soil moisture retention (SMR) and SMC means Soil Water Content consisting of antecedent moisture content (AMC) and the infiltrated water. The wetting front based on Green and Ampt [49] represents the boundary between two layers causing the retardation effect by the difference of infiltration and percolation rates.
According to Thomas et al. [50], there is a gap between the infiltration rate (92-267 mm/day) and the percolation rate (40-109 mm/day). If percolation is slower than the infiltration rate, a retardation boundary between two layers occurs because water is not readily flowing to the lower layer, as shown in Figure 9. The boundary could be described as the wetting front. Thus, when the upper layer is saturated and there is the retardation effect, it appears as if the total SMR is decreased without additional rainfall during the INR, even though the actual total SMR is increased by the drainage.
However, it is hard to conclude that a soil column has only one wetting front since the column consists of numerous soil layers. Percolation at the layers below the upper layer depends on the thickness and saturated hydraulic conductivity of each layer [29,30]. Generally, the upper layer has a loose texture, while the other layer below the upper layer has a relatively lower hydraulic conductivity [31]. Thus, it is highly possible that a soil column has several wetting fronts.
As mentioned above, this result is a limitation of the CN method since the method assumes that entire soil column is a single layer. To overcome this limitation and to derive the appropriate R-curve, a correction is needed.
This study examined the effects of the complex hydrograph separation and the initial abstraction (Ia) on the derivation of the R-curve. Figure 8b shows the derived R-curve results from various N-days and Ia. For applying the various N-days, this study changed ±10% of the estimated N-day, 58 h for the application watershed. Two Ia, 0.05 and 0.10, were used for the effect analysis. As a result, the changes of the N-days and the Ia have no obvious effect on the derivation of the R-curve. In the three sections of the Figure 9. Conceptual diagram of a soil column considering the soil layers in the first section: This study assumed that the soil column consists of two layers, an upper layer (saturated by infiltration) and a lower layer (saturated by percolation). S indicates the soil moisture retention (SMR) and SMC means Soil Water Content consisting of antecedent moisture content (AMC) and the infiltrated water. The wetting front based on Green and Ampt [49] represents the boundary between two layers causing the retardation effect by the difference of infiltration and percolation rates.
According to Thomas et al. [50], there is a gap between the infiltration rate (92-267 mm/day) and the percolation rate (40-109 mm/day). If percolation is slower than the infiltration rate, a retardation boundary between two layers occurs because water is not readily flowing to the lower layer, as shown in Figure 9. The boundary could be described as the wetting front. Thus, when the upper layer is saturated and there is the retardation effect, it appears as if the total SMR is decreased without additional rainfall during the INR, even though the actual total SMR is increased by the drainage.
However, it is hard to conclude that a soil column has only one wetting front since the column consists of numerous soil layers. Percolation at the layers below the upper layer depends on the thickness and saturated hydraulic conductivity of each layer [29,30]. Generally, the upper layer has a loose texture, while the other layer below the upper layer has a relatively lower hydraulic conductivity [31]. Thus, it is highly possible that a soil column has several wetting fronts.
As mentioned above, this result is a limitation of the CN method since the method assumes that entire soil column is a single layer. To overcome this limitation and to derive the appropriate R-curve, a correction is needed.
This study examined the effects of the complex hydrograph separation and the initial abstraction (I a ) on the derivation of the R-curve. Figure 8b shows the derived R-curve results from various N-days and I a . For applying the various N-days, this study changed ±10% of the estimated N-day, 58 h for the application watershed. Two I a , 0.05 and 0.10, were used for the effect analysis. As a result, the changes of the N-days and the I a have no obvious effect on the derivation of the R-curve. In the three sections of the R-curve, the recovery rates of R-curves with various conditions have only a little change under 5.0% against to the R-curve derived from the N-day with no change and I a of 0.20. The result is natural because the R-curve is based on the relation of INR and a change of the Ss. Since INRs are constant for all cases, a change of the Ss is important. A change of the Ss is defined as the difference between S b2 and S a1 in this study. When changing the N-days or the I a , S b2 and S a1 are changed together as the total runoff flows are increased or decreased. For example, the independent rainfall No. 1 in Table 3 has 1.0 mm of the recovery difference between the results from 0.20 and 0.05 I a . It is same when changing the N-days.
In this study, the R-curve was corrected assuming two conditions: (1) a trend of SMR between 0 to 10 INR follows the trend between 10 to 21 INR. Based on this assumption, the regression line in the first section is extended to the time zero, and (2) the difference between rates of infiltration and percolation, which cause the retardation effect, is constant from time zero. Given these assumptions, then the R-curve can be adjusted until ∆S is zero at the time zero. Figure 10 shows the result of the adjusted R-curve. In the Napa River case, the correction was 33.6 mm. which cause the retardation effect, is constant from time zero. Given these assumptions, then the R-curve can be adjusted until ΔS is zero at the time zero. Figure 10 shows the result of the adjusted R-curve. In the Napa River case, the correction was 33.6 mm. Figure 10. Diagram of the R-curve adjustment process: In this figure, the arrows represent the adjustment process of the uncorrected R-curve. Figure 11 shows the derived R-curve consisting of three sections with a hydrograph for event No. 16. The R-curve, recovered SMR, and CN curves are plotted on the hydrograph to show a potential application instance for flood forecasting. The first section for low INR periods is characterized by gradually recovering SMR values. For the Napa River watershed, the time period of this section is 0-21 h with 0.97 mm/h (23.28 mm/day). In this section, surface water is being infiltrated into soil, so the rate is slower than for the second section. Also, downward percolation is a primary contributor to the recovery of SMR in a short-term (Figure 12a).  Figure 11 shows the derived R-curve consisting of three sections with a hydrograph for event No. 16. The R-curve, recovered SMR, and CN curves are plotted on the hydrograph to show a potential application instance for flood forecasting. The first section for low INR periods is characterized by gradually recovering SMR values. For the Napa River watershed, the time period of this section is 0-21 h with 0.97 mm/h (23.28 mm/day). In this section, surface water is being infiltrated into soil, so the rate is slower than for the second section. Also, downward percolation is a primary contributor to the recovery of SMR in a short-term (Figure 12a). The second section of the R-curve involves more rapid recovery of SMR. For the Napa River events, the time period of this section is 21-36 h (Figure 11). The recovery rate is 2.11 mm/h (50.64 mm/day), and is the fastest rate among three sections of the R-curve (the rate is about two times faster than the first section and over seven times faster than the third section). From Figure 12, the main reasons for the faster recovery rate in Section 2 are: (1) infiltration is stopped, (2) the downward percolation rate in the second section is faster than the first section since the lower layer is more saturated than the first layer, and (3) the downward percolation continues with time.
The third R-curve section involves a gradually decreasing recovery rate, as shown in Figure 11. In this section, the downward movement of water is provisionally completed, and evapotranspiration  The second section of the R-curve involves more rapid recovery of SMR. For the Napa River events, the time period of this section is 21-36 h (Figure 11). The recovery rate is 2.11 mm/h (50.64 mm/day), and is the fastest rate among three sections of the R-curve (the rate is about two times faster than the first section and over seven times faster than the third section). From Figure 12, the main reasons for the faster recovery rate in Section 2 are: (1) infiltration is stopped, (2) the downward percolation rate in the second section is faster than the first section since the lower layer is more saturated than the first layer, and (3) the downward percolation continues with time.
The third R-curve section involves a gradually decreasing recovery rate, as shown in Figure 11. In this section, the downward movement of water is provisionally completed, and evapotranspiration The second section of the R-curve involves more rapid recovery of SMR. For the Napa River events, the time period of this section is 21-36 h (Figure 11). The recovery rate is 2.11 mm/h (50.64 mm/day), and is the fastest rate among three sections of the R-curve (the rate is about two times faster than the first section and over seven times faster than the third section). From Figure 12, the main reasons for the faster recovery rate in Section 2 are: (1) infiltration is stopped, (2) the downward percolation rate in the second section is faster than the first section since the lower layer is more saturated than the first layer, and (3) the downward percolation continues with time.
The third R-curve section involves a gradually decreasing recovery rate, as shown in Figure 11. In this section, the downward movement of water is provisionally completed, and evapotranspiration might be a dominant factor to recover SMR (Figure 12c). Thus, the recovery rate in the third section is 0.34 mm/h (8.16 mm/day), which is the slowest among the three sections. However, total capacity of SMR in this section is the largest of the three sections as the entire soil column is involved. The third section is initiated at 36 h and ends around 68 h. Ultimately, the recovery rate of SMR after the completion of third section is going to converge on a constant value that is close to zero because the evaporation approaches zero. Also, the remaining water in soil could be decreased by transpiration in the long-term [51][52][53].

Application Instance
A rainfall-runoff simulation was implemented as an example to demonstrate the usefulness of the new R-curve methodology for continuous runoff simulation. The HEC-HMS model developed by the US Army Corps of Engineers is used as the rainfall-runoff model. Among many of the rainfall-runoff methods in HEC-HMS, the ModClark approach that is based on Clark's unit hydrograph [54] was selected as a transform method for runoff routing. There are two reasons for selecting the ModClark method: (1) The model is a semi-distributed hydrologic model which can use gridded rainfall inputs such as a radar-based rainfall data; and, (2) It is convenient to apply the gridded CN method as the ModClark is based on a grid representation. For the detailed methodology of the ModClark, readers can see references [55][56][57]. Also, many previous studies have used the ModClark for various hydrologic applications [58][59][60][61][62][63][64]. In most of the studies cited above, the CN method was used with the ModClark approach.
The ModClark method requires two parameter inputs, time of concentration (Tc) defined as a time needed for runoff flow from the most remote location to the outlet in a watershed and storage coefficient (K) defined as a lag time due to storage effects. Using independent storm events in Table 2, the parameters are estimated. This study employs the HEC-HMS optimization tools [65] to estimate the proper parameters for each storm event. Table 4 shows the estimated parameters with error indices (peak flow, total runoff flow, correlation coefficient (CC), bias (BS), and Nash-Sutcliffe efficiency coefficient (NSE)). Except for two events (1 and 4), all of peak flow error (%) and total runoff flow error (%) were under 25%. The mean error values of peak flow, total runoff flow, and time to peak are 6.4%, −0.2%, and −0.1 h, respectively. CC, BS, and NSE, also, confirm that the simulation results using the estimated parameters are acceptable. For the application instance, the mean values from the estimated parameters are used as the representative parameters (Tc = 8.5 h and K = 6.5 h). Figure 13 shows a scatterplot of the observed runoff flows and the calibrated runoff flows using the representative parameters. The data are for the same sixteen storm events in Table 4. The error indices, NSE, BS, and CC, showed similar results, as above.

Figure 13
Scatterplot of the observed and calibrated runoff flows for the sixteen storm events in Table 4: In the figure, the acronyms are same as Table 4 and the values are in average.
In the application example, a different complex storm event was simulated, which was not used in the derivation of the R-curve and the estimation of the model parameters. The complex storm event occurred from 01-December-2014 to 11-December-2014. During the period, three independent rainfall events generated three corresponding streamflow responses. The first event is used as a spin-up period and as a starting point of estimating direct runoff using the R-curve. The MRMS data was used as the Figure 13. Scatterplot of the observed and calibrated runoff flows for the sixteen storm events in Table 4: In the figure, the acronyms are same as Table 4 and the values are in average.
In the application example, a different complex storm event was simulated, which was not used in the derivation of the R-curve and the estimation of the model parameters. The complex storm event occurred from 01-December-2014 to 11-December-2014. During the period, three independent rainfall events generated three corresponding streamflow responses. The first event is used as a spin-up period and as a starting point of estimating direct runoff using the R-curve. The MRMS data was used as the gridded rainfall data. An hourly time step was chosen as the temporal resolution of the model simulation. To apply the R-curve, a period for the model spin-up was needed. This process is common in continuous hydrologic models to set up initial conditions. During this spin-up period, a CN value was the initial condition for the soil moisture state, and the R-curve can start to estimate soil moisture retention with S values at the end of the spin-up period. In this application instance, the spin-up period encompasses the rainfall event from 01-December-2014 to 03-December-2014, and the two rainfall events on 06-December-2014 and 11-December-2014 were used to verify the R-curve performance with the hydrologic model. Figure 14 shows the rainfall-runoff simulation results with streamflow that was observed by the USGS stream station. In Figure 14, there are two simulations. One is the result of using the R-curve method, and the other one (antecedent condition) represents the conventional CN method by using AMC conditions [16]. The AMC conditions (AMC I, II, and III) depend on the accumulated precipitation amount over a five-day antecedent period (P5). For the 06-December-2014 event, the results show that the difference of total runoff flows between the observed and R-curve is only 11%. For the 11-December-2014 event, the results show that the R-curve method also performed well with errors in peak flow and total runoff flow of 23% and −6%, respectively. On the other hand, the simulation from the conventional CN method showed that the errors of peak flow and total runoff flow are −31% and −42%, respectively. According to the conventional CN method criteria used to decide the AMC conditions (I, II, and III), the 11-December-2014 was classified as AMC I (P5 < 35.5 mm in the rainy season [16]), which is a dry condition. This classification resulted in a CN value for the 11-December-2014 event that was underestimated. gridded rainfall data. An hourly time step was chosen as the temporal resolution of the model simulation. To apply the R-curve, a period for the model spin-up was needed. This process is common in continuous hydrologic models to set up initial conditions. During this spin-up period, a CN value was the initial condition for the soil moisture state, and the R-curve can start to estimate soil moisture retention with S values at the end of the spin-up period. In this application instance, the spin-up period encompasses the rainfall event from 01-December-2014 to 03-December-2014, and the two rainfall events on 06-December-2014 and 11-December-2014 were used to verify the R-curve performance with the hydrologic model. Figure 14 shows the rainfall-runoff simulation results with streamflow that was observed by the USGS stream station. In Figure 14, there are two simulations. One is the result of using the R-curve method, and the other one (antecedent condition) represents the conventional CN method by using AMC conditions [16]. The AMC conditions (AMC I, II, and III) depend on the accumulated precipitation amount over a five-day antecedent period (P5). For the 06-December-2014 event, the results show that the difference of total runoff flows between the observed and R-curve is only 11%. For the 11-December-2014 event, the results show that the R-curve method also performed well with errors in peak flow and total runoff flow of 23% and −6%, respectively. On the other hand, the simulation from the conventional CN method showed that the errors of peak flow and total runoff flow are −31% and −42%, respectively. According to the conventional CN method criteria used to decide the AMC conditions (I, II, and III), the 11-December-2014 was classified as AMC I (P5 < 35.5 mm in the rainy season [16]), which is a dry condition. This classification resulted in a CN value for the 11-December-2014 event that was underestimated.

Discussion and Conclusions
This study developed a new approach for deriving a soil moisture retention recovery curve (R-curve). The R-curve provides a means for tracking the time-variable state of soil moisture retention (SMR) at hourly time-scales. Its application can improve continuous flood forecasting when a complex storm event occurs in a short time frame. The R-curve method was derived from a change of SMR and the interval of no rain (INR). The data were obtained for the upper Napa River basin in California and the sixteen rainfall events from 2011 to 2012 were used in estimating the event-based SMR. The methodology that is proposed in this study was based on the CN method, and is considered appropriate to apply for the other watersheds using observed direct runoff. The salient results of the study include:

Discussion and Conclusions
This study developed a new approach for deriving a soil moisture retention recovery curve (R-curve). The R-curve provides a means for tracking the time-variable state of soil moisture retention (SMR) at hourly time-scales. Its application can improve continuous flood forecasting when a complex storm event occurs in a short time frame. The R-curve method was derived from a change of SMR and the interval of no rain (INR). The data were obtained for the upper Napa River basin in California and the sixteen rainfall events from 2011 to 2012 were used in estimating the event-based SMR. The methodology that is proposed in this study was based on the CN method, and is considered appropriate to apply for the other watersheds using observed direct runoff. The salient results of the study include: We addressed the importance of considering the soil profile with multi-layers since interaction between soil layers by infiltration and percolation influences the change of SMR. In this study, the recovery of SMR had a negative value at the beginning of INR since the soil column is assumed to be a single layer in the CN method. Thus, SMR resulting from the CN method cannot account for the negative retardation effect arising from the difference of downward flows by infiltration and percolation. However, the retardation effect concept is theoretically the same as the Green and Ampt infiltration model due to the wetting front. When the infiltration rate is faster than the percolation rate, the boundary (wetting front) between the upper layer and the lower layer could develop a temporary stagnant water layer when the retardation effect occurs. Thus, when the upper layer is saturated, but not the lower layer, it appears as if the total SMR is decreased without additional rainfall for the INR, even though the actual total SMR is recovered.
For the three sections of the R-curve, the first section is regarded as the infiltration period, so the recovery rate is relatively slower than the second section. Also, the downward percolation is the main contributor to recover SMR in this time period. The second R-curve section involves the rapid recovery of SMR. It is the fastest rate among three sections since most of soil layer is almost saturated. In the third section, the movement of water is provisionally complete, and evaporation is a potential dominant factor to the recovery of SMR. Thus, the recovery rate is the slowest among three sections. The recovery rate of this section converges on a constant value that is close to zero in the long-term.
The HEC-HMS model was used to comparison the performance of the CN and R-curve methods for a storm event with multiple rain and streamflow response pulses. It was shown that the R-curve method outperformed the traditional CN method in terms of characterizing the peak flow and the total runoff volume.
Based on the results above, this new R-curve approach could be used to improve the limitation of the traditional CN methodology, which regards the soil column as a single layer in accounting for the total value of SMR. Ultimately, the R-curve method could be used to estimate the state of SMR, which would improve flood forecasting during a complex storm event with multiple rainfall pulses without using complicated hydrologic models. Since the R-curve method can indicate the state of soil moisture using the recovery rates, it could be used to represent a critical time period that is vulnerable to flooding due to saturated soils, even by light rainfall.