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.

## I. INTRODUCTION

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 runoff^{6,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.

## II. DESCRIPTION OF MODEL

In this paper, we consider an ODE-based model for the local hillslope-link coupling^{4} 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=vw\u2212ponded/AH$, the volume per unit hillslope area of water stored in the ground surface (measured in meters); $v=\theta vunsat/AH$, the volume per unit hillslope area of moisture content in the unsaturated zone (measured in meters); *a* = *v _{sat}*/

*A*, the volume per unit hillslope area of the saturated zone (measured in meters); and

_{H}*q*=

*Q*/

*Q*, the (nondimensional) flow going downstream out of the channel at time $t\u0302$, where

_{r}*Q*is the runoff at the outlet (in m

^{3}/s) and

*Q*is the unit reference discharge (

_{r}*Q*= 1 m

_{r}^{3}/s). Time $t\u0302$ 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

and

Here, $E(t\u0302)$ represents the evapo-transpiration percentage loss function from ponded water; $p(t\u0302)$ is the precipitation input in [mm/h]; and parameters *c*_{1}, *c*_{2}, *c*_{3}, *c _{evap}*,

*γ*, and

*τ*are defined by $c1=Ksp60\u2009hb\u2009(m\u22121min\u22121),\u2009c2=10\u22126\u2009\alpha soil\u2009KSAT\u2009(2L)60\u2009AH\u2009(min\u22121),\u2009c3=10\u2212360\u2009(no\u2009unit)$, $cevap=Kevap60\u2009(min\u22121),\u2009\gamma =106\u2009AH60\u2009Qr\u2009(m\u22121\u2009min)$, and $\tau =(1\u2212\lambda 1)L60\u2009vr\u2009(Aupstream/Ar)\lambda 2\u2009(min)$. The incoming flux

*q*is taken here to be zero, given that the Shale Hills watershed is a first-order drainage basin, and the hillslope reference area

_{in}*A*= 1 km

_{r}^{2}is a normalization constant. The assumption is that

*h*represents an effective soil depth over the total hillslope area

_{b}*A*such that the total volume of the hillslope is

_{H}*V*=

_{tot}*h*. Therefore, variables

_{b}A_{H}*a*and

*v*are restricted to take values only in certain ranges, that is 0 ≤

*a*≤

*h*and 0 ≤

_{b}*v*≤

*h*−

_{b}*a*. Moreover, given the high density of tree roots in the Shale Hills watershed, the total volume of the hillslope effectively available for water (

*V*) is only a percentage of

_{T}*V*; this observation led to the introduction of a correction factor

_{tot}*β*such that

*V*is defined by $VT=\beta Vtot=\beta hbAH$ with $\beta \u2208(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.

_{T}Definition . | Physical interpretation . |
---|---|

$Qpl=c1\u2009sp(a\u2212ares+v\u2212vres)$ | Flux from ponded water to the channel link, per unit hillslope area |

$Qpu=c1\u2009sp(hb\u2212a\u2212v)$ | Flux from ponded water to the unsaturated zone, per unit hillslope area |

$Qus=d0(v\u2212vres)+$$d1(v\u2212vres)(a\u2212ares)2+d2(a\u2212ares)2$ | Flux exchange between the unsaturated and saturated zones, per unit hillslope area |

$Qsl=c2\u2009(a\u2212ares)\u2009\u2009exp(\alpha N\u2009a\u2212areshb)$ | Flux for subsurface runoff into the channel link, per unit hillslope area |

$Qevap=cevap\u2009(a\u2212ares)$ | Flux for the potential loss of groundwater through either evaporation or plant consumption, per unit hillslope area |

Definition . | Physical interpretation . |
---|---|

$Qpl=c1\u2009sp(a\u2212ares+v\u2212vres)$ | Flux from ponded water to the channel link, per unit hillslope area |

$Qpu=c1\u2009sp(hb\u2212a\u2212v)$ | Flux from ponded water to the unsaturated zone, per unit hillslope area |

$Qus=d0(v\u2212vres)+$$d1(v\u2212vres)(a\u2212ares)2+d2(a\u2212ares)2$ | Flux exchange between the unsaturated and saturated zones, per unit hillslope area |

$Qsl=c2\u2009(a\u2212ares)\u2009\u2009exp(\alpha N\u2009a\u2212areshb)$ | Flux for subsurface runoff into the channel link, per unit hillslope area |

$Qevap=cevap\u2009(a\u2212ares)$ | Flux for the potential loss of groundwater through either evaporation or plant consumption, per unit hillslope area |

Parameter . | Value . | Unit . | Physical interpretation . |
---|---|---|---|

A _{upstream} | 0.07736 | km^{2} | Total upstream area |

A _{H} | 0.07736 | km^{2} | Total hillslope area |

L | 420.15 | m | Channel (link) length |

h _{b} | 0.56 | m | 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 | |

v _{r} | 0.25 | m/s | Reference flow velocity |

K _{SAT} | 0.01 | m/h | Saturated hydraulic conductivity |

d_{0} | 0.00135 | min^{−1} | d_{0}, d_{1}, and d_{2} are parameters that describe the formof the hillslope-scale recharge relation coupling the statevariables |

d_{1} | 0.013 | $m\u22122min\u22121$ | |

d_{2} | 0.008 | $m\u22121min\u22121$ | |

a _{res} | 0.1 | m | a and _{res}v are the residual storage volumes for thegravity-drained hillslope _{res} |

v _{res} | 0.3 | m | |

β | 0.12 | Parameter that accounts for the heterogeneities in the soilmatrix porosity and areas inaccessible to water | |

α _{N} | 6 | Parameter that controls the recession | |

α _{soil} | 3.0 | Heterogeneity factor for soil saturated hydraulic conductivity | |

$Ksp$ | 1.75 | $h\u22121$ | Parameter that characterizes the overland flow to the channeland the infiltration to the soil |

K _{evap} | 0.025 | $h\u22121$ | Recession coefficient for the evaporation process andwater consumption by vegetation, depending on airtemperature, humidity and species of plants on the hillslope |

Parameter . | Value . | Unit . | Physical interpretation . |
---|---|---|---|

A _{upstream} | 0.07736 | km^{2} | Total upstream area |

A _{H} | 0.07736 | km^{2} | Total hillslope area |

L | 420.15 | m | Channel (link) length |

h _{b} | 0.56 | m | 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 | |

v _{r} | 0.25 | m/s | Reference flow velocity |

K _{SAT} | 0.01 | m/h | Saturated hydraulic conductivity |

d_{0} | 0.00135 | min^{−1} | d_{0}, d_{1}, and d_{2} are parameters that describe the formof the hillslope-scale recharge relation coupling the statevariables |

d_{1} | 0.013 | $m\u22122min\u22121$ | |

d_{2} | 0.008 | $m\u22121min\u22121$ | |

a _{res} | 0.1 | m | a and _{res}v are the residual storage volumes for thegravity-drained hillslope _{res} |

v _{res} | 0.3 | m | |

β | 0.12 | Parameter that accounts for the heterogeneities in the soilmatrix porosity and areas inaccessible to water | |

α _{N} | 6 | Parameter that controls the recession | |

α _{soil} | 3.0 | Heterogeneity factor for soil saturated hydraulic conductivity | |

$Ksp$ | 1.75 | $h\u22121$ | Parameter that characterizes the overland flow to the channeland the infiltration to the soil |

K _{evap} | 0.025 | $h\u22121$ | 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\u0302)$ 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\u0302)$ is multiplied by the factor $(1\u2212E(t\u0302))$. Contrarily, *Q _{evap}* 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

*Q*cannot take negative values as it is assumed that

_{evap}*a*is always above

*a*(see below for more explanations). These two methods of water loss are typically combined into the term “evapotranspiration.”

_{res}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 *a _{res}* and

*v*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 (

_{res}*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,

*a*,

_{res}*v*can be rather assumed at

_{res}*p*(

*t*) = 0; see Ref. 4 for an additional discussion and comparison with hydrological data.

In Secs. III–V, we will replace the term $[1\u2212E(t\u0302)]p(t\u0302)$ from (1) with the “effective precipitation” function $peff(t\u0302)$, 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 *d*_{0}, *d*_{1}, and *d*_{2} 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 *Q _{pu}*) or flows along the hillslope as surface runoff (through

*Q*). Then the runoff coefficient is defined by

_{pl}where *f* is the percentage of available hillslope storage volume calculated at the residual (gravity-drained hillslope) condition

## III. SYSTEM NONDIMENSIONALIZATION AND EQUILIBRIA

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\u0302)$. (In this study, we will keep fixed the value of parameter *K _{evap}* 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 *p _{eff}* and of the variables

*s*,

_{p}*v*,

*a*, and time $t\u0302$

leads to the nondimensional system

together with the equation for the runoff coefficient

The parameters *α* and *α _{evap}* and the flux $Q\u0303us$ are defined now by

For example, using the Definition *Q _{us}*(

*a*,

*v*) from Table I, we will work with the following formula for $Q\u0303us(A,V)$:

Note also that $Q\u0303us(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 $ares\u2264a\u2264a+v\u2264hb$. Therefore, *A* and *V* are restricted to the following range:

### A. Characterization of the equilibrium point

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\xaf$. From physical point of view, this might seem to be an unrealistic condition. However, it is reasonably valid if we interpret $P\xaf$ in a more general context; that is, we define $P\xaf$ as the averaged value of precipitation during a given interval of time, *T _{observe}* (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

*T*, taking the rain average value and using it as constant input to the system is an acceptable assumption.

_{observe}At constant (average) rain input of $P(t)=P\xaf$ with $P\xaf>0$, the system (5) reaches an equilibrium $(S*,A*,V*)$ that satisfies the following conditions: $S*=P\xaf/f,\u20090=S*\u2009(f\u2212A*\u2212V*)\u2212Q\u0303us(A*,V*)$ and $0=Q\u0303us(A*,V*)\u2212\alpha \u2009A*\u2009e\alpha N\u2009A*\u2212\alpha evap\u2009A*$. These allow for an explicit calculation of $V*$ according to the formula

and an implicit definition for $A*$ as a function of $P\xaf$ by

As already mentioned

We exclude here the case $P\xaf=0$. In that case, obviously, the only possible equilibrium in system (5) is (0, 0, 0) which corresponds to *s _{p}* = 0,

*v*=

*v*, and

_{res}*a*=

*a*.

_{res}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\xaf$, 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\xafeff$ is shown in Figure 1.

#### 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\u0303us(A,V)$ is defined by (8). The curves $(P\xaf,A*)$ and $(P\xaf,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\xafeff,a)$ and $(p\xafeff,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).

### B. How could the rain pattern prime the hillslope to the “pre-flooding” condition? The system's natural frequency and damping ratio

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 $(\lambda 2=\lambda ;\u2009\lambda 3=\lambda \xaf;\u2009Re\u2009(\lambda )<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*)$

where

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 $\lambda 1=\u2212f$ and *λ*_{2,3}, the roots of the quadratic equation $\lambda 2\u2212tr(J)\u2009\lambda +det(J)=0$, where

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)2\u22124det(J)<0$.

Under the circumstances described above ($\lambda 1<0$ and $\lambda 2,3\u2208C\R,\u2009Re\u2009(\lambda 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 *T _{observe}*. 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 $\lambda 2=\lambda $ and $\lambda 3=\lambda \xaf$; then

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 $\zeta <1/2\u22480.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 *λ*, $\lambda \xaf$ 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\xaf$ for precipitation that lead to complex eigenvalues $\lambda 2,3=\lambda ,\lambda \xaf$ 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\xaf$ falls in the range defined by

Moreover, since $tr(J)=2\u2009Re\u2009(\lambda )$ and $det(J)=|\lambda |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

and then $\lambda =\u2212\omega 0\zeta +i\omega 01\u2212\zeta 2$. Figure 2 illustrates the range of averaged effective precipitation (*p _{eff}* > 0.4 [mm/h]) for which condition (16) is valid, and the graphs of the damping ratio

*ζ*and the natural frequency $\omega \u03020=\omega 0\u2009Ksp$ in $[h\u22121]$ and period $T\u03020=2\pi /\omega \u03020$ in [h] versus $p\xafeff$. (Note that for

*p*< 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

_{eff}*ω*

_{0}nor

*ζ*make sense in this interval.)

### C. Estimation of the time to the equilibrium

We define the *time to the asymptotic solution with respect to the ponded water s _{p}* =

*h*as the minimum time necessary for

_{b}S*s*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

_{p}where *S*_{0} is any initial condition and $S*=P\xaf/f$ is the equilibrium value at $P(t)=P\xaf$. Therefore, the asymptotic time to the steady state $sp*$ can be computed according to the inequality

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

Consequently, the time to equilibrium *t _{eq}* is estimated by $teq=max(teq,S;teq,AV)$

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

An example that estimates the time to equilibrium $t\u0302eq$ when $p\xafeff$ varies is shown in Figure 3.

We calculate, similar to (19), and plot the time $t\u0302eq$ (in hours) required for the system to come within 1%, 5%, and 10% of its equilibrium value. For example, at $p\xafeff=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 *s*_{0} = 0, *a*_{0} = *a _{res}*, and

*v*

_{0}=

*v*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\omega 0\zeta $ and $1f$ determine which term will determine the maximum. For smaller values of average precipitation, the third and fourth expressions are larger because $1/\omega 0\zeta >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).

_{res}## IV. DYNAMICS OF THE RUNOFF COEFFICIENT UNDER FLUCTUATING RAIN PATTERNS

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 $[1\u2212E(t\u0302)]p(t\u0302)$ replaced by the effective precipitation function $peff(t\u0302)$. The fluctuations about the average value $P\xaf$ are assumed to have amplitude *M* with $0\u2264M\u2264P\xaf$ in order to always have $peff(t\u0302)\u22650$, so

As a reminder from (4), time $t\u0302$ is measured in minutes, *K _{sp}* is in $h\u22121$ and $t=t\u0302\u2009Ksp60$ is the nondimensional time. Then, obviously, the nondimensional input to system (5) is

Obviously, while the effective precipitation pattern is fluctuating, on average, its value is still $P\xaf$. We are interested to determine if there exits any particular choice of amplitude $M=M(P\xaf)$ and frequency $\omega =\omega (P\xaf)$ 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\xaf$ and the rainfall is oscillatory. Then, during the oscillation, the system receives precipitation input with values between 0 and $2P\xaf$. The steady state analysis predicts certain values for the runoff coefficients for each of $P\u2208[0,2P\xaf]$ (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\xaf$? 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 *p _{eff}* = 10 mm/h, the steady state $RC*$ is $RC*(10)=0.497$ while at

*p*= 0 and

_{eff}*p*= 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

_{eff}*p*varies between 0 and 20 mm/h in an organized (oscillatory) fashion, with certain frequency

_{eff}*ω*,

*RC*can reach significantly larger values. Thus,

*RC*may take values up to 0.8; see Figure 4. By comparison, in the steady state condition,

_{peak}*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.

Numerical simulations of the nonlinear model (1) with parameters from Table II indicate a nonlinear dependence of the *RC _{peak}* on the frequency of oscillation

*ω*(Figure 4). In this example, the average rainfall is

*p*= 10 mm/h. The peak of

_{eff}*RC*increases with

*ω*until it reaches a maximum at a certain

*ω*; then for larger frequencies

_{R}*ω*>

*ω*, the peak decreases with

_{R}*ω*.

To investigate the nonlinear relationship $RCpeak=RCpeak(\omega )$, 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

where *S*_{0} is any initial condition, $S*=P\xaf/f$ is the equilibrium value at $P(t)=P\xaf$, and $Sf(t)$ is the fluctuating (and non-damping) term

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

By assuming that a rain input of average magnitude $P\xaf$ but without any particular structure was applied in (5) for a time period at least equal to *t _{eq}*, 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\xaf$ 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\xaf+0.01P\xaf)/f$ and it is $tfluct,S=max(\u20090\u2009;\u20091fln|1+f\omega f2+\omega 2\u2009M0.01P\xaf|\u2009)$. 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

with $t\u0302eq$ 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 *t _{transient}* has passed (

*T*>

_{observe}*t*), the input to equation

_{transient}*V*is $S(t)\u2248P\xaf/f+Sf(t)$ and system (5) reduces to

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.

### A. Nonlinear response of the runoff coefficient to the frequency of the rain pattern

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

with linearization

where coefficients *a _{ij}* are defined by (13). Therefore, in first order approximation, we have $A(t)\u2248A*+y2(t)$ and $V(t)\u2248V*+y1(t)$. The runoff coefficient $RC=(A+V)/f$ is then approximated by $RC(t)\u2248RC*+(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=1T\u222b0TRC(t)dt$ with $T=2\pi \omega $ as a function of the frequency *ω* of the rain pattern *P*(*t*) (see Subsection 3 of the Appendix for more details). This is

with

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

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

Interestingly, the hillslope response (defined here as the peak of the runoff coefficient, *RC _{peak}*) to the fluctuating rain pattern $P(t)=P\xaf+M\u2009sin(\omega t)$ with same average value $P\xaf$, depends nonlinearly on the choice of frequency

*ω*. There exists a particular value

*ω*where the peak response is maximal. In fact, as shown in Figure 5, at

_{R}*ω*the peak response $RCmax=RCpeak(\omega 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\xafeff$ or $meff=p\xafeff$; Figure 5 (center) and (top)),

_{R}*RC*reaches at

_{max}*ω*even more significant levels. These show approximately 40% and 60% increases compared to nearby values of

_{R}*ω*or, to the steady state value $RC*$ at $p\xafeff$ ($RC*=0.497$ at $p\xafeff=10\u2009mm/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,

*ω*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

_{R}*ω*close to

*ω*is what we call a

_{R}*pre-flooding condition*.

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\xafeff=M/P\xaf=\u03f5$ approaches 1. This can be explained by the fact that mathematically, linear approximations to nonlinear system yield best results when $\u03f5\u21920$ given that $S(t)\u2248S*+Sf(t)=S*(1+O(\u03f5))$. 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(\omega )+A(\omega )$, to be safely omitted in *RC*. In fact, numerical comparisons do indicate an increased error close to the resonant frequency *ω _{R}* between the peak

*RC*of the runoff coefficient, and its corresponding value computed from the nonlinear system (see Figure 5). On the contrary, further away from

_{peak}*ω*and for small values of frequencies

_{R}*ω,*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\pi \omega R$. The unit-based counter-parts of the

*resonant frequency*and

*period*are then determined by the equations $\omega \u0302R=Ksp\u2009\omega R$ in $[h\u22121]$ and $T\u0302R=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

*ω*.

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\xafeff=10\u2009[mm/h]$ and its resonant frequency $\omega \u0302R=0.4242\u2009[1/h]$. The trajectory oscillates around the equilibrium point $(v*,a*)=(0.2,0.3)$ at $p\xafeff=10\u2009[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\xafeff$ 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.

## V. CONCLUSION

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}*,

*α*, and

_{soil}*K*, or coefficients

_{sp}*d*

_{0},

*d*

_{1},

*d*

_{2}from the definition of the flux

*Q*(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

_{us}*RC*is rather driven by the magnitude of the rain storm. As mentioned, by changing parameters such as

*α*,

_{N}*α*,

_{soil}*K*,

_{sp}*d*

_{0},

*d*

_{1}, or

*d*

_{2}, 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,

*K*and

_{sp}*α*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,

_{soil}*α*and

_{N}*d*

_{0},

*d*

_{1},

*d*

_{2}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.

## ACKNOWLEDGMENTS

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.

### APPENDIX: ANALYSIS OF THE LINEAR SYSTEM WITH PERIODIC COEFFICIENTS

##### 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\u0303=1$. The time-change is defined by $\omega t+\psi =:2\pi t\u0303$. Note also that in complex plane (here $i=\u22121$), $sin\u20092\pi t\u0303=12i(e2\pi i\u2009t\u0303\u2212e\u22122\pi i\u2009t\u0303)$.

The stationary matrix $(alj)$ has distinct complex conjugated eigenvalues *λ*_{1} and *λ*_{2} which are the roots of the characteristic equation $\lambda 2+2\zeta \omega 0\lambda +\omega 02=0$ with $\zeta ,\omega 0$ defined by (17), say, $\lambda 1=\u2212\zeta \omega 0+i\omega 01\u2212\zeta 2,\u2009\lambda 2=\lambda \xaf1$. We consider the matrix *Z* whose columns are the corresponding eigenvectors, and its inverse $Z\u22121$

and make the change of variables $X=(x1,x2)T=Z\u22121(y1,y2)T=Z\u22121Y$. (Here, *T* stands for the transpose.) This leads to the system

with *a*, *b*, *c, d* defined by $\Lambda 1=M\pi \beta \omega f2+\omega 2\u2009ia21det(Z),\u2009\Lambda 2=M\pi \beta \omega f2+\omega 2\u2009i(a11\u2212\lambda 1)det(Z),\u2009a=\Lambda 1\Lambda 2\u2009c,\u2009b=\Lambda 1\Lambda 2\u2009d,\u2009c=(\lambda 1\u2212a11+a12)\Lambda 2,\u2009d=\u2212(\lambda 1\u2212a11\u2212a21)\Lambda 2,\u2009e=\u2212(f\u2212A*\u2212V*)\u2009\Lambda 1$ and $g=\u2212(f\u2212A*\u2212V*)\u2009\Lambda 2$.

The order of magnitude of the above coefficients can be estimated with respect to parameter *β*, assumed to be small (e.g., $\beta =0.12$ as considered in Sec. II). Equations (13), (14), and (17) imply that $\beta alj=O(1)$ and $\beta \omega 0=O(1)$ while $\zeta =O(1)$; therefore $\beta \lambda 1,2=\u2212\zeta \beta \omega 0\xb1i\u2009\beta \omega 01\u2212\zeta 2=O(1)$ and $\beta 2det(Z)=O(1)$.

Then we have $\Lambda 1,\Lambda 2=2\pi \omega \u2009O(1)$ and $a,b,c,d=2\pi \omega \u2009O(\beta \u22121)$, and $e,g=2\pi \omega \u2009O(1)$. Consequently, the matrices *B*_{0}, *B*_{1}, and *F*_{0} have different order of magnitude with respect to *β* as they are defined by $B0=2\pi \beta \u22121\omega \u2009O(1),\u2009B1=2\pi \beta \u22121\omega \u2009O(1),\u2009F0=2\pi \beta \u22121\omega \u2009O(\beta )$ as $\beta \u226a1$.

##### 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\u0303)=C1e\mu 1t\u0303\u2009Q\mu 1(t\u0303)+C2e\mu 2t\u0303\u2009Q\mu 2(t\u0303)$, where *C*_{1}, *C*_{2} are constants that depend on the initial condition and $e\mu 1t\u0303\u2009Q\mu 1(t\u0303),\u2009e\mu 2t\u0303\u2009Q\mu 2(t\u0303)$ are two independent (nontrivial) solutions of the form

Indices *k* are integers. With notation *I*_{2} and *O*_{2} for the two-dimensional unity and null matrices, *Q _{k}* are two-dimensional constant vectors from the complex plane

^{11}satisfying

Then the solution of (A1) is the superposition of two functions: *X _{homog}* 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\mu 1,2t\u0303$ in the definition of

*X*shrink with time such that this part of the general solution is transient. This is equivalent to proving that $Re(\mu 1,2)<0$. However, unfortunately, $Re(\mu 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.

_{homog}According to the theory of approximating solutions for linear systems of differential equations with periodic coefficients,^{11} $z1=e\mu 1$ and $z2=e\mu 2$ are the two roots of the equation

with notation $r1=e2\pi \u2009\lambda 1\omega $ and $r2=e2\pi \u2009\lambda 2\omega $ and *K*_{1}, *K*_{2} defined by some complicated infinite determinants that are of the same order as *a*, *b* (for *K*_{1}), respectively, *c*, *d* (for *K*_{2}). Given that $a,b,c,d=2\pi \beta \u22121\omega \u2009O(1)$, we have $K1,K2=2\pi \beta \u22121\omega \u2009O(1)$ as well.

Note that the quadratic equation in *z* is equivalent to $(\u20091+K1+K22\u2009)z2\u2212(r1+r2)(\u20091\u2212(K1\u2212K2)(r1\u2212r2)2(r1+r2)\u2009)z+r1r2(\u20091\u2212K1+K22\u2009)=0$ and we have $r1,r2=eO(\beta \u22121)$ and $K1,K2=O(\beta \u22121)$ so $r1,r2\u226bK1,K2$ as $\beta \u226a1$. Then the roots *z*_{1}, *z*_{2} satisfy $z1\u2248r1(1\u2212K1)$ and $z2\u2248r2(1\u2212K2)$ so $\mu 1=2\pi \lambda 1\omega +2\pi \omega \u2009Re(\lambda 1,2)\u2009O(\beta )$ and $\mu 2=2\pi \lambda 2\omega +2\pi \omega \u2009Re(\lambda 1,2)\u2009O(\beta )$ with $Re(\lambda 1,2)=O(\beta \u22121)$.

As a consequence, the asymptotic behavior of $Xhomog(t\u0303)$ is prescribed by the limit of $eRe(\mu 1)t\u0303$ and $eRe(\mu 1)t\u0303$. In original time *t* ($\omega t+\psi =:2\pi t\u0303$) and coordinates *Y*, the dynamics of the transient term of the solution is given by the limit of $e\omega 2\pi Re(\mu 1)t=e(1+O(\beta ))\u2009Re(\lambda 1)t$ and $e\omega 2\pi Re(\mu 2)t=e(1+O(\beta ))\u2009Re(\lambda 2)t$, where *β* is small. Both exponentials tend to zero given that the dominant exponent is $Re(\lambda 1,2)=\u2212\zeta \omega 0=O(\beta \u22121)<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\u0303)$. If this can be found then the general solution of (A1) is $X(t\u0303)=Xhomog(t\u0303)+Xp(t\u0303)=C1e\mu 1t\u0303\u2009Q\mu 1(t\u0303)+C2e\mu 2t\u0303\u2009Q\mu 2(t\u0303)+Xp(t\u0303)$ and so $X(t\u0303)\u2192Xp(t\u0303)$ as $t\u0303\u2192\u221e$. The steady, non-decaying, oscillatory function $Xp(t\u0303)$ is a good approximation of the solution.

Given the form of Eq. (A1), we are searching for a function $Xp(t\u0303)$ such that $Xp(t\u0303)=\u2211k=\u2212\u221e\u221eRk\u2009e2\pi ik\u2009t\u0303$ with index *k* spanning the set of integers. Introducing *X _{p}* in (A1), it results that the two-dimensional vectors

*R*should satisfy

_{k}By considering the complex inner product between two dimensional vectors $\u27e8U,V\u27e9=(u1,u2)T\xb7(v1,v2)T=u1v\xaf1+u2v\xaf2$, we determine first the solution of the adjoint equation $B1*F=0$ and choose $F=(\lambda \xaf1\u2212a11,\u2009a21)T$ accordingly. Then, since $\u27e8B1(Rk+1\u2212Rk\u22121),F\u27e9=\u27e8Rk+1\u2212Rk\u22121,B1*F\u27e9=0$ and $\u27e8F0,F\u27e9=0$, we note that all vectors *R _{k}* satisfy the condition $\u27e8(B0\u22122\pi kiI)Rk,F\u27e9=0$. Moreover, since we are looking for a real solution $Yp(t)=ZXp(t)=\u2211k=\u2212\u221e\u221eZRk\u2009e2\pi ik\u2009t\u0303=ZR0+\u2211k=1\u221e(ZRk\u2009e2\pi ik\u2009t\u0303+ZR\u2212k\u2009e\u22122\pi ik\u2009t\u0303)$ with matrix

*Z*defined in Subsection 1 of the Appendix, we need

*ZR*

_{0}to be real and $ZR\u2212k=ZRk\xaf$ for $k\u22651$. That determines vectors $Rk=(v1,k,v2,k)T$ and $R\u2212k=(v1,\u2212k,v2,\u2212k)T$ with coordinates satisfying $v1,k=1detZa21(\lambda 2\u2212ik\omega )rk,\u2009v2,k=1detZ(a11\u2212\lambda 1)(\lambda 1\u2212ik\omega )rk$ and $v1,\u2212k=v1,k(\lambda 2+ik\omega )r\xafk(\lambda 2\u2212ik\omega )rk,\u2009v2,\u2212k=v2,k(\lambda 1+ik\omega )r\xafk(\lambda 1\u2212ik\omega )rk$, for some complex constant

*r*, and

_{k}However, $ZB0Z\u22121=2\pi \omega (a11a12a21a22),\u2009ZB1Z\u22121=i\pi M\omega \beta f2+\omega 2(1100),\u2009ZF0=\u2212i\pi M(f\u2212A*\u2212V*)\omega \beta f2+\omega 2(10)$. Since $(ZB0Z\u22121)(ZR0)=ZB0R0=ZB1(R1\u2212R\u22121)=(ZB1Z\u22121)(ZR1\u2212ZR\u22121)$ and with notation $r1=L+iH$, we obtain

Let us observe that we can determine *R*_{1} and $R\u22121$ by choosing $R2=R\u22122=0$. This is true because $(ZB0Z\u22121\u22122\pi iI)ZR1=Z(B0\u22122\pi iI)R1=\u2212(ZB1Z\u22121)ZR0\u2212ZF0$ and similar $(ZB0Z\u22121\u22122\pi iI)ZR1\xaf=(ZB0Z\u22121+2\pi iI)ZR\u22121=Z(B0+2\pi iI)R\u22121=(ZB1Z\u22121)ZR0+ZF0=\u2212(ZB1Z\u22121)ZR0\xaf\u2212ZF0\xaf$. These allow us to solve for *L* and *H* from the equation

which implies

and

For the first order approximation $Yp(t)\u2248ZR0+(ZR1\u2009e2\pi i\u2009t\u0303+ZR\u22121\u2009e\u22122\pi i\u2009t\u0303)$ and based on definition of the runoff coefficient $RC=RC*+(y1+y2)/f$, we get

where

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

or, written in an equivalent form, $RCavg(\omega )=RC*+\Omega (\omega )1+\Omega (\omega )$ with

Based on Eqs. (13) and (17), we note that $a21+a11=\u2212P\xaff\beta $ and $a21\u2212a22=(a21+a11)\u2212tr(J)=\u2212P\xaff\beta +2\omega 0\zeta $, so $\Omega (\omega )$ becomes

Then Eq. (A2) is equivalent to

with

and