T Cheviron

A Chriette

F Plestan

Generic Nonlinear Model of Reduced Scale UAVs

This paper proposes, through a survey of models of several UAV-Structures, a generic nonlinear model for reduced scale aerial robotic vehicles (6 DOF) 1 . Dynamics of an aircraft and some VTOL UAV (quadricopter, ducted fan and classical helicopter) are illustrated. This generic model focuses only on the key physical efforts acting on the dynamics in order to be sufficiently simple to design a controller. The Small Body Forces expression which can introduce a zero dynamics is then discussed.

I. INTRODUCTION

Usually, UAVs can be classified into close, short and long ranges according to their areas of mission. There exist numerous types of UAVs to fit for various demands about altitude, range and duration, in the same way than payload capability, volume capability and control capability.

With a range of applications in both civilian and military scenarios, the development of automated aerial robots are an increasingly important field of robotics research. Such vehicles have strong commercial potential in remote surveillance applications such as monitoring traffic congestion, regular inspection of infrastructure such as bridges, dam walls and power cables or investigation of hazardous environments, to name but a few of the possibilities. The development of such robotic vehicles states a number of problems in sensoring and control. A key challenge is to develop light aerial vehicles able of autonomous navigation.

To develop the flight control systems for maneuverable autonomous reduced scale aerial vehicles, dynamic models that are accurate for their flight envelope are needed. For example, in the case of a standard helicopter, dynamic models explicitly must take into account the effects such that the rotor/fuselage coupling. However, in order to design nonlinear control law, minimal complexity models are prefered.

The main difficulties (at a theoretical level) for design stable feedback controllers for such vehicles come from nonlinearities and couplings (for the solid mechanics part) and from the fact that inputs are not torques nor forces but displacements of some elements which enter the dynamics through aerodynamical forces/torques. Design, modeling and control of autonomous flying systems have now become very challenging areas of research, as shown by a large literature since 90's decade [START_REF] Munzinger | Development of a real-time flight simulator for an experimental model helicopter[END_REF], [START_REF] Chriette | Zero dynamics analysis for ibvs control of under-actuated rigid body dynamics[END_REF], [START_REF] Gravilets | Nonlinear model for a smallsize acrobatic helicopter[END_REF], [START_REF] Mettler | Identification modeling and characteristics of miniature rotorcraft[END_REF], [START_REF] Avila Vilchis | NonLinear modelling and control of helicopters[END_REF], [START_REF] Bouabdallah | Backstepping and sliding-mode techniques applied to an indoor micro quadrator[END_REF], [START_REF] Hauser | Non linear control design for slightly non minimum phase systems: Application to v/stol aircraft[END_REF], [START_REF] Koo | Output tracking control design of a helicopter model based on approximate linearization[END_REF], [START_REF] Mettler | Identification modeling and characteristics of miniature rotorcraft[END_REF] for small-size helicopter and [START_REF] Prouty | Helicopter, Performance, Stability, and Control[END_REF][13] [START_REF] Stepniewsky | Rotor wing aerodynamics, basic theories of rotor aerodynamics[END_REF] for nonlinear full-scale rotorcraft models. However, there is no work that has been made on the design of a general (generic) aerodynamic model valid for all autonomous flying system.

The main contribution of this paper is to present a generic nonlinear model of reduced scale UAVs. A panorama of dynamics of an aircraft and some VTOL UAV (quadricopter, ducted fan and classical helicopter) is then presented. A generic model focuses only on the key physical efforts acting on the dynamics in order to be simple sufficiently enough to design a controller. In addition, the Small Body Forces expression which can introduce a zero dynamics is then discussed in this general case.

The outline of the paper is as follows. Section II is dedicated to the Background of the study. In Section III different UAV structures are reviewed. The Generic 6DOF State Model and Control Strategy are the subjects of section IV. Finally, we present the conclusion of this work in section V.

II. PROBLEM STATEMENT

In this section, analytic expressions for the forces and moments on the rigid body are derived. The forces and moments are referred to a system of body-fixed axes, centered at the center of gravity.

There are in general two approaches in deriving equations of motion. One is to apply Newton's law and Euler's law which can give some physical insight through the derivation. The other one is more complicated, it provides the linkage between the classical framework and the Lagrangian or Hamiltonian framework.

In this paper, applying Newton's laws of motion relating the applied forces and moments to the resulting translational and rotational accelerations assembles the equations of motion for the 6 degrees of freedom.

We will make in the sequel some simplifying assumptions : the vehicle is rigid and the earth fixed reference frame is inertial, the gravitational field is constant, the aerial vehicle is supposed to be a rigid body, the density of air is supposed to be uniform.

A. Rigid body dynamics

We begin by define the reference frames are considered in the derivation of the kinematics and dynamical equations of motion 2 .

• Let I = {e i 1 , e i 2 , e i 3 } denote the frame attached to the earth can be assumed inertial. e i 3 denotes the downwards vertical direction and e i 1 points to the magnetic north.

• Let B = {e b
1 , e b 2 , e b 3 } be a body-fixed frame whose center coincides with the center of mass of the mobile.

• Let A = {e a
1 , e a 2 , e a 3 } denote the air-frame 3 attached to the body fixed frame, as the previous frame, its center coincides with the center of mass of the mobile.

• The attitude of the body-fixed frame is represented by a rotation matrix R : B → I. Let R : A → B denote the rotation matrix made up the angle of attack α and the sideslip angle β. The set of all rotation matrices is termed the Special Orthogonal group and denoted by SO(3). Let ξ i , v i , R and Ω b denote respectively the linear position and velocity of the center of mass, the attitude and the angular velocity of the aerial vehicle, i.e

ξ i = [x y z] T v i = [v x v y v z ] T R =   c θ c ψ s φ s θ c ψ -c φ s ψ c φ s θ c ψ + s φ s ψ c θ s ψ s φ s θ c ψ + c φ c ψ c φ s θ s ψ -s φ c ψ -s θ s φ c θ c φ c θ   Ω b = [p q r] T (1)
with η = [φ θ ψ] T is the Euler angles vector and c χ , s χ denoting, respectively, cos χ and sin χ .

Let m and J be the mass and the inertia matrix of the rigid object assumed to be constant. The motion of a rigid body subjected to forces and torques is then described by the Newton-Euler equations [START_REF] Goldstein | Classical mechanics 2nd ed[END_REF]

ξi = v i m vi = mge i 3 + RF b Ṙ = Rsk(Ω b ) J Ωb = -sk(Ω b )JΩ + C b (2)
where g ≈ 9.81 [m.s -2 ] the gravitational acceleration, F b and C b the resulting forces (excluding the gravity force) and torques acting on the rigid body.

The notation sk(Ω) denotes the skew-symmetric matrix such that sk(Ω)v = Ω × v for the vector cross-product and any vector v ∈ R.

The expressions of F b and C b are detailed in the next section for some classical UAV's architecture.

III. EXAMPLES OF UAV'S DYNAMICS

The complexity of a state model essentially depends on the expression of aerodynamic forces and torques. 2 In the sequel of the paper, a vector denoted w * means that w is expressed in {e *

1 , e * 2 , e * 3 } frame with * ∈ {i, a, b} 3 Also called Aerodynamic Reference Frame.

A. A reduced scale aircraft (Fig. 1)

In this case, the force F b is composed by F b A (produced by the airframe) and F b T (produced by the propeller thrust)

F b = F b A + F b T The propeller thrust F b T = F t e b 1 , F t > 0 generates a sufficient airspeed V t .
F b A = qS R   C X0 + C Xα α + C Xβ β C Y β β C Z0 + C Zα α   + ΣΓ Σ =   C Xδa C Xδe C Xδr C Y δa C Y δe C Y δr C Zδa C Zδe C Zδr   (3)
with C Z lift coefficient, C X and C Y logitudinal and lateral drag coefficients, q = 1 2 ρV 2 t the dynamic pressure (ρ is the air density), S the wing surface and Γ = [δ a δ e δ r ] T the control-surface deflection vector (aileron, elevator and rudder) used to generate the torque which turns the airplane the desired way.

C b A = qS   C Lβ β + C Lp p + C Lr r C M 0 + C M α α + C M q q C N β β + C N r r   + KΓ K =   C Lδa 0 0 0 C Lδe 0 0 0 C Lδr   (5)
where :

p = bp 2 Vt , q = cq 2 Vt , r = br 2 Vt (6)
with : b the wingspan (the distance from the left wingtip to the right wingtip), c the mean chord of wing.

B. VTOL vehicles

Different rotorcraft structures exist to perform hover or vertical take-off and landing. All these structures have, at least, one rotor such that its thrust directly confronts gravity force. By supposing that the vehicle has k rotors, the thrust induced by a rotor j (1 j k) may be written [START_REF] Prouty | Helicopter, Performance, Stability, and Control[END_REF] with ω j the angular velocity of the rotor, θ Tj the collective angle, α Tj the angle of attack of blades and the coefficient C Di the drag of the rotor. The constant C Mj = 1 4 R 3 j n j cj depends on the radius, the number and the mean chord of blades. Besides, the rotor induces a torque C b airj due to the air resistance and proportional to ω 2 j .

F Tj = ρ C Mj (θ Tj -α Tj ) -C Dj ω 2 j = b j ω 2 j (7)
To develop the flight control systems for maneuverable autonomous miniature helicopter, dynamic models that are accurate for their flight envelope are needed. However, in order to design nonlinear control law, minimal complexity models are preferred. This models, developped by using linear system identification for example, are only valid in the vicinity of the nominal operating point, for example at hover.

• Classical helicopter (Fig.

with a 1s and b 1s the longitudinal and lateral flapping angles assumed to be small. Under quasi-stationary flying conditions, these angles depend algebraically on cyclic angles (i.e.: mechanical control inputs) [START_REF] Prouty | Helicopter, Performance, Stability, and Control[END_REF], [START_REF] Mahony | Robust trajectory tracking for a scale model autonomous helicopter[END_REF].

The main rotor induces a torque C b mot acting on the fuselage which is compensated by the tail rotor. This smaller rotor controls the yaw motion which yields

F b Tt = F Tt e b 2 (9)
From ( 8) and ( 9

F b = -F Tm e b 3 + ΣΓ C b = C b air + C b damp + C b mot + F Tm K b 0 + KΓ Σ = LK -1 (11)
with

L =   1 0 0 0 1 1 0 0 0   , K b 0 = [l 2 Tm -l 1 Tm 0] T and K =   0 -l 3 Tm -l 3 Tt l 3 Tm 0 0 -l 2 Tm l 1 Tm l 1 Tt   , l i = 3 k=1 l k i E a
k (i stands for m, t) represents the distance between the center of gravity and the application point of the force. • Quadricopter (Fig. 3): This architecture composed by four identical rotors has some advantages with respect to conventional helicopters. In fact, gyroscopic effects and air resistance torques tend to cancel in trimmed flight because front and rear motors rotate counterclockwise while the other two rotate clockwise. Besides, the collective angle is constant because of the small radius of blades, the thrust of each rotor being then simply controlled by their angular velocities. The collective input is the sum of the thrusts of each rotor. Pitch motion is obtained by increasing (reducing) the speed of the rear motor while reducing (increasing) the speed of the front motor. Roll motion is obtained by a similar way by the lateral motors. The yaw movement is obtained by increasing (decreasing) the speed of the front and rear motors while decreasing (increasing) the speed of the lateral motors. This attitude motion is performed while keeping the total thrust constant. From this control principle, control input vector reads as [START_REF] Castillo | Real-Time stabilization and tracking of Four-Rotor Mini Rotorcraft[END_REF], [START_REF] Bouabdallah | Backstepping and sliding-mode techniques applied to an indoor micro quadrator[END_REF], [9] [START_REF] Hauser | Non linear control design for slightly non minimum phase systems: Application to v/stol aircraft[END_REF] with d the distance between the rotor shaft and the center of mass. Given the parameters d, b and C air0 , P is a full rank matrix. The gyroscopic torque induced by the four rotors reads as

    F T Γ 1 Γ 2 Γ 3     = P     ω 2 1 ω 2 2 ω 2 3 ω 2 4     P =     b 0 b 0 b 0 b 0 0 db 0 0 -db 0 db 0 0 -db 0 0 C air0 -C air0 C air0 -C air0     (
C b T gyros = - 4 i=1 J T ω i sk(Ω b )e b 3 ( 13
)
with J T the inertia matrix of a rotor and ω i the angular velocity of the rotor i. Consequently, F b and C b of system (2) may be written as [START_REF] Castillo | Real-Time stabilization and tracking of Four-Rotor Mini Rotorcraft[END_REF], [START_REF] Bouabdallah | Backstepping and sliding-mode techniques applied to an indoor micro quadrator[END_REF], [START_REF] Guenard | A practical Visual Servo Control for a Unmanned Aerial Vehicle[END_REF]:

F b = -F T e b 3 C b = Γ + C b T gyros (14)
Theoretically, translational and rotational dynamics (i.e.: ΣΓ = 0 3×1 , Σ = 0 3×3 ) are decoupled.

F b = -F T e b 3 -1 d sk(e b 3 )Γ C b = Γ (15)

IV. GENERIC 6DOF STATE MODEL

From the previous expressions of resulting forces and torques F b et C b , a generic nonlinear 6DOF state model is proposed covering a large class of reduced scale UAVs, i.e from aircraft to rotorcraft under quasi-stationary flight conditions.

A. Nonlinear model

Equation ( 2) describing the motion of a rigid body can be developped as follows

ξi = v i m vi = R F T e b T + mgR T e i 3 + F b air + ΣΓ + δF b Ṙ = Rsk(Ω b ) J Ωb = -sk(Ω b )JΩ + C b inter + C b air + C b T gyros + KΓ + δC b (

B. Zero dynamics

The study of PVTOL [START_REF] Hauser | Non linear control design for slightly non minimum phase systems: Application to v/stol aircraft[END_REF], [START_REF] Koo | Output tracking control design of a helicopter model based on approximate linearization[END_REF] brings into relief that the coupling matrix Σ between translational and rotational dynamics induce the Small Body Force ΣΓ from the torque control input Γ, which generates a zero-dynamics if Σ = 0 3×3 (Fig 5). In fact, assuming that the center of mass of the vehicle ξ i perfectly tracks the desired trajectory (ξ i ) d , one gets

ξ i -(ξ i ) d = 0 3×1 v i -(v i ) d = 0 3×1 vi -( vi ) d = 0 3×1 (17)
From ( 17) and the two first equations of ( 16), one gets

ΣΓ = -F T e b T -mgR T e i 3 -F b air -R T δF b (18)
By multiplying the last equation of ( 16) by Σ = ΣK -1 and substuting ΣΓ by equation ( 18), yields

ΣJ Ωb = Σ -sk(Ω b )JΩ b + C b inter + C b air +C b T gyros + δC b -F T e b T -mgR T e i 3 -F b air -δF b (19)
The attitude R d corresponding to the equilibrium position

Ω b = 0 3×1 , Ωb = 0 3×1 is given by mg(R d ) T e i 3 = -Σ C b air + C b inter + δC b +F T e b T + F b air + δF b (20)
R = R T R d denotes the error rotation matrix and may be linearized by

R = Id 3×3 + δαsk(∆e)
with ∆e ∈ IR 3 the rotation axis and δα ∈ [0 2π] the rotation angle. The third equation of ( 16) may be also linearized by ∆ = Ω b and ( 19) by ΣJ Ωb = -mgsk((R d ) T e i 3 )∆ with ∆ = δα∆e. Finally, one gets

S ∆ = -A ∆ ∆ = V ∆, A = mgU -1 ( ΣJ) T sk((R d ) T e i 3 )V -1 (21)
where the matrices U , S and V correspond to the singular value decomposition of the symmetrical matrix ( ΣJ) T ( ΣJ) (U and V are unity matrices and S is a diagonal matrix with nonnegative diagonal elements in decreasing order).

The matrix A depends exclusively of geometry and inertia of the UAV. If A is non Hurwitz, then the attitude diverges; else, the attitude periodically swings around the equilibrium attitude R d as illustrated for the helicopter [START_REF] Chriette | Zero dynamics analysis for ibvs control of under-actuated rigid body dynamics[END_REF]. This oscillation is not damped, undesirable but with a small amplitude, it may also be assumed to neglect Small Body Forces in the control design and to a posteriori verify its robustness [START_REF] Hauser | Non linear control design for slightly non minimum phase systems: Application to v/stol aircraft[END_REF], [START_REF] Koo | Output tracking control design of a helicopter model based on approximate linearization[END_REF].

C. Control strategy

F T ∈ IR and Γ ∈ IR 3 correspond respectively to the control inputs acting on the fuselage in order to control the position of the center of mass and the attitude of the rigid body. There is only one force to control the translational dynamics and three torques for the rotational dynamics. These vehicles may be also considered as underactuated.

• The guidance problem finally consists in the design of a control law for the translational dynamics

ξi = v i m vi = R(F T e b T + mgR T e i 3 + F b air + ΣΓ) + δF b (22)
F T Re T is considered as the control input so that ξ i converges to the desired position (ξ i ) d . A desired attitude R d is then deduced.

• The control problem consists in the design of a control law for the translational dynamics:

Ṙ = Rsk(Ω b ) J Ωb = -sk(Ω b )JΩ + C b inter + C b air + C T gyros + KΓ + δC b (23)
Γ is then the control input.

By neglecting the Small Body Forces (assumed stable), the state model ( 16) has a triangular structure adapted to the use of Backstepping techniques for instance. A first approach is to design a nonlinear controller with the full 4 th order state model by considering a dynamical extension of the thrust control input (i.e.: FT = u) [START_REF] Mahony | Robust trajectory tracking for a scale model autonomous helicopter[END_REF], [START_REF] Hamel | Adaptative estimation of aerodynamic forces for hover control of a reduced scale autonomous helicopter[END_REF], [START_REF] Chriette | Zero dynamics analysis for ibvs control of under-actuated rigid body dynamics[END_REF], [START_REF] Cheviron | Robust Control of a scale autonomous helicopter[END_REF]. An other way is to separate the full order state model into a slow and a fast timescale [START_REF] Kokotovic | A Control Engineer's Introduction to Singular Perturbation[END_REF], [START_REF] Krstic | Nonlinear and adaptative control design[END_REF]. Consequently, the translational dynamics represents the slow timescale and determines a desired attitude R d to reach which is viewed constant by the fast-timescale rotational dynamics. There is also no need of a dynamical extension of the thrust control input and the control design is simplified [START_REF] Mckli | Guidance and Control for Aerobatic Maneuvers of an Unmanned Airplane[END_REF], [START_REF] Pflimlin | A hierarchical control strategy for the autonomous navigation of a VTOL UAV[END_REF], [START_REF] Guenard | A practical Visual Servo Control for a Unmanned Aerial Vehicle[END_REF], [START_REF] Cheviron | Robust guidance and control of autonomous aerial vehicles in presence of wind gusts[END_REF] V. CONCLUSION

The main contribution of this paper was to present a generic nonlinear model of reduced scale UAVs in order to be simple sufficiently enough to design a controller. After a presentation of different architectures of some VTOL UAV (quadricopter, ducted fan and classical helicopter), a generic model focuses only on the key physical efforts acting on the dynamics is then proposed. In addition, the Small Body Forces expression, which can introduce a zero dynamics, was also studied.

1 ( 4 )

 14 The torque C b is composed by the gyroscopic torque C b P gyros induced by the propeller and the aerodynamic torque C b A C b = C b P gyros + C b A The gyroscopic torque induced by the propeller is expressed by C b P gyros = -J P ω P sk(Ω b )e bwith J P the inertia matrix of the propeller and ω P its angular velocity. As the aerodynamical force F b A , the aerodynamical torque C b A can be linearized which gives[START_REF] Mckli | Guidance and Control for Aerobatic Maneuvers of an Unmanned Airplane[END_REF]

Fig. 1 .

 1 Fig. 1. SIG Rascal 110.

  2): longitudinal and lateral cyclic angles control the main rotor flapping dynamics. When this dynamics tends to equilibrium, the total lift of the main rotor is tilted in comparison to the motor shaft[START_REF] Prouty | Helicopter, Performance, Stability, and Control[END_REF] a 1s b 1s -1] T

  ), one gets the expression of F b (Equation. 2). The resulting torque is composed by the propulsion momentum C b T = F Tm sk(l b Tm )e b Tm + F Tt sk(l b Tt )e b 2 induced by the lever arms l b T m(/t) = [l 1 T m(/t) l 2 T m(/t) l 3 T m(/t) ] T of the two rotors, the air resistance torque C b air = C airm e b 3 +C airt e b 2 , the damping torque C b damp = -κ Tm sk(e b Tm )e b 3 which resists the rotor disk distortion (κ Tm > 0 denote the elasticity constant of the main rotor) [26], the gyroscopic torque induced by the two rotors expressed by C b T gyros = -J Tm ω Tm sk(Ω b )e b 3 -J Tt ω Tt sk(Ω b )e b 2(10) with J Tm and J Tt , respectively, the inertia matrix of the main and the tail rotor, ω Tm and ω Tt , respectively, its angular velocity.By introducing Γ = [-F Tm a 1s F Tm b 1s F Tt ] T as control input vector, F b and C b may be written[START_REF] Chriette | Zero dynamics analysis for ibvs control of under-actuated rigid body dynamics[END_REF]

Fig. 2 .

 2 Fig. 2. VARIO Benzin Helicopter.

Fig. 3 .

 3 Fig. 3. Dragonflyer V TI-RC.

Fig. 4 .

 4 Fig. 4. HoverEye (Berthin Technology c ).

  16) with F T the resulting thrust, e b T the thrust direction, F b air the aerodynamic force, ΣΓ the Small Body Force induced by the torque control input Γ, δF b the disturbance force, C b inter the interaction torque between the different part of the vehicle, C b air the aerodynamic torque, C b T gyros the gyroscopic torque induced by the rotor, K the efficiency matrix of the torque control input Γ and δC b the disturbance torque. In appendix, the expression of each term on previous platforms is described.

Fig. 5 .

 5 Fig. 5. Scheme of UAV's dynamics.

4 C 1 L

 41 Aδa C Aδe C Aδr C Y δa C Y δe C Y δr C N δa C N δe C N δr sk(e 3 )

  The displacement of a streamlined surface in the air generates lift and drag forces. Assuming small angle of attack α and sideslip angle β, linearized aerodynamic forces F b A can be expressed as (with coefficients dimensioless C

* )

[START_REF] Mckli | Guidance and Control for Aerobatic Maneuvers of an Unmanned Airplane[END_REF]

Terms 'light' and 'reduced scale' implie that the flight range of aerial vehicles under interest is lower than 30mn.


E Perret
email: e.perret@cnr.tm.fr

M Lang

J Le Coz

A framework for detecting stage-discharge hysteresis due to flow unsteadiness : application to France's national hydrometry network

Keywords: rating curve, hysteresis, diagnosis, unsteady flow

A generic framework is proposed to evaluate the relative discharge error made when ignoring stage-discharge hysteresis due to transient flow over large gauging station networks. The diagnosis is conducted using the Jones equation, based on a simple hydraulic concept relating discharge to stage and its time-gradient. The main input data used for the method are the flow resistance coefficients, the temporal stage gradients, and the bed slopes. The hysteresis effect is quantified for each gauging station and mapped using the relative discharge error. The method was applied to 2618 gauging stations of France's national hydrometry network using observational data extracted from the national hydrological archive and from Digital Terrain Models. The diagnostic results highly depend on slope estimates used as inputs. Substantial hysteresis effects were found at stations with low bed slope combined with a fast flood regime. The France application shows the difficulty to provide a firm conclusion about stations prone to hysteresis due to the slope data uncertainty. This issue is not specific to France; slope estimates at a country-level is difficult to obtain in many countries. The use of local bed slope estimates is recommended to approach the slopes of the reaches controlling the station flow dynamics.

Introduction

Hysteresis in the stage-discharge relations (i.e. rating curves) of gauging stations can be observed during flood events [START_REF] Rantz | Measurement and computation of streamflow[END_REF][START_REF] Muste | Revisiting hysteresis of flow variables in monitoring unsteady streamflows[END_REF]. Figure 1 is an illustration of the effect of unsteady flows on different flow variables in time and on the stage-discharge relationship compared to the case of uniform steady flows. During the propagation of transient flows, the celerity of the pressure wave (stage) is smaller than the celerity of the velocity wave, hence smaller than the celerity of the discharge wave [START_REF] Graf | Flood hydrographs in open channels[END_REF]). In such case, for the same given stage, the discharge during the rising limb is higher than during the falling limb of the event, leading to a non-unique stage-discharge relation. Hysteresis creates a loop in the rating curve, which is more or less wide depending on the geometrical characteristics of the channel and on the type (intensity, gradient) of floods [START_REF] Lee | Evaluation of methodologies for continuous discharge monitoring in unsteady open-channel flows[END_REF]. Other phenomena can create looped rating curves, such as variable backwater (i.e. changes in the downstream conditions), or variable roughness [START_REF] Boyer | Streamflow measurement[END_REF][START_REF] Fenton | The calculation of streamflow from measurements of stage[END_REF][START_REF] Mansanarez | Non-unique stage-discharge relations: Bayesian analysis of complex rating curves and their uncertainties[END_REF]. In this study, the focus is on hysteresis induced by unsteady flow only.

In practice, this hysteresis effect is often neglected, partly because it is not captured by the occasional measurements of discharge and stage (a.k.a. gaugings).

Indeed, during flood events, gaugings are generally made after the flood peaks, for practical and safety reasons. Thus, the loop induced by the hysteresis phenomenon in the rating curve is often not observed. The evolution of the gauging techniques toward non-intrusive and less dangerous methods, e.g. radar [START_REF] Welber | Field assessment of noncontact stream gauging using portable surface velocity radars (SVR)[END_REF] and image velocimetry [START_REF] Dramais | Advantages of a mobile LSPIV method for measuring flood discharges and improving stagedischarge curves[END_REF], will certainly help to overcome this lack of information in the coming years. Substantial biases in flood prediction can arise if hysteresis is ignored, such as underestimation of discharge during the rising limb of the flood, including the peak discharge, time lag in the overall flood hydrograph and larger uncertainty of the discharge estimations due to the scatter of gauging data around the rating curve [START_REF] Holmes | River rating complexity[END_REF][START_REF] Mansanarez | Non-unique stage-discharge relations: Bayesian analysis of complex rating curves and their uncertainties[END_REF][START_REF] Muste | Revisiting hysteresis of flow variables in monitoring unsteady streamflows[END_REF].

Numerous methods exist to adjust the rating curves when unsteadiness is significant (cf. [START_REF] Lee | Evaluation of methodologies for continuous discharge monitoring in unsteady open-channel flows[END_REF]Dottori et al., 2009, for method reviews). They are usually based on a correction of the standard discharge, estimated with a unique rating curve valid for steady flow conditions, to account for hysteresis and compute the real discharge. [START_REF] Jones | A method of correcting river discharge for a changing stage[END_REF] method is the most used by the hydrometric community. It assumes that the flood wave propagates without any attenuation. Based on this kinetic wave approximation, the discharge in unsteady flows can be deduced from the stage and its time-gradient. More general expressions accounting for inertial forces have been proposed and used by [START_REF] Fread | Computation of stage-discharge relationship affected by unsteady flow[END_REF]; [START_REF] Fenton | The calculation of streamflow from measurements of stage[END_REF]; [START_REF] Perumal | Reproduction of hysteresis in rating curves[END_REF]; [START_REF] Petersen-Overleir | Modelling stage-discharge relationships affected by hysteresis using the Jones formula and nonlinear regression[END_REF]; [START_REF] Wolfs | Development of discharge-stage curves affected by hysteresis using time varying models, model trees and neural networks[END_REF]; [START_REF] Mansanarez | Non-unique stage-discharge relations: Bayesian analysis of complex rating curves and their uncertainties[END_REF]; [START_REF] Lee | Refinement of the fread method for improved tracking of stream discharges during unsteady flows[END_REF]; [START_REF] Muste | Revisiting hysteresis of flow variables in monitoring unsteady streamflows[END_REF] for example.

It is important to know where and when hysteresis can occur and how large it is to better evaluate the flood hazard [START_REF] Lee | Evaluation of methodologies for continuous discharge monitoring in unsteady open-channel flows[END_REF][START_REF] Muste | Revisiting hysteresis of flow variables in monitoring unsteady streamflows[END_REF]. The main objective of this paper is to propose a diagnostic approach to quantify the risk of hysteresis over a large set of hydrometric stations, typically over an entire national network.

The gauging stations prone to hysteresis are identified from the relative discharge error potentially made when ignoring the hysteresis effect. The parameter is defined as the relative bias between discharge estimates accounting (Q) or not (Q 0 ) for flow unsteadiness (Figure 1). It informs about the exposure of gauging stations to hysteresis effect during specific events. The parameter is similar to the PDIFF parameter introduced by Holmes (2016), which refers to the percent difference of the measured discharge from the discharge estimated using the unique standard rating curve. However, seems more suitable for massive diagnosis since it is based solely on discharge models and does not depend on available discharge measurements.

The diagnosis is intended to massive and large-scale deployment and to overview the areas influenced by hysteresis, before proceeding to a more accurate station-bystation analysis using available gaugings. The rest of the paper is organized as follows. First, the rating curve model for unsteady flows used in the diagnosis procedure is detailed in Section 2. Simple rating curve models such as the well-accepted [START_REF] Jones | A method of correcting river discharge for a changing stage[END_REF] equation are preferred for the diagnosis because it requires limited information that is available on large datasets and easily measurable on the field. The presented method is intended to be used for hydrometric purposes. Other more complex equations mentioned above require too much additional information and as shown by [START_REF] Mansanarez | Non-unique stage-discharge relations: Bayesian analysis of complex rating curves and their uncertainties[END_REF], they do not significantly improved nor reduced the uncertainty related to the estimation of the discharge during unsteady flows compared to [START_REF] Jones | A method of correcting river discharge for a changing stage[END_REF] equation. Then, the generic framework is described in Section 3 along with the required input data and the criteria retained for quantifying the hysteresis effect. The diagnosis is eventually applied to the gauging stations of the French national hydrometry network (Section 4). The limit of the method is pointed out using a sensitivity analysis of input data of the hysteresis model.

h Time Falling U,Q,h Q Rising Q max h max Unsteady flows Time U Q h Q 0,max ∆Q ∆Q Q 0 Q-Q 0 Q 0 = = ε s te a d y r a � n g c u r v e

Theory

Rating curve model for steady uniform flow

Hysteresis due to transient flows is known to be observed during intense floods.

In those conditions, it is usually relevant to approximate the hydraulic controls at gauging stations by a single channel control. The diagnosis is therefore performed with this assumption.

The Manning-Strickler equation is generally used for stage-discharge models of gauging stations with channel controls, for which the flow is mainly controlled by friction [START_REF] Rantz | Measurement and computation of streamflow[END_REF]World Meteorological Organization, 2010;[START_REF] Le Coz | Combining hydraulic knowledge and uncertain gaugings in the estimation of hydrometric rating curves: A Bayesian approach[END_REF] :

Q 0 ≈ KAR 2/3 h S 0 (1) with Q 0 [m 3 /s] the discharge in uniform steady flow, K [m 1/3 .s -1 ] the Strickler flow resistance coefficient, A [m 2 ] the wetted area, R h [m]
the hydraulic radius and S 0 [-] the channel slope. A unique relation between stage and discharge is thus obtained.

Equation 1 is valid only for uniform steady flow; although it is often used for nonuniform unsteady flows in the hydraulic community.

Rating curve model for unsteady flow -Jones equation

Equation 1 is based on the assumption that the energy slope S f can be approximated by the channel slope S 0 . It is not applicable in case of unsteady flow. Indeed, the water free-surface varies continuously during the flood wave propagation and the longitudinal water profile is not parallel to the river bed profile. Therefore, a correction of the steady-flow rating curve is required to capture the hysteresis effect. The energy slope S f can be expressed through the one-dimensional momentum equation of Saint-Venant, which describes the full dynamics of a flood wave propagation : 2) determines the type of wave occurring at specific sites and particular events [START_REF] Muste | Revisiting hysteresis of flow variables in monitoring unsteady streamflows[END_REF]. For example, neglecting the inertia terms compared to the pressure term and gravity force leads to the so-called diffusion wave (Equation 3):

S f = S 0 - ∂h ∂x - 1 g ∂U ∂t + U ∂U ∂x (
S f ≈ S 0 - ∂h ∂x (diffusion wave assumption) (3)
The diffusion wave assumption is generally accepted for low-gradient channels.

Combining Equation 1 and Equation 3, the stage-discharge model for unsteady flow becomes:

Q = Q 0 1 - 1 S 0 ∂h ∂x (4)
The discharge is therefore expressed as a steady-flow reference discharge Q 0 multiplied by a corrective term accounting for flow unsteadiness.

The [START_REF] Jones | A method of correcting river discharge for a changing stage[END_REF] approximation avoids estimating the longitudinal gradient term ∂h/∂x, which is rarely measured at gauging stations, except in the case of twin gauge stations [START_REF] Petersen-Øverleir | Bayesian analysis of stage-fall-discharge models for gauging stations affected by variable backwater[END_REF][START_REF] Mansanarez | Non-unique stage-discharge relations: Bayesian analysis of complex rating curves and their uncertainties[END_REF]. As a substitute, the temporal variation of stage ∂h/∂t is used, which is always available from stage records. The Jones approximation is based on the kinematic wave assumption, which assumes that the wave propagates with no attenuation along the channel. The flood wave celerity c can therefore be expressed as :

c = ∂x ∂t = ∂Q ∂A (5)
Assuming a prismatic channel and vertical river banks over the range of stage variation, the flood wave celerity becomes:

c ≈ 1 B ∂Q 0 ∂h (6)
where B [m] is the channel width.

The continuity equation of Saint Venant for quasi-steady flows can be rear-ranged to show the relationship between the longitudinal and temporal gradients (see [START_REF] Mansanarez | Non-unique stage-discharge relations: Bayesian analysis of complex rating curves and their uncertainties[END_REF] for details):

∂h ∂x = - 1 c ∂h ∂t (7)
With such approximation, the rating curve model can be expressed as follows,

where S f is expressed from the bed slope, the flood wave celerity and the stage time-gradient:

Q = Q 0 1 + 1 cS 0 ∂h ∂t (8)
Equation 8 is referred to as the Jones equation in the following.

In the rest of the paper, the channels at gauging stations are approximated by wide and rectangular channels with equivalent conveyance in order to perform the hysteresis diagnosis at large-scale with a minimum of information on the channel geometry. This assumption is acceptable in conditions of floods inducing hysteresis [START_REF] Le Coz | Combining hydraulic knowledge and uncertain gaugings in the estimation of hydrometric rating curves: A Bayesian approach[END_REF]. Then, the hydraulic radius can be approximated by the flow depth : R h ≈ (h -b), where b is the offset of the channel control. Equation 8becomes:

Q = KB(h -b) 5/3 S 0 1 + 1 5 3 K(h -b) 2/3 S 3/2 0 ∂h ∂t (9)

Quantification of the hysteresis effect

The hysteresis effect is quantified based on the relative discharge error between the discharge calculated considering the flow as unsteady (Q) and the discharge calculated assuming a steady flow regime (Q 0 ). Combining Equations 1 and 9, is expressed as follows :

= Q -Q 0 Q 0 = 1 + 1 5 3 K(h -b) 2/3 S 3/2 0 ∂h ∂t -1 (10)
Analysing Equation 10, the main requirements for detecting stations with hys-teresis due to flow unsteadiness are : 1) low S 0 , 2) high stage temporal gradients ∂h/∂t (i.e. stations subject to intense and rapid floods), and 3) rough beds inducing high flow resistances (i.e. low K).

Assuming that << 1 and after a first-order Taylor expansion, Equation 10 can be rearranged as follows :

≈ 0.3K -1 (h -b) -2/3 S -3/2 0 ∂h ∂t (
u 2 = u 2 K + 4 9 u 2 (h-b) + 9 4 u 2 S 0 + u 2 ∂h /∂t (12)
where u X = u X /X is the relative standard uncertainty of the variable X. The sensitivity coefficients in Equation 12 inform about the relative importance of each factor for the hysteresis effect quantification. The bed slope is the most sensitive parameter, with a sensitivity coefficient of 9/4. To a lesser extent, K and ∂h/∂t are also important; they have a sensitivy coefficient equal to 1. Equation 12indicates that (h -b) is the less sensitive factor, with a sensitivity coefficient of 4/9.

3 The diagnostic approach

Required data

The diagnosis consists in applying Equation 10 to each hydrometric station, which requires estimates of the flow resistance coefficient K, the bed slope S 0 , the bed elevation b, the stage time series and the time-gradient ∂h/∂t. rough gravel-beds, or beds with bedforms or presence of vegetation [START_REF] Coon | Estimation of roughness coefficients for natural stream channels with vegetated banks[END_REF]).

An accurate measure or evaluation of the bed slope S 0 is crucial for a good performance of the method (see Equation 12). It may be available at some gauging stations from hydraulic studies or specific topographic surveys. But over a national network, we need to use S 0 estimates from external datasets, e.g. based on Digital Terrain Models (DTM).

The time-gradient ∂h/∂t is estimated based on the stage time series h(t) that are automatically and continuously recorded at gauging stations (cf. Figure 2). Those records are then registered in large databases and converted into discharge time series Q(t) using rating curves established from discharge measurements (gaugings) at the stations. The higher gradients ∂h/∂t are assumed to be found during the most intense floods. The method scans the data to estimate the largest five time-gradients ∂h/∂t for each specific station. First, the largest five flood peaks are detected from the discharge time series Q(t) for each station over the entire period of available data. Samples of stage time series h(t) are extracted around those extreme events.

Data starting 10 days before the flood peak to 20 days after are kept. This duration was defined assuming it was sufficient to capture the flood dynamics at all stations.

The resulting stage data are smoothed using a spline function to remove noise and

Mapping hysteresis effect

Once the parameters K, b, h, S 0 and ∂h/∂t are identified (see Section 3.1), the discharge error can be computed using Equation 10. The procedure gives five estimates of for each stations, because the greatest five flood events are analysed (Figure 2).

The hysteresis effect is qualitatively assessed based on the resulting that we choose to present and map in two ways. Firstly, the focus is made on the most critical event per stations inducing the higher relative error max . The maximal discharge errors max are mapped by classes to facilitate the qualitative assessment.

Four classes are chosen according to usual criteria in hydrometry accounting for the uncertainties of the rating curves and gaugings : 1) negligible deviation between Q and Q 0 (i.e. max ≤ 1%), 2) low deviation (i.e. 1% < max ≤ 5%), 3) intermediate deviation (i.e. 5% < max ≤ 10%), and 4) high deviation (i.e. max > 10%), i.e.

larger than typical gauging uncertainties, which are generally evaluated to range 7-10% depending on the measurement technique. Secondly, the hysteresis effect can be assessed for each station by the number n representing the number of floods over the selected five for which the discharge error is larger than an arbitrary threshold, set as 10% in our application. The indicator n is divided into three classes to facilitate the hysteresis assessment using the map: 0 -station not prone to hysteresis, 1station with low exposure to this effect, from 2 to 5 -station prone to hysteresis effect.

The first map produces a general diagnosis and enables to detect all the potentially affected stations. The second map is complementary and gives a refined diagnosis showing the stations that are frequently affected by hysteresis. Those maps can be used as a guide for hydrological services to identify the affected stations and to better manage their stations. The thresholds for the assessment of hysteresis influence can be reviewed to match the requirements of each hydrological service.

4 Application to France's national hydrometry net- (2018), as it exhibits the best correlations to their validation dataset. The reach lengths vary from 156 m to 50 km with a mean length of ∼17 km in the world dataset. The smallest reach slope detected is 5.3×10 -5 ; smaller slopes are set to 0. In France, a specific correspondence is made between the reach and the location of the gauging station on the basis of the closest distance between both objects. The station slope is then set equal to the slope of the associated reach. If the distance exceeds 200 m between station and reach, the station is excluded. In France, the smallest slope is 2.5×10 -4 .

2. the theoretical hydrographical network (RHT) [START_REF] Pella | Le réseau hydrographique théorique (RHT) français et ses attributs environnementaux[END_REF]) (Figure 3b), which is a joined and oriented digital river network (see https://ecoflows .inrae.fr/software/). It is derived from the BD Alti ® 50 m digital elevation model of the French Geographic National Institute (IGN) and from the extended hydrological network (RHE) [START_REF] Pella | Réseau hydrographique naturel étendu (RHE) construit à partir de la BD Carthage[END_REF] <10 -5 10 -5 to 10 -4 10 -4 to 10 -3 10 -3 to 10 -2 10 -2 to 10 -1 >10 -1 The GloRS provides slope estimates only for the main rivers, whereas the two other datasets also include S 0 estimates for small rivers and headwater catchments.

The values of S 0 differ according to the chosen dataset, but the same general trend is observed, i.e. high and low slopes are detected approximately in the same areas (Figure 3). The slope estimates can vary for more than one magnitude order (Figure 4). Slopes from RHT and J2000 are in closer agreement than slopes from GloRS, which might be explained by differences in the DTM resolution.

Results over France

The diagnostic maps for France's national hydrometry network differ according to the chosen source of S 0 estimates (Figure 5). Slope estimates were available for the 2618 gauging stations using the RHT dataset whereas the GloRS and J2000-Rhône According to the diagnosis, 186 out of 2618 stations (i.e. 7.1%) are prone to substantial hysteresis effect ( >10%) based on S 0 estimates from the RHT (Figure 6a), whereas only 29 stations out of 1053 (i.e. 2.8%) are prone to hysteresis based on S 0 estimates from the GloRS (Figure 6b). Even if the diagnostic results are different, we can safely say that the hysteresis effect is low in France. Note that some lowgradient stations prone to hysteresis might have been missed because their slopes were lower than the threshold value imposed by the used datasets. In addition, further investigation is required to conclude if the detected stations are really prone to hysteresis, such as a specific analysis of gaugings and rating curves at individual stations.

Profiles of stations prone to hysteresis

The framework results can be used to determine the typical profile of stations prone to hysteresis. Figure 7 presents the needed geometrical and hydraulic conditions (i.e S 0 and ∂h/∂t) to observe hysteresis at a station. As expected, low bed slope combined with high temporal stage variation are the critical conditions. This is in accordance with the uncertainty analysis made for Equation 10. As shown in

Sensitivity to slope data

Testing diverse sources of bed slope data enables to evaluate the sensitivity of the diagnosis to those input data. The application over the Rhône basin using the J2000-Rhône data stresses the importance of using such accurate estimates of bed river slopes. Figure 8 shows the three diagnostic maps for the 99 stations shared by the three datasets in the Rhône basin. If the GloRS slope dataset is used as input, no stations are detected as prone to hysteresis; the relative discharge error being always lower than ∼8% (see Figure 9a). Conclusions are significantly different if RHT or J2000-Rhône slope datasets are used. Only 6 and 7 stations are affected by hysteresis according to the diagnosis made with S 0 estimates from RHT and J2000-Rhône datasets, respectively. Two of these stations are common to the two datasets.

Relative discharge error 0%<ϵ≤1% 1%<ϵ≤5% 5%<ϵ≤10% ϵ>10%

Figure 8: Diagnostic maps of hysteresis influence on the 99 gauging stations located in the Rhône basin, France, using slopes estimated from different DTM: GloRS, RHT and J2000-Rhône.

The resulting discharge errors highly depend on the S 0 estimates used for their computation. The discharge errors calculated using the GloRS slope dataset are mostly underestimated compared to those calculated using the RHT slope dataset (see regression line in Figure 9a). There is no obvious bias between deduced from RHT and J2000-Rhône slope datasets. A large scatter of values is nevertheless observed in Figure 9b. The high values of are realistic, in particular those calculated using J2000-Rhône slope estimates; in that case does not exceed 48% (Figure 9b).

Only one station out of the 99 stations in the Rhône basin has a discharge error that exceeds 100% ( =445% according to the diagnosis made with the RHT slope dataset). For clarity, this point is not presented in Figure 9. This station is also identified as prone to hysteresis by the diagnosis deduced from the J2000-Rhône slope dataset ( ≈15%). The out of range value is probably the consequence of a wrong S 0 estimate. Indeed, S 0 differs by more than one order of magnitude com-pared to S 0 from J2000-Rhône and GloRS datasets. In addition, the consistency of the other input data such as the time stage gradient was verified. Regarding the entire France network, there is only 15 out of 2618 stations (less than 0.6%) with a discharge error higher than 100% and all of these result from an estimation based on RHT slope dataset. This is encouraging with regard to the use of such a diagnostic approach.

Limit of the diagnostic approach

The main limit of the proposed framework is therefore related to the difficulty to estimate the bed slope massively and accurately. Slope calculation is scale-dependent and is thus sensitive to the spatial resolution of the DTM as well as the reach length.

The key to have a good diagnostic performance is to use local estimation of bed slope.

In the best case scenario, the slope should be evaluated over the length where the channel characteristics control the stage-discharge relation, but such a length is hard to identify and especially is site-specific. Such diagnosis is thus challenging, because we need to perform a diagnosis at large-scale (e.g. over country) of individual local objects (stations). No dataset at the local scale for the slope estimates exists, so we extract them from large-scale datasets that were not developed for local studies.

The splitting procedure for the reaches in the datasets is most often not suitable for our purpose, though accurate for hydrological studies. The reach length is often too long to capture the accurate value of the bed slope at the station. For example, the average reach length associated to the 1053 stations in the GloRS dataset is 27 km, whereas it is equal to 4 km and 6.1 km, for the 2618 and 298 stations in the RHT and J2000-Rhône datasets, respectively. Using such input data in the framework gives nonetheless a good overview of the global dynamics around the station in conditions of intense floods. In our framework, we mix data from different scales: averaged bed slopes over variable-length reach and local information at the station such as, the stage time series and geometrical characteristics. If there is any doubt about the value of at a specific station, we recommend measuring directly and locally the river bed slope around the station. The discharge error can then be re-calculated with good S 0 estimate. If the station is detected as prone to hysteresis effect, an adjusted stage-discharge relation should be used for discharge prediction rather than the standard rating curve, such as for example Equation 8, or variant models as those reported in [START_REF] Petersen-Overleir | Modelling stage-discharge relationships affected by hysteresis using the Jones formula and nonlinear regression[END_REF]; [START_REF] Mansanarez | Non-unique stage-discharge relations: Bayesian analysis of complex rating curves and their uncertainties[END_REF].

Conclusion

A diagnostic method is proposed to detect the gauging stations prone to hysteresis due to flow unsteadiness over a given hydrometry network. It uses well-known hydraulic concepts related to unsteady flows, such as the Jones equation, and simple data, i.e. easily measurable or estimable parameters on the field, such as those generally present in national hydrological databases. The final output is a map assessing the hysteresis effect at all the stations from the studied area through the relative discharge error between the discharge calculated with the steady flow rating curve and the actual discharge. The application of the method to France's national hydrometry network highlights the importance of using reliable bed slope estimates to produce a robust diagnosis. The major limit of the method is to have access to accurate local bed slopes. The diagnostic maps obtained using three different inputs of bed slope estimates show that in France the hysteresis effect due to flow unsteadiness is low for the majority of the gauging stations; 7.1% of the 2618 gauging stations are potentially affected by hysteresis according to the diagnosis made with the slope estimates from RHT dataset. However, the diagnostic results differ depending on the source (and accuracy) of the slope estimates. In the future, particular attention should be paid to measuring accurately the river bed slopes near the gauging stations and therefore refining the rating curves in case of unsteady flows with specific gauging campaigns.

Figure 1 :

 1 Figure 1: Discharge error due to flow unsteadiness hysteresis: temporal time series (U : average velocity, Q : discharge, h : stage), rating curves (dashed line for steady uniform flow regime and solid line for unsteady flow regime) and the relative discharge error comparing the two discharges of the two different regimes (Q: unsteady discharge, Q 0 : steady discharge).

  2) with g [m.s -2 ] the gravitational acceleration, U [m.s -1 ] the cross-sectional average water velocity, t [s] the time, x [m] the streamwise distance and h [m] the water surface elevation (a.k.a. the stage). The relative importance of the terms detailing the full dynamics of the flood wave propagation (Equation

  11) To investigate the relative importance of the different factors of Equation 11 in this hysteresis phenomena, i.e. K, S 0 , (h -b) and ∂h/∂t, the relative uncertainty of is assessed applying the uncertainty propagation law of JCGM (2008) (Joint Committee for Guides in Metrology) to Equation 11 and assuming independent measurement errors :

Figure 2 :

 2 Figure 2: The framework for diagnosing gauging stations influenced by hysteresis effect due to flow unsteadiness.

  keep only the stage variations reflecting the flood propagation. When the amount of missing data within the recorded stage time series is too large, it is excluded from the diagnosis. The gradient ∂h/∂t time series are then calculated from the smoothed h(t) and the five maximum values of the gradient over the five events are kept for computing .It is difficult to set the offset b of the channel control, mostly because it can evolve due to erosion/deposition, and in particular after floods inducing large sediment transport. As a consequence, b is seldom recorded in databases. If not available, we suggest setting b = h min in the framework, where h min corresponds to the lowest stage value of the stage time series of the studied event. Indeed, the cease-to-flow level b of the channel should not be higher than the lowest stage recorded.

work 4 . 1

 41 Presentation of the dataThe framework is applied to France's national hydrometry network. A number of 2618 gauging stations are analysed, which corresponds to stations with discharge rating (not stage only), including closed or discontinued stations. The stage and discharge times series required for the method to work come from the national hydrological archive (HYDRO2 database). Estimates for K, b, and S 0 are not available in HYDRO2 database. The flow resistance coefficient was assumed equal to 25 for all the stations as a conservative assumption. Approximates of channel bed slopes are deduced from three datasets using different DTMs for slope extraction with different resolutions and spatial coverages (Figure3) :1. the Global River Slope dataset, GloRS (Figure3a), which is a worldwide geospatial dataset detailed in[START_REF] Cohen | Global river slope: a new geospatial dataset and global-scale analysis[END_REF], where the slope is simply calculated from the elevation depression over the length of a river reach. The slope extraction is made based on the 15arc-sec resolution (∼ 460×460 m) SHuttle Elevation Derivatives at multiple Scales (HydroSHEDS) DTM and streamnetwork (see http://hydrosheds.cr.usgs.gov/index.php). The reaches are defined according to the confluence points and to an additional featuresplitting procedure, which splits the river segments that are longer than a user-defined distance. A 50-km splitting interval was selected inCohen et al.

  corresponding to an oriented simplification of the referential hydrographical network of the IGN, BD Carthage ® . The final RHT represents a network of 283 639 km with a total of 114 601 reaches associated with different topographic, climatic and hydrologic attributes. Various environmental attributes are available and presented in[START_REF] Pella | Le réseau hydrographique théorique (RHT) français et ses attributs environnementaux[END_REF], such as slope estimates, catchment areas, mean discharge, reach width and length for example. The slope is calculated as the difference in elevation between the upstream and downstream ends of the reach divided by the reach length. The lowest reach slope detected in the RHT is equal to 10 -4 ; smaller slopes are set to 0. The average reach length is 2.5 km and varies from a few meters to more than 40 km. The reach-station association is made by geographical proximity. On average, the distance between the two objects does not exceed 36 m.3. the slopes from the hydrological distributed model J2000-Rhône[START_REF] Branger | A distributed hydrological model to assess the impact of global change on water resources in the Rhône catchment[END_REF] (Figure3c), covering the Rhône River basin (97 800 km 2 ), located Southeast of France. The tool called HRU-delin (Hydrological Response Unit -delineation) combined with a DTM is used to prepare the mesh and generate several inputs for the hydrological modeling, including the creation of river reaches and calculation of their slopes (see https://forge.irstea.fr/ projects/hru-delin). The slope dataset results from the SRTM (Shuttle Radar Topographic Mission) digital elevation data, with a resolution of 90 m. The reach lengths are defined by the confluences and the locations of the gauging stations, hence a wide range of lengths, from 90 m to more than 44 km with a mean length of ∼6 km. The slopes are extracted for each reach by computing the difference in altitude between the upstream and downstream pixel of the reach. The minimum slope detected in the dataset is equal to 8.3×10 -5 ; smaller slopes are set to 0. A specific correspondence is made between the reach and the location of the gauging station to obtain the S 0 value at the station. The stations are first relocated on the stream network by proximity and similarity of drainage area. The reaches are then cut at the station position on the stream, so that the end of the reach corresponds to the station location.

Figure 3 :Figure 4 :

 34 Figure 3: River slope estimates in France from different datasets : a) GloRS, b) RHT and c) J2000-Rhône.

  datasets provide slope estimates for only 1053 and 298 gauging stations, respectively. The smallest French rivers are not represented in GloRS dataset, leaving some gauging stations without a slope estimation during the reach-station matching. The diagnostic maps based on the RHT and J2000-Rhône slope estimates reflect what was initially expected : a greater risk of hysteresis near estuaries (e.g. estuary of the Rhône River) and on large rivers with low-gradient slopes (e.g. Saône River); see Figure5d. In addition, stations located in mountains (e.g. Alps) are identified as not prone to hysteresis effect.

Figure 5 :

 5 Figure 5: Diagnostic maps of hysteresis influence on 2618, 1053 or 298 gauging stations in France, using a) RHT, b) GloRS or c) J2000-Rhône dataset as input data for river bed slope, respectively. d) Overview of the main rivers and mountains in France.

Figure 6 :Figure 7 :

 67 Figure 6: Relationship between classes of relative errors due to hysteresis effect and bed slopes S 0 deduced from a) the RHT and b) from the GloRS datasets. The red numbers indicate the number of stations within the class and the percentage with respect to the total number of stations, respectively.

Figure 9 :

 9 Figure 9: Comparison of discharge errors calculated using the slope estimates from different datasets for 99 gauging stations located in the Rhône basin: a) using GloRS versus RHT, and b) using J2000-Rhône versus RHT. Marginal boxplots detail the distribution of the data. Red solid and dotted lines represent the perfect agreement and the threshold of = 10%, respectively. Blue line refers to a linear regression line with its grey envelop corresponding to the confidence interval.

Acknowledgements

This work was funded by the SCHAPI, France's national hydrological service (SRNH 2018 and 2019 contracts about natural risks between INRAE and the Ministry of Ecology). The authors would like to gratefully thank Benjamin Renard (INRAE) for his contribution to the method development and for all the interesting discussions we had about the hysteresis subject. We also thank Flora Branger (INRAE) and Hervé Pella (INRAE) for their help for collecting the bed slope data from the J2000-Rhône model and the RHT, respectively.

Competing interests

The authors declare that they have no competing interests.


Guillaume Artigue
email: guillaume.artigue@mines-ales.fr

Anne Johannet

Valérie Borell

Séverin Pistre

Neural network flash flood forecasting : generalizing to ungauged basins

Neural networks for hydrology : adapting to available data

Statistical modeling on ungauged basins : how to ? 4. Results on ungauged basins of the Cévennes range

• Application of the model to 15 other "ungauged" basins • Bad results expected : taking advantage of the meaningfulness of errors ? • Building a correction method based on the basin characteristics.

Dimensional (left) and non dimensional (right) relationship between specific discharge and surface.

• For more than 30 sq.km basins and when representative rainfall are available, the model proved it was able to generalize to ungauged basins. • Its ability to generalize is limited to hydro-climatic similar areas but could be transposed to other areas.

• Acknowledgements to the French flood forecasting service (SCHAPI) for funding and data, to OHMCV for data and to Dominique Bertin from Geonosis for the development and constant improvement of the RNFPro Software.

•

  Numerous heavy rainfall episodes leading to flash floods, • More than 100 fatalities and billions of damage costs in the past years,• Low predictability and great difficulty to understand hydrological processes.

•

  Toukourou, M.S., Johannet, A., Dreyfus, G. and Ayral P.-A., 2011, Rainfall-runoff modeling of flash floods in the absence of rainfall forecasts: the case of "Cévenol flash floods", Journal of Applied Intelligence vol. 35, 2, pp. 1078-189. • Kong-A-Siou, L., 2011, Modélisation des crues de bassins karstiques par réseaux de neurones, cas de l'aquifère du Lez (Montpellier), thèse de l'Université Montpellier 2, 230p. • Artigue, G., Johannet, A., Borrell, V., Pistre, S., 2012 : Flash flood forecasting in poorly gauged basins using neural networks: case study of the Gardon de Mialet basin (southern France), Natural Hazards and Earth System Sciences 12, pp. 3307-3324, doi : 10.5194/nhess-12-3307-2012.