The runoff coefficient of a hillslope is a reliable measure for changes in the streamflow response at the river link outlet. A high runoff coefficient is a good indicator of the possibility of flash floods. Although the relationship between runoff coefficient and streamflow has been the subject of much study, the physical mechanisms affecting runoff coefficient including the dependence on precipitation pattern remain open topics for investigation. In this paper, we analyze a rainfall-runoff model at the hillslope scale as that hillslope is forced with different rain patterns: constant rain and fluctuating rain with different frequencies and amplitudes. When an oscillatory precipitation pattern is applied, although the same amount of water may enter the system, its response (measured by the runoff coefficient) will be maximum for a certain frequency of precipitation. The significant increase in runoff coefficient after a certain pattern of rainfall can be a potential explanation for the conditions preceding flash-floods.

During past three decades, significant advances have been made in understanding the correlation between antecedent soil moisture and runoff coefficient for simulated rain events, with the goal of predicting floods. Several methods have been proposed by which water movement through the soil causes changes in the runoff coefficient and the channel runoff. Nevertheless, despite this concentrated effort, the underlying mechanisms of the inter-dependence between soil dynamics and rain storm characteristics remain insufficiently explained. In this paper, we analyze the nonlinear dynamics of a recently proposed hydrological model at the hillslope-scale, under fluctuating precipitation input. We show that, under certain conditions that depend both on the hillslope's properties and past rain intensity and frequency, the runoff coefficient can reach an amplified peak. This finding suggests a possible explanation of how the structure of the rainfall, even though in average, same precipitation volume is applied, is important for the setting of the hillslope state to conditions optimal to produce flash-flooding.

Although floods typically occur following high rates of precipitation, the reverse is not necessarily true. In fact, the same amount of precipitation when applied in different patterns can lead to significantly different hydrographs at the hillslope scale. The possibility then exists of the same volume of precipitation causing flooding in one situation but not another.

Physically, the method by which the hillslope maintains information from previous storm events is through the water that remains in the soil. Antecedent soil moisture represents a state of the hillslope system and is well known to affect surface runoff by influencing the amount of infiltration that can occur during a storm. Soil moisture has been identified to be a more dominant factor in determining surface runoff for a rain event at the hillslope scale than event rainfall depth.13 Several studies have investigated the contribution from soil moisture on runoff coefficient under semi-arid climate conditions for various rainfall events.1,5,16 Considering hillslopes with and without plant cover, Castillo et al.1 discovered that, for semi-arid climates, antecedent soil moisture is a controlling factor influencing soil runoff specifically when the rain pattern has low to medium intensity but not for high intensity, low frequency storms. The sensitivity of the soil-moisture influence depends upon the dominant runoff mechanisms in these situations. Fitzjohn et al.5 investigated the impacts of spatial variability in soil moisture in the semi-arid climate and focused on the effects of soil moisture on runoff and erosion. To contribute to predictability, Zhang et al.16 investigated the relationship between the soil moisture and its predictability on the influence of soil moisture on runoff. Continued effort is being applied to make appropriate estimations of soil moisture conditions in locations without measurements with the acceptance that soil moisture can help to predict flash floods.7,15 Much of this previous effort has focused on the correlation between antecedent soil moisture and runoff coefficient for simulated rain events to be used for prediction. The study involved Monte Carlo methods and statistical sensitivity analyses of current models. Although the authors suggest methods by which water moves through the soil that can cause this influence to be so strong, the underlying mechanisms in the soil that influence this relationship need to be more thoroughly analyzed.

The purpose of this paper is to investigate how fluctuating patterns of precipitation affect the hillslope's response, potentially leading to pre-flooding conditions. We are particularly interested in determining sufficient conditions on the hillslope's properties and rainfall structure (intensity and frequency) for large variations (increases) in the runoff coefficient to occur, in a relatively short period of time, as usually observed in flash floods. The runoff coefficient is defined as the ratio of surface runoff to precipitation effectively entering the system (water infiltrating into the soil and water being transported to the channel). It is an important hydrological measure for determining the amount of water that exits the hillslope as surface runoff.10 It takes subunitary values (RC ∈ [0,1]) with the following interpretations: a low value of the runoff coefficient (RC ≈ 0) corresponds to a hillslope with dry, permeable soil so that most rainfall will be infiltrated into the soil; contrarily, a large value of the runoff coefficient (RC ≈ 1) is associated with saturated soil or an impermeable surface that cannot absorb any water so that precipitation quickly flows to the river as surface runoff. The latter is typically reached during floods.

We propose to address this question in the context of a rainfall-runoff model recently proposed in Ref. 4. (See also Refs. 2, 3, and 9 for similar modeling approaches.) There are several advantages to this approach. First, the model in Ref. 4 is an integral-balance nonlinear system that describes the simultaneous dynamics of the soil moisture, ground-water, and surface runoff at a hillslope scale. Second, it has been already compared with hydrological data, and it was shown to provide both good qualitative and quantitative approximations of the latter (e.g., soil and flow data from the Shale Hills watershed in Pennsylvania, during a controlled experiment on the role of soil moisture and runoff6,8). Third, the input to the model can be easily manipulated allowing us to compare its response to fluctuating rain patterns of different frequencies and intensities. In this paper, we show that same average volume of precipitation, however applied in different patterns (fluctuating rain patterns of different frequencies), can lead to differential responses (significant increases in magnitude) of the runoff coefficient, which is a good predictor of changes in channel-link runoff.10 

The paper is organized as follows: Sec. II introduces the model at the hillslope scale as a system of nonlinear ordinary differential equations (ODEs). The variables, fluxes, and parameters have units and physical meaning and can be directly related to field experiments. The system is nondimensionalized and then two scenarios are investigated: (i) the dependence of the hillslope's response to different constant rainfalls—in Sec. III; and (ii) the dependence of the hillslope's response to fluctuating rain patterns—in Sec. IV. Conclusions are then summarized in Sec. V.

While it is much easier mathematically to work with a nondimensionalized equivalent system, the theoretical results still need to be anchored into hydrological terms. For example, we are interested in working within an admissible range for the rainfall intensity (in millimeters per hour), period of oscillation for the rain pattern (in hours), and volumes per unit hillslope area of water in the saturated and unsaturated zones (in meters). Thus, to ensure that variables and parameters take realistic values, and to allow for straightforward interpretations, we adopt here the following convention: all calculations and mathematical proofs are done in the nondimensional system from Sec. III; then, the scaled variables are transformed back into their physically based counterparts prior to being graphically displayed in figures. The visual representation of the results is always achieved by unit-based plots, and the convention applies to the entire paper.

In order to determine the extent to which oscillatory rainfall pattern could create conditions for large changes in the runoff coefficient that can be associated with flash-floods, we follow a series of analytic steps. First, we imagine the scenario in which the system is forced into a steady state by rainfall of constant intensity. This assumption seems reasonable if the constant input is interpreted loosely as the average value of rainfall occurring randomly over a sufficiently long period of time. The rainfall intensity level and the hillslope specific parameter values determine the type of the equilibrium. It can be shown that given a choice of realistic hillslope parameters, the eigenvalues of the equilibrium are complex within an entire range of precipitation. The equilibrium is attractive and trajectories starting in its vicinity approach it in an oscillatory fashion. The absolute value of those complex numbers is called the system's natural frequency, given that it would be the frequency of the oscillation around the equilibrium in the absence of damping. The positive ratio between the real part and the absolute value is called the damping coefficient and estimates the strength of the decay. Both the natural frequency and the damping coefficient affect the time, it takes a point from the neighborhood to come close enough to equilibrium.

Once it has settled to equilibrium, the system is then subjected to an oscillatory precipitation pattern. The natural frequency and damping coefficient are then used to determine which frequency of periodic precipitation might mostly amplify the soil response, defining the conditions optimal to produce floods. This is called the resonant frequency. To demonstrate the significance of this “feature” of the input related to the hillslope properties, we consider the following experiment: we fix an average precipitation level and calculate its corresponding resonant frequency; we then simulate the system using fluctuating precipitations with same average value as above but with different frequencies (including the resonant frequency); finally, we compute and examine the runoff coefficient. It results that indeed, as predicted by the theory, the amplification of the hillslope's runoff coefficient for same average volume of water is maximal at the resonant frequency. The results are also confirmed by numerical simulations of the nonlinear model.

In this paper, we consider an ODE-based model for the local hillslope-link coupling4 and use it to determined the dynamics of the runoff coefficient when rainfall inputs have different intensities and/or frequencies of occurrence.

The model consists of four ordinary differential equations that account for the interactions between atmospheric inputs and landscape properties, and their implication to the runoff transport dynamics. Four physical control volumes are considered with the following corresponding variables in the model: sp=vwponded/AH, the volume per unit hillslope area of water stored in the ground surface (measured in meters); v=θvunsat/AH, the volume per unit hillslope area of moisture content in the unsaturated zone (measured in meters); a = vsat/AH, the volume per unit hillslope area of the saturated zone (measured in meters); and q = Q/Qr, the (nondimensional) flow going downstream out of the channel at time t̂, where Q is the runoff at the outlet (in m3/s) and Qr is the unit reference discharge (Qr = 1 m3/s). Time t̂ is considered to be measured in minutes. Fluxes between these four control volumes are limited to those listed in Table I with parameters listed in Table II, and the model is reduced to

dspdt̂=c3[1E(t̂)]p(t̂)QplQpu,βdvdt̂=QpuQus,βdadt̂=QusQslQevap,
(1)

and

τdqdt̂=qλ1(q+qin+γQpl+γQsl).
(2)

Here, E(t̂) represents the evapo-transpiration percentage loss function from ponded water; p(t̂) is the precipitation input in [mm/h]; and parameters c1, c2, c3, cevap, γ, and τ are defined by c1=Ksp60hb(m1min1),c2=106αsoilKSAT(2L)60AH(min1),c3=10360(nounit), cevap=Kevap60(min1),γ=106AH60Qr(m1min), and τ=(1λ1)L60vr(Aupstream/Ar)λ2(min). The incoming flux qin is taken here to be zero, given that the Shale Hills watershed is a first-order drainage basin, and the hillslope reference area Ar = 1 km2 is a normalization constant. The assumption is that hb represents an effective soil depth over the total hillslope area AH such that the total volume of the hillslope is Vtot = hbAH. Therefore, variables a and v are restricted to take values only in certain ranges, that is 0 ≤ a ≤ hb and 0 ≤ v ≤ hba. Moreover, given the high density of tree roots in the Shale Hills watershed, the total volume of the hillslope effectively available for water (VT) is only a percentage of Vtot; this observation led to the introduction of a correction factor β such that VT is defined by VT=βVtot=βhbAH with β(0,1]. The values from Table II of all physical parameters of the hydrologic system are the same or very close to those considered in Ref. 4.

TABLE I.

Fluxes in model (1) are defined per unit hillslope area; measured in [m/min].

DefinitionPhysical interpretation
Qpl=c1sp(aares+vvres) Flux from ponded water to the channel link, per unit hillslope area 
Qpu=c1sp(hbav) Flux from ponded water to the unsaturated zone, per unit hillslope area 
Qus=d0(vvres)+d1(vvres)(aares)2+d2(aares)2 Flux exchange between the unsaturated and saturated zones, per unit hillslope area 
Qsl=c2(aares)exp(αNaareshb) Flux for subsurface runoff into the channel link, per unit hillslope area 
Qevap=cevap(aares) Flux for the potential loss of groundwater through either evaporation or plant consumption, per unit hillslope area 
DefinitionPhysical interpretation
Qpl=c1sp(aares+vvres) Flux from ponded water to the channel link, per unit hillslope area 
Qpu=c1sp(hbav) Flux from ponded water to the unsaturated zone, per unit hillslope area 
Qus=d0(vvres)+d1(vvres)(aares)2+d2(aares)2 Flux exchange between the unsaturated and saturated zones, per unit hillslope area 
Qsl=c2(aares)exp(αNaareshb) Flux for subsurface runoff into the channel link, per unit hillslope area 
Qevap=cevap(aares) Flux for the potential loss of groundwater through either evaporation or plant consumption, per unit hillslope area 
TABLE II.

Parameter values used in the numerical simulations of the nonlinear ODE model.

ParameterValueUnitPhysical interpretation
Aupstream 0.07736 km2 Total upstream area 
AH 0.07736 km2 Total hillslope area 
L 420.15 Channel (link) length 
hb 0.56 Effective soil depth to the impermeable layer 
λ1 0.25  Nonlinear exponent for flow velocity function-discharge 
λ2 −0.1  Nonlinear exponent for flow velocity function-upstream area 
vr 0.25 m/s Reference flow velocity 
KSAT 0.01 m/h Saturated hydraulic conductivity 
d0 0.00135 min−1 d0, d1, and d2 are parameters that describe the formof the hillslope-scale recharge relation coupling the statevariables 
d1 0.013 m2min1 
d2 0.008 m1min1 
ares 0.1 ares and vres are the residual storage volumes for thegravity-drained hillslope 
vres 0.3 
β 0.12  Parameter that accounts for the heterogeneities in the soilmatrix porosity and areas inaccessible to water 
αN  Parameter that controls the recession 
αsoil 3.0  Heterogeneity factor for soil saturated hydraulic conductivity 
Ksp 1.75 h1 Parameter that characterizes the overland flow to the channeland the infiltration to the soil 
Kevap 0.025 h1 Recession coefficient for the evaporation process andwater consumption by vegetation, depending on airtemperature, humidity and species of plants on the hillslope 
ParameterValueUnitPhysical interpretation
Aupstream 0.07736 km2 Total upstream area 
AH 0.07736 km2 Total hillslope area 
L 420.15 Channel (link) length 
hb 0.56 Effective soil depth to the impermeable layer 
λ1 0.25  Nonlinear exponent for flow velocity function-discharge 
λ2 −0.1  Nonlinear exponent for flow velocity function-upstream area 
vr 0.25 m/s Reference flow velocity 
KSAT 0.01 m/h Saturated hydraulic conductivity 
d0 0.00135 min−1 d0, d1, and d2 are parameters that describe the formof the hillslope-scale recharge relation coupling the statevariables 
d1 0.013 m2min1 
d2 0.008 m1min1 
ares 0.1 ares and vres are the residual storage volumes for thegravity-drained hillslope 
vres 0.3 
β 0.12  Parameter that accounts for the heterogeneities in the soilmatrix porosity and areas inaccessible to water 
αN  Parameter that controls the recession 
αsoil 3.0  Heterogeneity factor for soil saturated hydraulic conductivity 
Ksp 1.75 h1 Parameter that characterizes the overland flow to the channeland the infiltration to the soil 
Kevap 0.025 h1 Recession coefficient for the evaporation process andwater consumption by vegetation, depending on airtemperature, humidity and species of plants on the hillslope 

Note that E(t̂) represents water loss from the ponded control volume. This loss is mostly due to surface evaporation. As a simplification, we choose to remove this water before it even enters the system by taking only a percentage of precipitation; thus, in the first equation of (1), the precipitation function p(t̂) is multiplied by the factor (1E(t̂)). Contrarily, Qevap represents water loss from the saturated zone of soil. This loss is mostly due to transpiration as plants pull water out of the soil volume through the roots that reach the saturated zone of soil. Here, the flux Qevap cannot take negative values as it is assumed that a is always above ares (see below for more explanations). These two methods of water loss are typically combined into the term “evapotranspiration.”

Note that model (1) has been developed in parallel to a calibration-free, more parsimonious integral-balance ODE system (Ref. 2 and Appendix B in Ref. 4). The latter, however, while showing a good approximation of the channel discharge at the outlet, it fails to reproduce the soil moisture data for the hydrological site under investigation. On the other hand, model (1) was successful in simulating both the flow and soil data. For this reason, we chose it as a “toy-model” for the present study. An explanation of the differences in the definition of fluxes between model (1) and that in Ref. 2, in particular, the way the evapotranspiration terms are chosen, can be found in Ref. 4.

Another peculiarity of the model (1) is the introduction of residual storage volumes ares and vres at zero precipitation. We justify the approach as follows: the definition of the fluxes in the model is phenomenological and they apply during extended periods of time up to centuries. Therefore, when there is no precipitation for such a long time (p(t) = 0), the ODE system will stabilize at a steady state with variables a and v zero. However, the events investigated in this paper occur at much shorter timescale during which the soil is not completely dry even if there is no precipitation for weeks; instead, in this situation, certain nonzero steady state values, say, ares, vres can be rather assumed at p(t) = 0; see Ref. 4 for an additional discussion and comparison with hydrological data.

In Secs. III–V, we will replace the term [1E(t̂)]p(t̂) from (1) with the “effective precipitation” function peff(t̂), assumed to be measured in [mm/h]. The flux between the unsaturated and saturated zones Qus=Qus(a,v) can take a general functional form but it needs to satisfy Qus(ares,vres)=0. However, in this paper, whenever we exemplify our theoretical results, we do that by considering its definition from Table I with coefficients d0, d1, and d2 as in Table II. All other fluxes (per unit hillslope area) involved in the model are kept fixed.

Note that the runoff coefficient RC for system (1) can be easily computed since all water entering the system is either infiltrated into the soil (through Qpu) or flows along the hillslope as surface runoff (through Qpl). Then the runoff coefficient is defined by

RC=QplQpl+Qpu=11f(1a+vhb),

where f is the percentage of available hillslope storage volume calculated at the residual (gravity-drained hillslope) condition

f=1ares+vreshb(0,1).
(3)

In Secs. III–IV, we propose a general approach for a systematic analysis of the dynamics of this hillslope-river local coupling model as a function of the system's intrinsic properties and the external forcing, specifically, the fluctuating patterns of input precipitation peff(t̂). (In this study, we will keep fixed the value of parameter Kevap that controls evaporation/plant water consumption from the soil.) As a first step, we rewrite system (1) in a nondimensional form, and focus mainly on the dynamics of its first three equations and consequently, on the dynamics of the runoff coefficient RC.

The transformation of the input function peff and of the variables sp, v, a, and time t̂

P(t)=103hbKsppeff(t̂),t=t̂Ksp60,S=sphb,V=vvreshb,A=aareshb,
(4)

leads to the nondimensional system

dSdt=P(t)fS,βdVdt=S[fAV]Q̃us(A,V),βdAdt=Q̃us(A,V)αAeαNAαevapA,
(5)

together with the equation for the runoff coefficient

RC=(A+V)/f.
(6)

The parameters α and αevap and the flux Q̃us are defined now by

α=106αsoilKSAT(2L)KspAH,αevap=KevapKsp,Q̃us(A,V)=60hbKspQus(a,v).
(7)

For example, using the Definition Qus(a,v) from Table I, we will work with the following formula for Q̃us(A,V):

Q̃us(A,V)=60d0KspV+60d1hb2KspVA2+60d2hbKspA2=d̃0V+d̃1VA2+d̃2A2.
(8)

Note also that Q̃us(0,0)=0 given that Qus(ares,vres)=0. Moreover, the variables A and V should satisfy a set of mathematical constraints due to the fact that: (i) the water level in the soil cannot exceed the effective soil depth and (ii) there is a residual value for water level in the soil at zero precipitation. This implies that aresaa+vhb. Therefore, A and V are restricted to the following range:

0A1ares/hb,vres/hbVfA.
(9)

Since our goal is to investigate how intrinsic properties of the hillslope-channel link coupling interact with patterns of rain (external signal) to produce flash-flooding conditions, we start by analyzing the steady state condition of the system (5). For that, we will assume initially that the effective precipitation level is constant over time, P(t)=P¯. From physical point of view, this might seem to be an unrealistic condition. However, it is reasonably valid if we interpret P¯ in a more general context; that is, we define P¯ as the averaged value of precipitation during a given interval of time, Tobserve (up to several weeks/from mathematical point of view, this is an asymptotical time). If the rain occurs in a random (non-correlated, non-structured) fashion during the observational time Tobserve, taking the rain average value and using it as constant input to the system is an acceptable assumption.

At constant (average) rain input of P(t)=P¯ with P¯>0, the system (5) reaches an equilibrium (S*,A*,V*) that satisfies the following conditions: S*=P¯/f,0=S*(fA*V*)Q̃us(A*,V*) and 0=Q̃us(A*,V*)αA*eαNA*αevapA*. These allow for an explicit calculation of V* according to the formula

V*=fA*fA*P¯(αeαNA*+αevap),
(10)

and an implicit definition for A* as a function of P¯ by

A*(αeαNA*+αevap)=Q̃us(A*,fA*fA*P¯(αeαNA*+αevap)).
(11)

As already mentioned

S*=P¯/f.
(12)

We exclude here the case P¯=0. In that case, obviously, the only possible equilibrium in system (5) is (0, 0, 0) which corresponds to sp = 0, v = vres, and a = ares.

When computing a solution of Eq. (11), we take into account the range for A* and V* according to (9). This may restrict the range of effective precipitation values P¯, which have physical meaning as the volumes of water a and v need to be positive and under certain bounds. Nevertheless, for all precipitation values within the physical acceptable range, the results presented in the paper do apply. An example showing the dependence of a* and v* (the unit-based counterparts of the equilibria A*, V*), as well as of the runoff coefficient RC* on the constant effective precipitation p¯eff is shown in Figure 1.

FIG. 1.

Coordinates of the equilibrium point (a*,v*) as well as the runoff coefficient RC* as a function of constant, effective precipitation p¯eff. For easier comparison with results in Sec. IV, two runoff coefficient equilibrium values are shown at p¯eff=10 and 20 [mm/h].

FIG. 1.

Coordinates of the equilibrium point (a*,v*) as well as the runoff coefficient RC* as a function of constant, effective precipitation p¯eff. For easier comparison with results in Sec. IV, two runoff coefficient equilibrium values are shown at p¯eff=10 and 20 [mm/h].

Close modal

1. Numerical example

We illustrate these theoretical results by choosing the flux Qus(a,v) as in Table I and parameter values from Table II. Therefore, its corresponding Q̃us(A,V) is defined by (8). The curves (P¯,A*) and (P¯,V*) show the coordinates of the equilibrium point for any given constant precipitation level and are computed according to (9)(11). Their unit-based counterparts (p¯eff,a) and (p¯eff,v) are then plotted in Figure 1 (upper panels). Although it seems counterintuitive that v, the amount of water in the unsaturated zone, should decrease as more water is added, this can be explained by the increase in the amount of water in the saturated zone. That is to say, as more water enters the system, the soil becomes saturated, so the unsaturated zone, and the amount of water in it, becomes smaller. The graph of the runoff coefficient RC at the equilibrium offers a better measure of the water dynamics on the hillslope and it is included in Figure 1 (lower panel).

The eigenvalues of the equilibrium point (S*,V*,A*), their type (real versus complex), and sign of their real part (positive versus negative) are important indicators of the local dynamics of (5). This is because the nonlinear system is topologically conjugated to its linearization, in a sufficiently small neighborhood of any hyperbolic equilibrium (Ref. 12; the Hartman-Grobman theorem). Note that an equilibrium is said to be hyperbolic if all its eigenvalues have non-zero real part. In other words, nonlinear systems “look alike” their linearization near hyperbolic equilibria. For example, if the equilibrium has a negative real eigenvalue λ1 then, in its neighborhood, all trajectories collapse to the equilibrium in an exponentially decaying fashion along the (eigen)direction associated with λ1. On the other hand, if the equilibrium has a pair of complex conjugated eigenvalues λ2,3 with negative real part (λ2=λ;λ3=λ¯;Re(λ)<0) then, in its neighborhood, the trajectories tend to the equilibrium in a decaying oscillatory (spiraling down) fashion along the two-dimensional manifold associated with λ2,3.

The eigenvalues of the equilibrium point (S*,V*,A*) are the eigenvalues of the Jacobian matrix evaluated at (S*,V*,A*)

(f00fA*V*βa11a120a21a22),

where

a11=1βQ̃usV(A*,V*)P¯fβ,a12=1βQ̃usA(A*,V*)P¯fβ,a21=1βQ̃usV(A*,V*),a22=1βQ̃usA(A*,V*)1β[αeαNA*(1+αNA*)+αevap].
(13)

From here on, we will use the notation J for the 2-by-2 sub-block matrix J=(aij).

Therefore, the eigenvalues of the equilibrium point are λ1=f and λ2,3, the roots of the quadratic equation λ2tr(J)λ+det(J)=0, where

tr(J)=a11+a22;det(J)=a11a22a12a21,
(14)

are the trace and the determinant of the sub-block J of the Jacobian matrix. Note that both eigenvalues λ2,3 have negative real part if and only if tr(J) < 0 and det(J) > 0. Moreover, they are complex conjugated if and only if tr(J)24det(J)<0.

Under the circumstances described above (λ1<0 and λ2,3C\R,Re(λ2,3)<0), the system tends to an attractive equilibrium with trajectories approaching it in a spiral fashion when projected on the (V, A)-plane, during the observational time Tobserve. For the dynamics of the system, two quantities are important: the absolute value of the complex eigenvalues and the positive ratio between their real part and the absolute value. Say λ2=λ and λ3=λ¯; then

ω0=|λ|=Re(λ)2+Im(λ)2,ζ=|Re(λ)||λ|=Re(λ)|λ|=Re(λ)Re(λ)2+Im(λ)2.
(15)

They are called the system's natural frequency (ω0) and the damping coefficient (ζ). The natural frequency is the frequency at which a trajectory of initial condition close enough to the equilibrium would oscillate around it if no damping were present. On the other hand, the damping coefficient is a reliable measure for how fast the spiralling down trajectory approaches the equilibrium under the damping condition (i.e., when the complex eigenvalues have negative real part). In other words, ζ represents the degree to which the response of the system is diminished. Obviously, both ω0 and ζ depend upon the parameters of the system. A large damping ratio (ζ close to 1) means that the system will quickly reduce the magnitude of its response around the equilibrium point. We say, that the system works to “dampen” the effects of the applied forcing. To the contrary, a small damping ratio (typically ζ<1/20.7) means that the response of the system maintains a greater portion of the magnitude of the applied forcing for a longer period of time.

We say, that under the above circumstances (complex eigenvalues λ, λ¯ with negative real part) and ζ satisfying given inequalities to be determined (see the inequalities (16) and formulas (17)), the system is primed to the pre-flooding condition. We introduce this terminology according to the following logic: the pre-flooding condition described above is an equivalent to the resonance condition from classical models of forced linear systems.14 As a reminder, the resonance condition is defined as a situation corresponding to a critical frequency value at which the amplitude of the system's response signal is significantly amplified in comparison to those obtained at neighboring frequency values. Therefore, by similarity, one would expect that a fluctuating rain pattern of a certain frequency ω might lead to a maximal hillslope-channel link response or, equivalently, to larger than usual peak amplitudes in the runoff coefficient, as seen in flash-floods. The question is: “How do the intrinsic properties of the hillslope contribute to such phenomenon?” We show in Sec. IV how this is possible, and exemplify the answer with results obtained for system (1).

Given that we are interested to study the dynamics of system (5) in what we called as being “the pre-flooding” condition, we focus now only on those values P¯ for precipitation that lead to complex eigenvalues λ2,3=λ,λ¯ with negative real part. (A condition that determines a choice interval for rain frequency ω will be then determined in Sec. IV A.) Therefore, for a given parameter set that characterizes the hillslope-channel link coupling, we investigate the flow response only when the average rain level P¯ falls in the range defined by

{tr(J)<0,det(J)14tr(J)2>0.
(16)

Moreover, since tr(J)=2Re(λ) and det(J)=|λ|2, the system's natural frequency ω0 and the damping ratio ζ defined by (15) can be written in terms of the trace and determinant of the sub-block matrix J as

ω0=det(J)andζ=tr(J)2det(J),
(17)

and then λ=ω0ζ+iω01ζ2. Figure 2 illustrates the range of averaged effective precipitation (peff > 0.4 [mm/h]) for which condition (16) is valid, and the graphs of the damping ratio ζ and the natural frequency ω̂0=ω0Ksp in [h1] and period T̂0=2π/ω̂0 in [h] versus p¯eff. (Note that for peff < 0.4 [mm/h], the eigenvalues of the equilibrium are real and negative, and all trajectories starting in the vicinity of the equilibrium decay to it exponentially fast. So neither ω0 nor ζ make sense in this interval.)

FIG. 2.

The damping ratio ζ(top), the natural frequency ω̂0 in units [h1](center), and the period associated with the natural frequency (bottom) for the range of averaged effective precipitation peff where the equilibrium is an attractive focus.

FIG. 2.

The damping ratio ζ(top), the natural frequency ω̂0 in units [h1](center), and the period associated with the natural frequency (bottom) for the range of averaged effective precipitation peff where the equilibrium is an attractive focus.

Close modal

We define the time to the asymptotic solution with respect to the ponded water sp = hbS as the minimum time necessary for sp to reach a level very close to its steady solution. The “very close level” is assumed here to be under 1% relative error. Given that the first equation in (5) is linear and decoupled, its solution can be determined explicitly, and it is

S(t)=S*+(S0S*)eftS*ast,

where S0 is any initial condition and S*=P¯/f is the equilibrium value at P(t)=P¯. Therefore, the asymptotic time to the steady state sp* can be computed according to the inequality

|spsp*sp*|=|SS*S*|=|S0S*|S*eft0.01
teq,S=max(0;1fln|fS0P¯0.01P¯|).
(18)

We will assume the following scenario: (i) the hillslope-link coupled system is subject to rain input of average magnitude P¯ but without any particular structure (random variability about P¯) for a time period at least equal to the time to equilibrium S*. This is computed according to (18) for some initial condition S0. Then, (ii) input P¯ is maintained in order to allow for variables V and A to reach their equilibrium (V*,A*) as well, with exponential decay eω0ζt; we define this time-period by

teq,AV=max(0;1ω0ζln|A0A*0.01A*|;1ω0ζln|V0V*0.01V*|).

Consequently, the time to equilibrium teq is estimated by teq=max(teq,S;teq,AV)

teq=max(0;1fln|fS0P¯0.01P¯|;1ω0ζln|A0A*0.01A*|;1ω0ζln|V0V*0.01V*|),

or, written in units (hours), for initial conditions (sp,0,v0,a0) measured in meters

t̂eq=1Kspmax(0;1fln|103Kspfsp,0p¯eff0.01p¯eff|;
1ω0ζln|a0a*0.01(a*ares)|;1ω0ζln|v0v*0.01(v*vres)|)[h].
(19)

An example that estimates the time to equilibrium t̂eq when p¯eff varies is shown in Figure 3.

FIG. 3.

The quantities 1/(ζω0) and 1/f (unitless; 1/f = 3.5) (top). The time for the system to reach within 1%, 5%, and 10% of the equilibrium conditions under constant precipitation as determined in Eq. (19) (bottom).

FIG. 3.

The quantities 1/(ζω0) and 1/f (unitless; 1/f = 3.5) (top). The time for the system to reach within 1%, 5%, and 10% of the equilibrium conditions under constant precipitation as determined in Eq. (19) (bottom).

Close modal

We calculate, similar to (19), and plot the time t̂eq (in hours) required for the system to come within 1%, 5%, and 10% of its equilibrium value. For example, at p¯eff=10 mm/h, the system takes roughly 10 h to come to equilibrium, within 10%. Better approximations of the equilibrium value can be obtained if waiting approximately 14 h (for an error within 5%) or, say, 22 h (for an error within 1%). Because the difference between the times associated with 1% and 5% is not significantly larger than the difference between times associated with 5% and 10%, we can say, that the system reached equilibrium in a steady way. (This is opposed to a possible, alternative situation in which the system would have approached the equilibrium very quickly, within 5%, then slowed down and taken a very long time to get under 1% difference to it.)

The time to equilibrium is found using the values s0 = 0, a0 = ares, and v0 = vres in Eq. (19). This formula takes the maximum of four expressions. Because of these initial conditions, the natural logs are constant for any fixed percentage 1%, 5%, or 10%. Then the values of 1ω0ζ and 1f determine which term will determine the maximum. For smaller values of average precipitation, the third and fourth expressions are larger because 1/ω0ζ>1/f (see top panel in Figure 3). For average precipitation values larger than 17 mm/h, however, the second term dominates in Eq. (19), and the time to equilibrium becomes constant for all precipitation values (lower panel in Figure 3).

We investigate now the consequences of applying a fluctuating rain pattern to the hillslope-link coupled system (1) primed to the “pre-flooding” condition as defined in Sec. III. We consider the term [1E(t̂)]p(t̂) replaced by the effective precipitation function peff(t̂). The fluctuations about the average value P¯ are assumed to have amplitude M with 0MP¯ in order to always have peff(t̂)0, so

peff(t̂)=p¯eff+meffsin(ωKspt̂60)=hbKsp103(P¯+Msin(ωt))[mm/h].

As a reminder from (4), time t̂ is measured in minutes, Ksp is in h1 and t=t̂Ksp60 is the nondimensional time. Then, obviously, the nondimensional input to system (5) is

P(t)=P¯+Msin(ωt)with|P(t)P¯P¯|1.

Obviously, while the effective precipitation pattern is fluctuating, on average, its value is still P¯. We are interested to determine if there exits any particular choice of amplitude M=M(P¯) and frequency ω=ω(P¯) that leads to larger responses in the runoff coefficient RC (i.e., considerably larger values reached during oscillation) than those anticipated by the equilibrium value. For example, say, M=P¯ and the rainfall is oscillatory. Then, during the oscillation, the system receives precipitation input with values between 0 and 2P¯. The steady state analysis predicts certain values for the runoff coefficients for each of P[0,2P¯] (see Figure 1, lower panel). The questions are: Is it possible to exceed the RC*(P) values if controlling P to change in an organized way (regular oscillation) between 0 and 2P¯? If yes, how big are those differences, and at which frequency they occur?

In order to address such questions, we start by investigating numerically the dynamics of the runoff coefficient in the nonlinear system (1). Note that at peff = 10 mm/h, the steady state RC* is RC*(10)=0.497 while at peff = 0 and peff = 20 mm/h we find RC*(0)=0 and RC*(20)=0.501, respectively (Figure 1). Moreover, on the entire interval spanning [0, 20] mm/h, the runoff coefficients RC* can reach only values up to approximately 0.51. However, if peff varies between 0 and 20 mm/h in an organized (oscillatory) fashion, with certain frequency ω, RC can reach significantly larger values. Thus, RCpeak may take values up to 0.8; see Figure 4. By comparison, in the steady state condition, RC does not come close to 0.8, not even for very high constant precipitation rates such as 35 mm/h; Figure 1. Those observations justify a more thorough analysis of RC under fluctuating rain patterns.

FIG. 4.

Simulated RC values using p¯eff=10 [mm/h] and ω̂=0.1 [1/h] (top left), ω̂=0.275 [1/h] (top right), ω̂=0.424 [1/h] (this is near the resonant frequency, ω̂R=0.4242 [1/h]) (center left), ω̂=0.625 [1/h] (center right), and ω̂=0.8 [1/h] (bottom). The corresponding periods (T̂=2πω̂) in these examples are 62.8 [h] (top left), 22.8 [h] (top right), 14.8 [h] (center left), 10.1 [h] (center right), and 7.9 [h] (bottom), respectively.

FIG. 4.

Simulated RC values using p¯eff=10 [mm/h] and ω̂=0.1 [1/h] (top left), ω̂=0.275 [1/h] (top right), ω̂=0.424 [1/h] (this is near the resonant frequency, ω̂R=0.4242 [1/h]) (center left), ω̂=0.625 [1/h] (center right), and ω̂=0.8 [1/h] (bottom). The corresponding periods (T̂=2πω̂) in these examples are 62.8 [h] (top left), 22.8 [h] (top right), 14.8 [h] (center left), 10.1 [h] (center right), and 7.9 [h] (bottom), respectively.

Close modal

Numerical simulations of the nonlinear model (1) with parameters from Table II indicate a nonlinear dependence of the RCpeak on the frequency of oscillation ω (Figure 4). In this example, the average rainfall is peff = 10 mm/h. The peak of RC increases with ω until it reaches a maximum at a certain ωR; then for larger frequencies ω > ωR, the peak decreases with ω.

To investigate the nonlinear relationship RCpeak=RCpeak(ω), we employ results from the linearization theory. Note that the first equation in (5) being linear and decoupled, it allows us to explicitly compute its solution

S(t)=S*+Sf(t)+(S0S*+ωMf2+ω2)eft[S*+Sf(t)]ast,
(20)

where S0 is any initial condition, S*=P¯/f is the equilibrium value at P(t)=P¯, and Sf(t) is the fluctuating (and non-damping) term

Sf(t)=Mf2+ω2sin(ωt+ψ),ψ=arctan(ω/f).
(21)

The function S*+Sf(t) is called the asymptotic solution for S.

By assuming that a rain input of average magnitude P¯ but without any particular structure was applied in (5) for a time period at least equal to teq, the system's trajectory is in 1% error from the system's equilibrium. Then we assume that the rain pattern gathers structure, with fluctuations about P¯ of certain amplitude M and frequency ω. This should occur for longer time than tfluct,S in order to stabilizes the solution S to its fluctuating non-damping part S*+Sf(t). Time tfluct,S is computed now by starting with initial condition S0=(P¯+0.01P¯)/f and it is tfluct,S=max(0;1fln|1+fωf2+ω2M0.01P¯|). Consequently, the total transient time to the fluctuating dynamics that we aim to investigate will be ttransient=teq+tfluct,S or, written in units (hours), for initial conditions (sp,0,v0,a0) measured in meters

t̂transient=t̂eq+1Kspmax(0;1fln|1+fωf2+ω2meff0.01p¯eff|)[h],
(22)

with t̂eq given by (19).

Under such circumstances, we can use formulas (20) and (21) to reduce (5) to a two-dimensional nonlinear system of differential equations in V and A. After time ttransient has passed (Tobserve > ttransient), the input to equation V is S(t)P¯/f+Sf(t) and system (5) reduces to

βdVdt=(P¯/f+Sf(t))(fAV)Q̃us(A,V),
βdAdt=Q̃us(A,V)αAeαNAαevapA.

Moreover, if the point (V, A) along the trajectory is in the neighborhood of (V*,A*) at t>ttransient, then it is reasonable to consider the linearization of the two-dimensional system about (V*,A*) and to approximate the dynamics of the non-linear system with that of the linearized counter-part. We take this approach in Sec. IV A.

We start our analysis by finding the expansion of the nonlinear (V, A)-system about its equilibrium (V*,A*) at P¯. Through a linear change of variables, say, V=V*+y1 and A=A*+y2, the system becomes

βdy1dt=(fA*V*)Sf(t)(P¯f+Sf(t)+Q̃usV(A*,V*))y1(P¯f+Sf(t)+Q̃usA(A*,V*))y2+O(y12,y22,y1y2),βdy2dt=Q̃usV(A*,V*)y1+(Q̃usA(A*,V*)αeαNA*(1+αNA*)αevap)y2+O(y12,y22,y1y2),

with linearization

dy1dt=(fA*V*)1βSf(t)+(a111βSf(t))y1+(a121βSf(t))y2,dy2dt=a21y1+a22y2,
(23)

where coefficients aij are defined by (13). Therefore, in first order approximation, we have A(t)A*+y2(t) and V(t)V*+y1(t). The runoff coefficient RC=(A+V)/f is then approximated by RC(t)RC*+(y1(t)+y2(t))/f.

The construction of an approximate solution for system (23) is provided in the Appendix. As a consequence, we determine a formula for the average value of the runoff coefficient, defined by RCavg=1T0TRC(t)dt with T=2πω as a function of the frequency ω of the rain pattern P(t) (see Subsection 3 of the Appendix for more details). This is

RCavg(ω)=RC*+Ω(ω)1+Ω(ω),
(24)

with

Ω(ω)=M2P¯2(12ω0ζfβP¯)[ω2ω02(12ω0ζfβP¯)]2f2β4ω02(f2+ω2)[(ω02ω2)2+4ω2ω02ζ2].
(25)

Moreover, for an oscillatory input pattern with frequency ω, the runoff coefficient RC can be approximated (by excluding the higher order harmonics) by

RCRCavg(ω)+A(ω)sin(ωt+ψ+θ),
(26)

with ψ and θ defined by (21) and (A5), and

A(ω)=M(1RC*)β[1+Ω(ω)]ω2+(P¯fβ2ω0ζ)2(f2+ω2)[(ω02ω2)2+4ω2ω02ζ2].
(27)

Therefore, when an oscillatory rain pattern with average value P¯ and frequency ω is applied, the runoff coefficient reaches its peak value RCpeak(ω)RCavg(ω)+A(ω) at certain moment during the cycle.

Interestingly, the hillslope response (defined here as the peak of the runoff coefficient, RCpeak) to the fluctuating rain pattern P(t)=P¯+Msin(ωt) with same average value P¯, depends nonlinearly on the choice of frequency ω. There exists a particular value ωR where the peak response is maximal. In fact, as shown in Figure 5, at ωR the peak response RCmax=RCpeak(ωR) shows an increase up to 20% by comparison with other values of potential frequencies in a range close to it, when small amplitude variations occur in precipitation intensity (in a range of 25% of average value; Figure 5 (bottom)). When larger variations in amplitude are considered (meff=0.5p¯eff or meff=p¯eff; Figure 5 (center) and (top)), RCmax reaches at ωR even more significant levels. These show approximately 40% and 60% increases compared to nearby values of ω or, to the steady state value RC* at p¯eff (RC*=0.497 at p¯eff=10mm/h in Figure 5). For this reason, we use the term resonant frequency when referring to this particular choice of the frequency of the oscillatory input, ωR and argue that flooding has higher chance to occur under those circumstances. The existence of ωR and presence of a fluctuating precipitation pattern of frequency ω close to ωR is what we call a pre-flooding condition.

FIG. 5.

Comparing RCpeak from the linearized system (solid line) and the original nonlinear system (dashed line) for meff=p¯eff(top), meff=0.5p¯eff(center), and meff=0.25p¯eff(bottom). Here, p¯eff=10 [mm/h] at which the corresponding steady state runoff coefficient is RC*=0.497. The maximal percentage increases of the runoff coefficient from simulations of the original nonlinear system and from the linearized system are labeled on the plots. The right panels depict the rain patterns being applied to both situations, with ω̂ in the interval [0.1, 0.8] 1/h.

FIG. 5.

Comparing RCpeak from the linearized system (solid line) and the original nonlinear system (dashed line) for meff=p¯eff(top), meff=0.5p¯eff(center), and meff=0.25p¯eff(bottom). Here, p¯eff=10 [mm/h] at which the corresponding steady state runoff coefficient is RC*=0.497. The maximal percentage increases of the runoff coefficient from simulations of the original nonlinear system and from the linearized system are labeled on the plots. The right panels depict the rain patterns being applied to both situations, with ω̂ in the interval [0.1, 0.8] 1/h.

Close modal

Note that the error between the numerical results of the nonlinear system and the numerical results of the linearized, theoretical, approximation increases when meff/p¯eff=M/P¯=ϵ approaches 1. This can be explained by the fact that mathematically, linear approximations to nonlinear system yield best results when ϵ0 given that S(t)S*+Sf(t)=S*(1+O(ϵ)). Another source of error comes from neglecting the higher order harmonics in the approximating solution (A4) of the linearized system. We were not able to show that those terms are sufficiently small relative to the amplitude RCavg(ω)+A(ω), to be safely omitted in RC. In fact, numerical comparisons do indicate an increased error close to the resonant frequency ωR between the peak RCpeak of the runoff coefficient, and its corresponding value computed from the nonlinear system (see Figure 5). On the contrary, further away from ωR and for small values of frequencies ω, there is no significant error between those peaks. Interestingly, in this case, the higher order harmonics seem to rather influence the shape of the RC signal than the peak itself (see, for example, first panel in Figure 4—here the time series of RC computed from the nonlinear system is less sinusoidal than predicted by formula (26)). Nevertheless, as illustrated by Figure 5, even for ϵ large (ϵ = 1) and by dropping the higher harmonics in the approximating solution of the linearized system, the latter still provides a good qualitative (if not quantitative) approximation of the runoff coefficient's peak dynamics.

Given ωR, the corresponding period is TR=2πωR. The unit-based counter-parts of the resonant frequency and period are then determined by the equations ω̂R=KspωR in [h1] and T̂R=TRKsp in [h], and they are plotted in Figure 6. Note that at small constant precipitation rates (≤1.6 [mm/h]), there is no local maximum in RC for either the linearized system or the simulated nonlinear model over different values of ω.

FIG. 6.

Theoretical resonant frequency, ω̂R, and the corresponding period, T̂R=2πω̂R, for different average rainfall intensity.

FIG. 6.

Theoretical resonant frequency, ω̂R, and the corresponding period, T̂R=2πω̂R, for different average rainfall intensity.

Close modal

To illustrate the dynamics of the soil variables in the nonlinear model (1) when a fluctuating rain pattern is applied, we draw the phase-plane (v, a) for p¯eff=10[mm/h] and its resonant frequency ω̂R=0.4242[1/h]. The trajectory oscillates around the equilibrium point (v*,a*)=(0.2,0.3) at p¯eff=10[mm/h] which corresponds to RC*=0.497. For comparison, we also plotted the projection on (v, a) of the curve of steady state points for other constant precipitation intensities p¯eff in the range [0 mm/h, 35 mm/h]. These are represented in Figure 7. As observed, by introducing structure to the rainfall, similar (higher) levels of groundwater (a) as seen for higher rain intensity can be obtained with reduced rain intensity as well.

FIG. 7.

The dotted line represents the steady state points for the given precipitation rates in the range [0 mm/h, 35 mm/h]. The solid line is the trajectory for p¯eff=10 [mm/h] and ω̂R using meff=p¯eff(top), meff=0.5p¯eff(center), and meff=0.25p¯eff(bottom). Here, the equilibrium values at p¯effmeff,p¯eff, and p¯eff+meff are represented by the symbols △, ○, and ◇, respectively.

FIG. 7.

The dotted line represents the steady state points for the given precipitation rates in the range [0 mm/h, 35 mm/h]. The solid line is the trajectory for p¯eff=10 [mm/h] and ω̂R using meff=p¯eff(top), meff=0.5p¯eff(center), and meff=0.25p¯eff(bottom). Here, the equilibrium values at p¯effmeff,p¯eff, and p¯eff+meff are represented by the symbols △, ○, and ◇, respectively.

Close modal

In this paper, we analyzed the hydrological ordinary differential equation model presented in Ref. 4, and its ability to predict conditions that could lead to flash-flooding while undergoing an oscillating precipitation pattern. The measure to determine if pre-flooding conditions occur is the runoff coefficient, which is directly related to the amount of antecedent soil moisture.

To isolate the effects of rainfall fluctuations on the hillslope response, we propose the following “experiment”: first, to bring the system to equilibrium by applying constant rainfall; second, to perturb its state by an oscillatory input that in average introduces the same volume of water into the system; third, to monitor the dynamics of the runoff coefficient, in particular, the maximum value it can reach.

Note that, although constant rainfall may be a nonrealistic scenario, it is acceptable if rainfall intensity is interpreted on average over a sufficiently long period of time. After reaching a steady state under constant rainfall, the system is linearized about the equilibrium, and eigenvalues of the Jacobian are determined. Since these eigenvalues are complex for all precipitation values considered here (in a realistic range), the equilibrium is attractive with trajectories that spiral toward it. The spiraling nature of the equilibrium introduces a natural frequency and a damping ratio inherent to the system that describe how quickly the trajectories reach the steady state.

We then compare theoretical results of the linearized system with numerical simulations of the nonlinear system when the rainfall is assumed to oscillate around the same average value as before but having certain frequencies. We found in both cases that there exists a particular frequency of precipitation that maximizes the runoff coefficient response. We call it the resonant frequency. This finding shows that the pattern of the rainfall, even though, say, same average precipitation volume is applied, has significant consequences for the setting of the hillslope state to conditions optimal for flash-flooding.

An interesting observation stemming from our analysis is that man can relatively easily intervene at the hillslope scale to eliminate the occurrence of the “pre-flooding conditions.” According to our definition of pre-flooding conditions, those are possible only if the system's equilibrium is an attractive focus. However, the type of equilibrium (attractive focus versus attractive node) depends strongly on the system's parameters, some of them allowing for human intervention. For example, parameters such as αN, αsoil, and Ksp, or coefficients d0, d1, d2 from the definition of the flux Qus (see Table II) play a significant role in setting the type of the equilibria. Numerical investigation of model (1) shows indeed that changes in values of the above mentioned parameters lead to equilibria with only real eigenvalues (not shown here). In the corresponding neighborhood, the amplitude of the runoff coefficient cannot be amplified in this case simply because of the changes in the rainfall frequency (as seen for pre-flooding conditions). An increase in RC is rather driven by the magnitude of the rain storm. As mentioned, by changing parameters such as αN, αsoil, Ksp, d0, d1, or d2, the equilibrium point can be made a node for all precipitation values in the realistic range and so pre-flooding conditions are avoided. Conversely, reverse manipulations of these parameters may prime the hillslope for pre-flooding. In particular, Ksp and αsoil depend on the soil properties. Modifying the soil permeability will affect water transport at the hillslope scale and so its response to rain events. Dramatic changes in local vegetation or intensive agriculture activity are classical examples when those parameters are affected. On the other hand, αN and d0, d1, d2 are determined by the geometry of the hillslope, and they can also be influenced by human exploitation of land. While mathematically, this change in the system's predicted dynamics is not a bifurcation per se (since the equilibrium is still hyperbolic and attractive); in practice, it can be interpreted as such and is of general hydrological interest.

As a final note, we point out that model (1) undergoes Hopf bifurcations in certain parameter ranges. However, this mathematically interesting case is quite difficult to interpret from hydrological point of view and it was postponed for future study.

The authors would like to thank their collaborators Witold Krajewski, Ricardo Mantilla, and Scott Small from the Iowa Flood Center for suggestions and helpful discussions that led to this project. This work was supported by the National Science Foundation Award No. DMS-1025483.

1. Reduction of the linearized system to its Jordan form

System (23) defined by (13) and periodic function Sf(t) defined by (21) is transformed to its Jordan normal form. In addition, we make a time-change that brings the periodic coefficients of the linear system to period T̃=1. The time-change is defined by ωt+ψ=:2πt̃. Note also that in complex plane (here i=1), sin2πt̃=12i(e2πit̃e2πit̃).

The stationary matrix (alj) has distinct complex conjugated eigenvalues λ1 and λ2 which are the roots of the characteristic equation λ2+2ζω0λ+ω02=0 with ζ,ω0 defined by (17), say, λ1=ζω0+iω01ζ2,λ2=λ¯1. We consider the matrix Z whose columns are the corresponding eigenvectors, and its inverse Z1

Z=(a12λ2a22λ1a11a21),Z1=1det(Z)(a21a22λ2a11λ1a12),det(Z)=a12a21(λ1a11)(λ2a22)=2(λ1a11)iω01ζ2,

and make the change of variables X=(x1,x2)T=Z1(y1,y2)T=Z1Y. (Here, T stands for the transpose.) This leads to the system

Ẋ=(2πλ1ω002πλ2ω)X+(abcd)Xe2πit̃(abcd)Xe2πit̃+(eg)e2πit̃(eg)e2πit̃=:B0X+B1Xe2πit̃B1Xe2πit̃+F0e2πit̃F0e2πit̃,
(A1)

with a, b, c, d defined by Λ1=Mπβωf2+ω2ia21det(Z),Λ2=Mπβωf2+ω2i(a11λ1)det(Z),a=Λ1Λ2c,b=Λ1Λ2d,c=(λ1a11+a12)Λ2,d=(λ1a11a21)Λ2,e=(fA*V*)Λ1 and g=(fA*V*)Λ2.

The order of magnitude of the above coefficients can be estimated with respect to parameter β, assumed to be small (e.g., β=0.12 as considered in Sec. II). Equations (13), (14), and (17) imply that βalj=O(1) and βω0=O(1) while ζ=O(1); therefore βλ1,2=ζβω0±iβω01ζ2=O(1) and β2det(Z)=O(1).

Then we have Λ1,Λ2=2πωO(1) and a,b,c,d=2πωO(β1), and e,g=2πωO(1). Consequently, the matrices B0, B1, and F0 have different order of magnitude with respect to β as they are defined by B0=2πβ1ωO(1),B1=2πβ1ωO(1),F0=2πβ1ωO(β) as β1.

2. Description of the transient term of the solution

The solution of Eq. (A1) in the homogeneous case (F0=0) takes the form Xhomog(t̃)=C1eμ1t̃Qμ1(t̃)+C2eμ2t̃Qμ2(t̃), where C1, C2 are constants that depend on the initial condition and eμ1t̃Qμ1(t̃),eμ2t̃Qμ2(t̃) are two independent (nontrivial) solutions of the form

eμt̃Q(t̃)=eμt̃k=Qke2πikt̃.

Indices k are integers. With notation I2 and O2 for the two-dimensional unity and null matrices, Qk are two-dimensional constant vectors from the complex plane11 satisfying

B1Qk1+[B0(μ+2πik)I2]QkB1Qk+1=O2,kinteger.

Then the solution of (A1) is the superposition of two functions: Xhomog and a particular solution of the nonhomogenous equation. We assume here that the oscillatory rainfall covers a long enough time interval to allow for the solution to stabilize. So we are interested only in the asymptotic behavior of the runoff coefficient. It is therefore important to verify that the exponentials eμ1,2t̃ in the definition of Xhomog shrink with time such that this part of the general solution is transient. This is equivalent to proving that Re(μ1,2)<0. However, unfortunately, Re(μ1,2) depend on the peculiarities (equations and parameters) of the system and there is no general theorem to predict their sign. So we do need to estimate them.

According to the theory of approximating solutions for linear systems of differential equations with periodic coefficients,11z1=eμ1 and z2=eμ2 are the two roots of the equation

(zr1)(zr2)+K12(z+r1)(zr2)+K22(zr1)(z+r2)=0,

with notation r1=e2πλ1ω and r2=e2πλ2ω and K1, K2 defined by some complicated infinite determinants that are of the same order as a, b (for K1), respectively, c, d (for K2). Given that a,b,c,d=2πβ1ωO(1), we have K1,K2=2πβ1ωO(1) as well.

Note that the quadratic equation in z is equivalent to (1+K1+K22)z2(r1+r2)(1(K1K2)(r1r2)2(r1+r2))z+r1r2(1K1+K22)=0 and we have r1,r2=eO(β1) and K1,K2=O(β1) so r1,r2K1,K2 as β1. Then the roots z1, z2 satisfy z1r1(1K1) and z2r2(1K2) so μ1=2πλ1ω+2πωRe(λ1,2)O(β) and μ2=2πλ2ω+2πωRe(λ1,2)O(β) with Re(λ1,2)=O(β1).

As a consequence, the asymptotic behavior of Xhomog(t̃) is prescribed by the limit of eRe(μ1)t̃ and eRe(μ1)t̃. In original time t (ωt+ψ=:2πt̃) and coordinates Y, the dynamics of the transient term of the solution is given by the limit of eω2πRe(μ1)t=e(1+O(β))Re(λ1)t and eω2πRe(μ2)t=e(1+O(β))Re(λ2)t, where β is small. Both exponentials tend to zero given that the dominant exponent is Re(λ1,2)=ζω0=O(β1)<0.

3. A particular solution for the nonhomogenous equation

We are interested to determine now a particular solution for the nonhomogenous equation (A1), say, Xp(t̃). If this can be found then the general solution of (A1) is X(t̃)=Xhomog(t̃)+Xp(t̃)=C1eμ1t̃Qμ1(t̃)+C2eμ2t̃Qμ2(t̃)+Xp(t̃) and so X(t̃)Xp(t̃) as t̃. The steady, non-decaying, oscillatory function Xp(t̃) is a good approximation of the solution.

Given the form of Eq. (A1), we are searching for a function Xp(t̃) such that Xp(t̃)=k=Rke2πikt̃ with index k spanning the set of integers. Introducing Xp in (A1), it results that the two-dimensional vectors Rk should satisfy

B0Rk2πikRk+B1Rk1B1Rk+1=0,k1,1,
B0R12πiR1+B1R0B1R2=F0,
B0R1+2πiR1+B1R2B1R0=F0.

By considering the complex inner product between two dimensional vectors U,V=(u1,u2)T·(v1,v2)T=u1v¯1+u2v¯2, we determine first the solution of the adjoint equation B1*F=0 and choose F=(λ¯1a11,a21)T accordingly. Then, since B1(Rk+1Rk1),F=Rk+1Rk1,B1*F=0 and F0,F=0, we note that all vectors Rk satisfy the condition (B02πkiI)Rk,F=0. Moreover, since we are looking for a real solution Yp(t)=ZXp(t)=k=ZRke2πikt̃=ZR0+k=1(ZRke2πikt̃+ZRke2πikt̃) with matrix Z defined in Subsection 1 of the Appendix, we need ZR0 to be real and ZRk=ZRk¯ for k1. That determines vectors Rk=(v1,k,v2,k)T and Rk=(v1,k,v2,k)T with coordinates satisfying v1,k=1detZa21(λ2ikω)rk,v2,k=1detZ(a11λ1)(λ1ikω)rk and v1,k=v1,k(λ2+ikω)r¯k(λ2ikω)rk,v2,k=v2,k(λ1+ikω)r¯k(λ1ikω)rk, for some complex constant rk, and

ZRk=rk(a22ikωa21),ZRk=r¯k(a22+ikωa21).

However, ZB0Z1=2πω(a11a12a21a22),ZB1Z1=iπMωβf2+ω2(1100),ZF0=iπM(fA*V*)ωβf2+ω2(10). Since (ZB0Z1)(ZR0)=ZB0R0=ZB1(R1R1)=(ZB1Z1)(ZR1ZR1) and with notation r1=L+iH, we obtain

ZR0=M[Lω+H(a21a22)]det(J)βf2+ω2(a22a21).

Let us observe that we can determine R1 and R1 by choosing R2=R2=0. This is true because (ZB0Z12πiI)ZR1=Z(B02πiI)R1=(ZB1Z1)ZR0ZF0 and similar (ZB0Z12πiI)ZR1¯=(ZB0Z1+2πiI)ZR1=Z(B0+2πiI)R1=(ZB1Z1)ZR0+ZF0=(ZB1Z1)ZR0¯ZF0¯. These allow us to solve for L and H from the equation

L(ω02ω2)H2ωω0ζ+i[H(ω02ω2)+L2ωω0ζ)]=iM(fA*V*)2βf2+ω2iM(a22a21)2βf2+ω2M[Lω+H(a21a22)]det(J)βf2+ω2,

which implies

H=ω02ω22ωω0ζL,

and

L=(a11+a22)(fA*V*)βωf2+ω2M(a21+a11)(a21a22)(ω2ω02a21a22a21+a11)(1+2β2(f2+ω2)[(ω02ω2)2+4ω2ω02ζ2]M2(a21+a11)(a21a22)(ω2ω02a21a22a21+a11)).

For the first order approximation Yp(t)ZR0+(ZR1e2πit̃+ZR1e2πit̃) and based on definition of the runoff coefficient RC=RC*+(y1+y2)/f, we get

RCRC*+LM(a11+a21)(a22a21)(a11+a22)fβωf2+ω2(ω2ω02a21a22a21+a11)+L2ω2+(a22a21)2f1+(ω02ω22ωω0ζ)2sin(ωt+ψ+θ)+higherorderharmonics(k3),
(A2)

where

tanθ=(a22a21)+ω02ω22ω0ζω(a22a21)ω02ω22ωω0ζ.

In this case, the average value of the runoff coefficient is shifted from the equilibrium value RC* with a quantity depending on ω—see the second term (non-oscillating term) in the sum (A2). Using the definition of L, it results that the average value of the runoff coefficient is

RCavg(ω)=RC*+1RC*1+2β2(f2+ω2)[(ω02ω2)2+4ω2ω02ζ2]M2(a21a22)(a21+a11)(ω2ω02a21a22a21+a11),

or, written in an equivalent form, RCavg(ω)=RC*+Ω(ω)1+Ω(ω) with

Ω(ω)=M2(a21a22)(a21+a11)(ω2ω02a21a22a21+a11)2β2(f2+ω2)[(ω02ω2)2+4ω2ω02ζ2].
(A3)

Based on Eqs. (13) and (17), we note that a21+a11=P¯fβ and a21a22=(a21+a11)tr(J)=P¯fβ+2ω0ζ, so Ω(ω) becomes

Ω(ω)=M2P¯2(12ω0ζfβP¯)[ω2ω02(12ω0ζfβP¯)]2f2β4ω02(f2+ω2)[(ω02ω2)2+4ω2ω02ζ2].

Then Eq. (A2) is equivalent to

RCRCavg(ω)+A(ω)sin(ωt+ψ+θ)+higherorderharmonics(k3),
(A4)

with

A(ω)=M(1RC*)β[1+Ω(ω)]ω2+(P¯fβ2ω0ζ)2(f2+ω2)[(ω02ω2)2+4ω2ω02ζ2],

and

tanθ=P¯fβ2ω0ζ+ω02ω22ω0ζω(P¯fβ2ω0ζ)ω02ω22ωω0ζ.
(A5)
1.
V. M.
Castillo
,
A.
Gomez-Plaza
, and
M.
Martinez-Mena
, “
The role of antecedent soil water content in the runoff response of semiarid catchments: A simulation approach
,”
J. Hydrol.
284
,
114
130
(
2003
).
2.
L. K.
Cunha
, “
Exploring the benefits of satellite remote sensing for flood prediction across scales
,” Ph.D. thesis (
Department of Civil and Environmental Engineering, The University of Iowa
,
2012
), see http://ir.uiowa.edu/cqi/viewcontent.cgi?article=321&&context=etd.
3.
L. K.
Cunha
,
P. V.
Mandapaka
,
W. F.
Krajewski
,
R.
Mantilla
, and
A. A.
Bradley
, “
Impact of radar-rainfall error structure on estimated flood magnitude across scales: An investigation based on a parsimonious distributed hydrological model
,”
Water Resour. Res.
48
,
W10515
, doi: (
2012
).
4.
R.
Curtu
,
R.
Mantilla
,
M.
Fonley
,
L. K.
Cunha
,
S. J.
Small
,
L. O.
Jay
, and
W. F.
Krajewski
, “
An integral-balance nonlinear model to simulate changes in soil moisture, groundwater, and surface runoff dynamics at the hillslope scale
,”
Adv. Water Resour.
71
,
125
139
(
2014
).
5.
C.
Fitzjohn
,
J. L.
Ternan
, and
A. G.
Williams
, “
Soil moisture variability in a semi-arid gully catchment: Implications for runoff and erosion control
,”
Catena
32
,
55
70
(
1998
).
6.
See http://criticalzone.org/shale-hills/ for data and research from the Shale Hills Critical Zone Observatory in Pennsylvania
2012
.
7.
P.
Javelle
,
C.
Fouchier
,
P.
Armaud
, and
J.
Lavabre
, “
Flash flood warning at ungauged locations using radar rainfall and antecedent soil moisture estimations
,”
J. Hydrol.
394
,
267
274
(
2010
).
8.
J. A.
Lynch
and
E. S.
Corbett
, “
Source area variability during peakflow: A function of antecedent soil moisture content
,” in
Watershed Management in the Eighties Conference
, April 30–May 1, 1985, Denver, CO, edited by
E. B.
Jones
and
T. J.
Ward
(
American Society of Civil Engineers
,
New York
,
1985
), pp.
300
307
.
9.
R.
Mantilla
,
V. K.
Gupta
, and
O. J.
Mesa
, “
Role of coupled flow dynamics and real network structures on Hortonian scaling of peak flows
,”
J. Hydrol.
322
,
155
167
(
2006
).
10.
D.
Mackenzie
and
S.
Masten
,
Principles of Environmental Engineering and Science
(
McGraw-Hill Science/Engineering/Math
,
2003
).
11.
H. E.
Meadows
, “
Solution of systems of linear ordinary differential equations with periodic coefficients
,”
Bell Syst. Tech. J.
41
,
1275
1294
(
1962
).
12.
J. D.
Meiss
,
Differential Dynamical Systems
(
SIAM
,
Philadelphia
,
2007
).
13.
R.
Merz
and
G.
Bloeschl
, “
A regional analysis of event runoff coefficients with respect to climate and catchment characteristics in Austria
,”
Water Resour. Res.
45
,
W01405
, doi: (
2009
).
14.
S. H.
Strogatz
,
Nonlinear Dynamics and Chaos
(
Perseus Books Publishing
,
Cambridge, MA
,
1994
).
15.
E.
Zehe
and
G.
Bloeschl
, “
Predictability of hydrologic response at the plot and catchment scales: Role of initial conditions
,”
Water Resour. Res.
40
,
W10202
, doi: (
2004
).
16.
Y.
Zhang
,
H.
Wei
, and
M. A.
Nearing
, “
Effects of antecedent soil moisture on runoff modeling in small semiarid watersheds of southeastern Arizona
,”
Hydrol. Earth Syst. Sci.
8
,
6227
(
2011
).