Giovanni De Nunzio
email: giovanni.de-nunzio@ifpen.fr

Laurent Thibault
email: laurent.thibault@ifpen.fr

Antonio Sciarretta
email: antonio.sciarretta@ifpen.fr

A Model-Based Eco-Routing Strategy for Electric Vehicles in Large Urban Networks

Keywords: Eco-routing, energy consumption estimation, electric vehicles, adjoint graph

A novel eco-routing navigation strategy and energy consumption modeling approach for electric vehicles are presented in this work. Speed fluctuations and road network infrastructure have a large impact on vehicular energy consumption. Neglecting these effects may lead to large errors in eco-routing navigation, which could trivially select the route with the lowest average speed. We propose an energy consumption model that considers both accelerations and impact of the road infrastructure, also separating the costs of all the possible turning movements in the transportation network by means of the adjoint graph. The approach makes use of only aggregated data available from existing map services. We demonstrate that the proposed strategy is more effective and reliable than the existing approaches in predicting vehicle energy consumption and in suggesting an energy-efficient route.

I. INTRODUCTION

While research on energy-efficient driving trajectories (i.e. eco-driving) has already witnessed numerous efforts [START_REF] Sciarretta | Optimal Ecodriving Control: Energy-Efficient Driving of Road Vehicles as an Optimal Control Problem[END_REF], the impact analysis of route choice on energy consumption and the development of energy-efficient navigations systems (i.e. eco-routing) are relatively new subjects. Current available navigation systems suggest either the fastest or the shortest route. However, these route choices are not always the most energy-efficient ones [START_REF] Ahn | The Effects of Route Choice Decisions on Vehicle Energy Consumption and Emissions[END_REF], [START_REF] Ericsson | Optimizing Route Choice for Lowest Fuel Consumption Potential Effects of a New Driver Support Tool[END_REF].

The main challenge of an eco-routing navigation system lies in the accurate estimation of the vehicle energy consumption, as demonstrated in [START_REF] Kubička | Performance of Current Eco-Routing Methods[END_REF], [START_REF] Guo | EcoMark 2.0: Empowering Eco-Routing with Vehicular Environmental Models and Actual Vehicle Fuel Consumption Data[END_REF]. Such energy consumption should be estimated on each of the different road segments composing the routing network of the area under analysis. Models for estimating vehicle energy consumption can be broadly divided into macroscopic and microscopic models, depending on how vehicular activities are aggregated over time and space. In a comprehensive study in [START_REF] Guo | EcoMark 2.0: Empowering Eco-Routing with Vehicular Environmental Models and Actual Vehicle Fuel Consumption Data[END_REF], a comparison of several existing energy consumption and emission models for combustion-engine vehicles was conducted. In general, microscopic models are considered more precise because they take as input instantaneous velocities and accelerations, which are difficult to obtain without equipped probe vehicles recording GPS or CAN-bus data. Eco-routing navigation systems based on this type of models need to acquire a large amount of driving data in order to determine a statistical cost on each road segment of the routing model. For instance, in [START_REF] Minett | Eco-Routing: Comparing the Fuel Consumption of Different Routes between an Origin and Destination Using Field Test Speed Profiles and Synthetic Speed Profiles[END_REF] the authors propose a strategy to reconstruct synthetic speed profiles from historical speed data to calculate an instantaneous power demand and an energy cost for the road segments. In [START_REF] Guo | EcoSky : Reducing Vehicular Environmental Impact Through Eco-Routing[END_REF], analogously, the authors propose to compute eco-weights based on GPS data, and to periodically maintain these weights as new information is collected. Intuitively, this approach requires a long time to have a reliable and dense mapping of the energy costs of all the road segments in the routing network.

On the other hand, macroscopic models take as input mean velocities, mean travel time and road grade, which are typically easier to obtain through free or commercial historical databases. Such aggregated models can be further classified into regression-based models and load-based models. The regression-based models are generally difficult to calibrate in order for all the parameters of the model to be significant, and they may lead to random errors, which are not easy to explain due to the little physical foundation of these models [START_REF] Ahn | Network-Wide Impacts of Eco-Routing Strategies : A Large-Scale Case Study[END_REF], [START_REF] Boriboonsomsin | Eco-Routing Navigation System Based on Multisource Historical and Real-Time Traffic Information[END_REF], [START_REF] Yao | Study on Eco-Route Planning Algorithm and Environmental Impact Assessment[END_REF]. The load-based models rely upon the longitudinal dynamics of the vehicle and are easier to calibrate using the vehicle construction parameters [START_REF] Richter | Comparison of Eco and Time Efficient Routing of ICEVs, BEVs and PHEVs in Inner City Traffic[END_REF], [START_REF] Jurik | Energy Optimal Real-Time Navigation System[END_REF].

All these models make strong assumptions for the energy consumption estimation, by neglecting accelerations, or road grade, or the impact of the road network topology and signalization. Only few recent works address the importance of considering accelerations even when using macroscopic models, as well as considering the elements of the road infrastructure that may cause vehicular speed disruptions. In [START_REF] Nie | An Eco-Routing Model Considering Microscopic Vehicle Operating Conditions[END_REF], the authors propose a method to consider interlink accelerations induced by the different adjacent average speeds. The costs for all turning movements available at a certain intersection are separated. However, the scalability issue of the routing problem is not addressed.

In addition, since electric vehicles and energy recovery phenomena are considered, the eco-routing strategy has to deal with the challenge of negative energy weights. Optimal path-searching algorithms on graphs with negative weights (e.g. Bellman-Ford algorithm [START_REF] Bellman | On a Routing Problem[END_REF]) are often discarded in literature because of their computational cost. Alternative approaches have been proposed to tackle this issue, mostly based on heuristics for graph weights shifting [START_REF] Jurik | Energy Optimal Real-Time Navigation System[END_REF] or nonoptimal path searching [START_REF] Abousleiman | Energy-Efficient Routing for Electric Vehicles using Metaheuristic Optimization Frameworks[END_REF].

In this work, we propose a novel macroscopic energy consumption model and a novel eco-routing strategy. The novelty of the proposed approach is summarized as follows. First, our macroscopic load-based energy consumption model considers accelerations to move from one road segment to the adjacent ones. Also, additional accelerations induced by the road topology and infrastructure (e.g. traffic lights, stop signs, turns, etc.) are considered in the energy consumption model. Then, since considering accelerations introduces a correlation between the energy cost of adjacent links, and the cost may not be unique due to the presence of multiple turning movements, the routing strategy is deployed on the weighted adjoint graph. This solves the size issues mentioned in [START_REF] Nie | An Eco-Routing Model Considering Microscopic Vehicle Operating Conditions[END_REF], because all the turning movements are decoupled avoiding the exponential growth of the routing graph. Finally, an optimal path-searching algorithm for graphs with negative weights is used, the computational time is showed to remain acceptable, and a strategy for a real-time use of the ecorouting navigation system is suggested.

The problem is presented and formulated in Section II. The simulation and experimental results are illustrated in Section III. Concluding remarks are given in Section IV.

II. PROBLEM FORMULATION

Eco-routing navigation aims to find the most energyefficient route in a vehicular road network to travel from an origin to a destination. It may be seen as an optimization problem where the function to be minimized is the energy consumption of the vehicle.

In the following we will introduce and detail the different elements of the proposed eco-routing strategy.

A. Vehicle Energy Consumption Model

Our analysis focuses mainly on the eco-routing for electric vehicles, therefore the energy consumption model is required to capture regenerative braking, as well as electric drive efficiency.

The vehicle longitudinal dynamical model may be generally written as [START_REF] Guzzella | Vehicle Propulsion Systems[END_REF]:

m dv(t) dt = F w -F aero -F friction -F slope ( 1
)
where m is the vehicle mass, v(t) the vehicle speed, F w the force at the wheels, F aero the aerodynamic force, F friction the rolling resistance force, F slope the gravity force. Therefore, the vehicle model shall be written as:

ẋ(t) = v(t) m v(t) = F w -1 2 ρ a A f c d v(t) 2 -mgc r -mg sin(α(x)) (2
) where ρ a is the external air density, A f is the vehicle frontal area, c d is the aerodynamic drag coefficient, c r is the rolling resistance coefficient, α(x) is the road slope as a function of the position, and g is the gravity.

The sum of the aerodynamic and rolling frictions, namely the road load force, is often approximated as a second order polynomial in the speed v:

F res = F aero + F friction = a 2 v(t) 2 + a 1 v(t) + a 0 (3)
where a 0 , a 1 and a 2 are parameters identified for a considered vehicle. Therefore, the force at the wheels can be expressed as:

F w = m v(t) + a 2 v(t) 2 + a 1 v(t) + a 0 + mg sin(α(x)) (4)
The torque requested from the electric motor to meet the force demand at the wheels is defined as:

T m =        F w r ρ t η t , if F w ≥ 0 F w rη t ρ t , if F w < 0 (5)
where r is the wheel radius, ρ t and η t are the transmission ratio and efficiency, respectively. The torque generated by the electric motor is saturated by a maximum torque T m,max and a minimum torque T m,min . In particular, during braking phases, if the motor torque is less negative than the negative saturation T m,min , then the vehicle is slowed down only by the regenerative brake. Otherwise, the mechanical brake works together with the regenerative brake.

The power available at the electric motor shaft, in the presence of regenerative braking mechanism, shall be written as follows:

P m =      T m,max • ω(t), if T m ≥ T m,max T m • ω(t), if T m,min < T m < T m,max T m,min • ω(t), if T m ≤ T m,min (6)
where ω(t) is the motor rotational regime and is defined as:

ω(t) = v(t)ρ t r (7)
Finally, the power demand at the battery of the electric vehicle, can be expressed as:

P b =    P m η b , if P m ≥ 0 P m η b , if P m < 0 (8)
where η b is the electric drive overall efficiency, here assumed to be the same both in motor and generator mode. Therefore, the battery energy consumption over the generic travel time T is obtained as:

E b = T 0 P b dt (9)

B. Eco-Routing Model

The road transportation network can be modeled as a graph. Let G = (V, A) be such a graph, where V is the set of vertices (or nodes) of cardinality n, and A is the set of feasible arcs (or links) of cardinality m connecting the nodes of the graph. Let us define a weighting function w : A → W , which associates each link of the graph with a weight.

In conventional routing graphs, the weight associated with each arc is either the length of the arc or its travel time. In the eco-routing framework, each link of the graph is assigned a weight that represents the travel energy expenditure.

The objective of this work is to design a reliable ecorouting strategy solely based on statistical and topological information about the road transportation network, without any use of real driving data. This kind of information, although often incomplete and/or imprecise, is generally available through commercial maps web-services. For each link i ∈ A of the graph, it is possible to know the link length l i , a prevailing average traffic speed vi which depends on the time of the day, and the road grade α i (x) which varies within the considered link depending on the position.

Note that a time-variant speed or acceleration profile is not available through the maps web-services, and the energy consumption model described in Section II-A cannot be directly used to assign the weights to each link of the graph.

1) Modified Energy Consumption Model: The model equations will be rewritten by replacing the time-variant speed v(t) with the average traffic speed v. All the vehicles on link i are supposed to travel at speed vi . Therefore, the force expression in ( 4) is modified for each link i as follows:

Fw,i = a 2 v2 i + a 1 vi + a 0 + mg sin(α i (x)) (10)
with no acceleration term. The motor torque becomes:

Tm,i =        Fm,i r ρ t η t , if Fw,i ≥ 0 Fm,i rη t ρ t , if Fw,i < 0 (11)
The electric motor rotational regime is also constant over time if constant speed is assumed:

ωi = vi ρ t r (12)
The mechanical power available at the electric motor is rewritten as:

Pm,i =      T m,max • ωi , if Tm,i ≥ T m,max Tm,i • ωi , if T m,min < Tm,i < T m,max T m,min • ωi , if Tm,i ≤ T m,min (13)
In the following, it is assumed that the saturation torque is independent from the motor regime. Finally, the power demand at the battery of the electric vehicle, considering the electric drive efficiency η b constant, can be rewritten as:

Pb,i =    Pm,i η b , if Pm,i ≥ 0 Pm,i η b , if Pm,i < 0 (14)
and ultimately, the battery energy consumption is:

Ēb,i = Pb,i T i ( 15
)
where T i = l i /v i is the travel time on link i when traveling at the average traffic speed vi .

2) Interface Accelerations: In order to take into account the acceleration effects and improve the model and the energy consumption estimation, the trip on each link is supposed to be composed of two phases: a cruising phase at the constant speed vi and an acceleration phase to go from vi-1 to vi .

The energy consumption associated with the change of velocity at a node between two links is defined as:

E jump,i = tjump,i 0 Pb,i dt (16)
where the battery power demand Pb,i is derived analogously from an interface force at the wheels defined as:

Fw = m • a + a 2 v(t) 2 + a 1 v(t) + a 0 (17)
The time-varying velocity v(t) in every transient is linearly modeled as v(t) = vi-1 + sign(v i -vi-1 ) • a • t, where vi-1 is the constant velocity on the incoming link, vi the constant velocity on the outgoing link, and a the constant acceleration to perform the jump. The speed change is performed in:

t jump,i = vi -vi-1 sign(v i -vi-1 ) • a (18)
Therefore, the total energy consumption on link i can be rewritten as:

E b,i = Ēb,i + E jump,i (19)
As mentioned before, a limitation of this modeling solution is that speed fluctuations within a link are neglected, which could be an issue in traffic. However, modeling the acceleration within a link seems unrealistic at a route planning level, because it highly depends on local traffic conditions, and as such is often spontaneous and unpredictable [START_REF] Nie | An Eco-Routing Model Considering Microscopic Vehicle Operating Conditions[END_REF].

3) Impact of the Infrastructure: The consideration of interface accelerations allows the energy consumption model to be more precise and realistic. However, as mentioned before, the available a-priori information about the average traffic speed is not always complete nor updated. In particular, it is very difficult and unlikely to have accurate speed information on secondary streets. This may lead to underestimations of the energy consumption and, therefore, lead the eco-routing algorithm to blindly favor this type of streets.

Disruptions in the speed profiles and accelerations are caused not only by traffic, but also by the infrastructure. In particular, critical elements of the road infrastructure, such as traffic lights, intersections, and turning movements are very likely to induce stops or significant decelerations.

Therefore, when the information about the position of these critical elements is available, the energy consumption term associated with acceleration in ( 16) should be modified.

For instance, the speed change between two links connected by a stop sign will be modeled as two distinct transients: the first transient from vi-1 to 0, the second transient from 0 to vi . Therefore, the energy consumption associated with the speed change at the node can be split in:

E jump,i = tjump1,i 0 Pb1,i dt + tjump2,i 0 Pb2,i dt (20)
where the speed variation in the first term is modeled as v 1 (t) = vi-1 -a • t, the time to perform the first transient is t jump1,i = vi-1 /a, the speed variation in the second term is modeled as v 2 (t) = a • t, and the time to perform the second transient is t jump2,i = vi /a. 4) Adjoint Graph: The consideration of the interface accelerations between adjacent links arises a major issue in the graph modeling of the road network. In particular, every node of the graph with two or more incoming links is critical because vi-1 and, consequently, E jump,i are not unique. Evidently, this prevents from assigning unique energy costs to the links of the graph. Hence, the graph G is not adequate for the proposed energy modeling approach. This difficulty is resolved by introducing the adjoint graph and by solving the routing problem on the adjoint graph instead of graph G.

The adjoint graph L(G) of a graph G has as nodes the links of G, and two nodes of L(G) are adjacent whenever the corresponding links of G are adjacent. The in-degree id(n j ) of a node n j ∈ V is the number of links entering n j , while its out-degree od(n j ) is the number of links leaving n j .

Let G = (V, A) be the graph with n nodes and m arcs as before. Then L(G) = (V * , A * ) is the adjoint graph with n * = m nodes, and

m * = n j=1 [id(n j ) • od(n j )] links.
Each link i ∈ A of the original graph G becomes a node i ∈ V * of the adjoint graph L(G). We can define as before an in-degree id(i) and an out-degree od(i) for node i ∈ V * .

Let us then define a weighting function for the links of the adjoint graph w * : A * → W * , where the weight for each link k ∈ A * is:

W * k = Ēb,i + E jump,i , if id(i -1) = 0 Ēb,i + E jump,i + Ēb,i-1 , if id(i -1) = 0 (21)
By using the adjoint graph L(G) as routing graph, the energy weights are uniquely assigned to each link of the graph by decoupling all the possible turning movements modeled in the original graph G.

Note that our objective in this framework is to address solely energy aspects, and therefore the weight assigned to each link of the adjoint graph represents only the associated energy consumption. The framework can be easily extended to consider also travel time, and the optimization would search then for a trade-off solution between energy consumption and travel time minimization.

C. Shortest Path Algorithm

Finally, the last stage of the problem formulation consists in the actual minimization of the energy consumption to travel from a selected origin to a destination in the road network. Minimization of a cost in a graph can be solved by means of a standard shortest path algorithm.

Such algorithm will be run on the adjoint graph. Note that the weights on the arcs can be negative, since regenerative braking is considered. Therefore, the Bellman-Ford algorithm is a necessary choice to find the global optimum [START_REF] Bellman | On a Routing Problem[END_REF]. The design and functioning scheme of the proposed eco-routing strategy is summarized in Fig. 1.

III. EXPERIMENTS

The simulation study and the experimental campaign were conducted on the road network of the city center of Lyon (i.e. about 5,400 nodes and 9,500 links), in France. The webservice HERE Maps [START_REF]HERE Maps[END_REF] was used in the experiments for obtaining information about network topology, road grade, current traffic conditions, and critical elements of the road infrastructure as previously discussed. The road grade and the road network topology, along with traffic lights position, are supposed not to vary with time. The traffic conditions, on the contrary, are extremely timedependent. An ideal real-time query on the web-service to retrieve the traffic conditions of the entire road network is not suitable for online use of the eco-routing navigation system.

A solution, more suitable for online purposes, may be to store offline the global information about historical traffic conditions for different days of the week and for different times of the day. The real-time adaptation is carried out only after the driver selects the desired origin, destination and departure time. The total cost of only the historical Kbest eco-routes is updated according to the current traffic conditions, and compared in order to determine the current best itinerary in terms of energy consumption.

The vehicle parameters are reported in Table I. The first set of simulation experiments aims to compare the eco-route suggested by our navigation system with the shortest and the fastest routes. In particular, the performance analysis has been conducted in terms of energy consumption computed by the proposed model and travel time. An energy consumption ratio (i.e. ratio between the energy consumption of the eco-route and the shortest or fastest route) and a travel interface accelerations and infrastructure critical points that induce stops and additional accelerations.

IV. CONCLUSIONS

A novel model-based eco-routing strategy for electric vehicles in large urban networks was proposed. The energy consumption model takes into account accelerations at the interface between adjacent road segments, as well as additional acceleration terms induced by critical road infrastructure elements. The introduction of the adjoint graph allowed for a correct and accurate energy modeling of the routing network. The proposed strategy exploits only information available online, without any model refinement by means of userrecorded experimental data. Future integration of travel time terms in the optimization framework, in order to discard impractical and unappealing eco-routes, is under study.

Fig. 1 :

 1 Fig. 1: Overview of the proposed eco-routing strategy.

Fig. 4 :

 4 Fig. 4: Test 1: driving speed profiles and mean traffic speed.

Fig. 5 :

 5 Fig. 5: Test 1: comparison between the real cumulative energy consumption and the energy consumption estimation.

TABLE I :

 I Vehicle parameters

	Description	Parameter	Value
	Vehicle mass	m	1190 kg
	Tire radius	r	0.2848 m
	Transmission ratio	ρt	5.763
	Transmission efficiency	ηt	0.95
	Acceleration	a	1.5 m/s 2
	Coefficient	a 0	125.73 N
	Coefficient	a 1	1.72 N/(m/s)
	Coefficient	a 2	0.58 N/(m/s) 2
	Motor minimum torque	T m,min	-50 Nm
	Motor maximum torque	Tm,max	200 Nm
	Electric drive efficiency	η b	0.85
	A. Simulation Study: Eco-Route Performance

ACKNOWLEDGMENTS

This project has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No. 653288 -OPTEMUS.

time ratio (i.e. ratio between the travel time of the fastest route and the eco or shortest route) are evaluated.

For a random selection of 1,000 different origindestination (O/D) pairs in the considered road network, the shortest-path algorithm has been run on the adjoint graph weighted with the energy consumption to compute the ecoroute, and on the original graph weighted with link length and travel time to compute the shortest-route and fastestroute. The Bellman-Ford algorithm convergence time was 0.02 seconds in average on a laptop computer equipped with a CPU at 2.8 GHz and 16 GB of RAM. For the sake of comparison, traffic conditions have been fixed at a specific date and time, namely at peak hour of a working day.

In Fig. 2 the probability distribution of the energy consumption ratio is reported. The energy consumption predicted by the model for the eco-route has been compared to the energy consumption of the shortest and fastest itinerary, respectively. The probability distribution of the energy gain provided by the eco-route with respect to the shortest route is represented in Fig. 2a. The mean value of the distribution is µ = 0.94, and the standard deviation is σ = 0.07. Therefore, the eco-route shows an average energy saving of about 6% with respect to the shortest route. The probability distribution of the energy gain with respect to the fastestroute is represented in the histogram in Fig. 2b. The mean value of the distribution is µ = 0.9, and the standard deviation is σ = 0.1. So, the eco-route shows an average energy saving of about 10% with respect to the fastest route.

An analogous analysis has been conducted on the travel time performance of the three itineraries calculated for each O/D pair, as reported in Fig. 3. The probability distribution of the travel time performance provided by the shortest-route with respect to the eco-route is represented in Fig. 3a. The mean value of the distribution is µ = 0.95, and the standard deviation is σ = 0.07. Therefore, the shortest-route shows an average travel time improvement of about 5% with respect to the eco-route. The probability distribution of the travel time performance provided by the fastest-route with respect to the eco-route is represented in Fig. 3b. The mean value of the distribution is µ = 0.93, and the standard deviation is σ = 0.07. Therefore, the fastest-route shows an average travel time improvement of about 7% with respect to the ecoroute. Note that the shortest-route can also present a worse performance than the eco-route in terms of travel time.

B. Experimental Campaign: Validation Results

A validation analysis was carried out in order to verify the accuracy of the proposed energy consumption model and, therefore, the reliability of the eco-routing navigation results.

In order to perform the validation analysis, several real driving profiles, consisting of GPS coordinates and driving speed, were recorded during an experimental campaign in the city center of Lyon. The driving profiles used for validation were all recorded at the same time of the day on different dates. The energy consumption associated with the real driving profiles was calculated using the instantaneous model described in Section II-A. The eco-routing model was constructed using the average traffic speed of a sample working day at the same time of the tests.

In the following, the validation results of two driving tests are reported. The first driving test was conducted on a combination of urban roads and highway. The second driving test was conducted on a secondary urban road. In the first test, as shown in Fig. 4, the real driving speed profiles in the dotted blue lines show good repeatability characteristics, though recorded on different days. This translates into a similar trip energy consumption, as shown in Fig. 5, where the cumulative energy expenditure is represented as a function of the traveled distance. The average value of the real energy consumption (i.e. dotted blue lines) at the end of the trip is 296.4 Wh, and the standard deviation is 12.5 Wh. The proposed energy consumption model, in green in Fig. 5, overestimates the total energy consumption by 7% with respect to the average of the real consumption. In order to further assess the performance and the advantages of the proposed strategy, the current state-of-the-art eco-routing approach based only on the average traffic speed was also tested on the validation test. As shown in magenta in Fig. 5, the energy consumption model based solely on the average traffic speed makes an estimation error of 32%.

Similar results can be observed on the second validation test, as shown in Fig. 6 and7. The average value of the real energy consumption is 246.3 Wh, and the standard deviation is 7.8 Wh. The proposed energy consumption model, in green in Fig. 7, underestimates the total energy consumption by 9%. On the contrary, the energy consumption model based solely on the average traffic speed makes a significantly larger error of 38%. Such model is not able to properly capture all the stops and accelerations of the real driving profiles in Fig. 6. The proposed model is more robust to imperfect information, thanks to the consideration of


Ian Ondo

Janice Burns

Christian Piedallu
email: christian.piedallu@agroparistech.fr

Including the lateral redistribution of soil moisture in a supra regional water balance model to better identify suitable areas for tree species

Keywords: Water balance, lateral flux, soil moisture, topography, digital map, species distribution model

To assess suitable areas for species, plant ecologists need accurate spatial information about available water for plants. Despite the recognized importance of topography in controlling soil moisture patterns, existing maps do not account for the redistribution of water through lateral fluxes. We included lateral fluxes in a GIS-based soil water balance model with the aim of evaluating the influence of lateral fluxes on soil moisture patterns and their importance to explain tree species distribution at regional scale. We used hydrological knowledge about lateral fluxes to map the distribution of monthly average soil moisture over the 1961-1990 period, for a 43,000-km² area in northeastern France. We then compared the ability of soil water estimated with or without lateral fluxes to explain the distribution of 19 common tree species. Spatial patterns significantly change when lateral fluxes are included in the model, with both large-scale effects due to variations in climate and soil properties, and local effects due to topography. The lateral redistribution given by the model revealed from 5% to 25% less water on the crests compared to in the valleys for metamorphic, sand and sedimentary bedrocks. Most of the tree species distributions studied were better explained when lateral fluxes were taken into account. Estimating soil moisture dynamics improves the ability to determine suitable areas for species at the landscape scale. It has major implications in the current climate change context owing to the potential to delineate topographic refugia or areas where species could colonize.

Highlights

• Contribution of lateral redistribution of water was included in soil moisture maps.

• Lateral runoff leads to 5 to 25% less water available for plants on the crests compared to the valleys. • Tree species distribution was better explained when taking lateral fluxes into account.

Introduction

Soil moisture is both recognized as one of the major determinants for plant composition and ecological processes, and one of the most difficult to estimate due to its high variability in space and time [START_REF] Porporato | Soil water balance and ecosystem response to climate change[END_REF]. The evaluation of the fine-scale spatial variability of available water for plants over large geographic areas and for long periods of time is crucial for plant ecologists in order to improve their understanding of species ecology and to adapt vegetation management to local conditions [START_REF] Barbour | North American terrestrial vegetation[END_REF][START_REF] Botkin | Environmental Science: Earth as a Living Planet[END_REF][START_REF] Chabot | Physiological ecology of North American plant communities[END_REF]. This knowledge is particularly important in the current climate change context, with an expected decrease in water availability in large parts of the world [START_REF] Bates | Climate Change and Water[END_REF]. The important spatial variability of soil moisture and the difficulties to obtain relevant datasets at the landscape scale make its estimation particularly difficult. It is often evaluated using the soil water balance (SWB, see table 1 for abbreviations), which estimates the amount of plant available water (PAW) for a defined period. Its calculation, based on the principle of the conservation of water contained in a volume of soil [START_REF] Choisnel | Le calcul du bilan hydrique du sol : options de modélisation et niveaux de complexité[END_REF], states that the amount of water entering is equal to the amount of water leaving, plus the change in the amount of water stored.

The variables involved are related to climate, soils and vegetation [START_REF] Dyck | Overview on the present status of the concepts of water balance models[END_REF][START_REF] Saxton | Soil water hydrology: simulation for water balance computations[END_REF]. Climate includes both precipitation (P) and potential evapotranspiration (PET), defined as the water demand of the atmosphere that would be possible under ideal conditions of moisture supply [START_REF] Thornthwaite | An approach toward a rational classification of climate[END_REF]. The soil-related components are linked to soil water holding capacity (SWHC), actual evapotranspiration (AET), and soil water runoff processes. The SWHC represents the maximum amount of water that plants can extract from the soil [START_REF] Granier | A lumped water balance model to evaluate duration and intensity of drought constraints in forest stands[END_REF], corresponding to the difference between the water contents at field capacity (Θfc) and the permanent wilting point (Θpwp). It depends on soil physical properties such as soil depth, texture, density and organic matter content, and the prospectable soil volume [START_REF] Bruand | Use of class pedotransfer functions based on texture and bulk density of clods to generate water retention curves[END_REF]. AET represents the amount of soil water delivered to the atmosphere both by evaporation and transpiration. Soil water runoff processes concern the surface runoff (water that flows at the ground surface), subsurface lateral fluxes (soil water that moves laterally) and percolation (soil water that flows downwards). Vegetation plays a significant role in the processes of interception and influences evapotranspiration through the transpiration of plants and the evaporation of soil and foliage [START_REF] Thornthwaite | An approach toward a rational classification of climate[END_REF].

Numerous water balance models have been developed at various time scales (e.g., hourly, daily, monthly and yearly) and to varying degrees of complexity [START_REF] Xu | A Review on Monthly Water Balance Models for Water Resources Investigations[END_REF]Singh, 1998, Schwärzel et al 2011). Most existing PAW maps available over broad areas have been comprised using simplified equations [START_REF] Van Der Schrier | Summer moisture variability across Europe[END_REF][START_REF] Zierl | A water balance model to simulate drought in forested ecosystems and its application to the entire forested area in Switzerland[END_REF]. They are often based on models at the monthly scale pioneered in the middle of the last century by Thornthwaite and Palmer [START_REF] Palmer | Meteorological drought[END_REF]Thornthwaite and Mather, 1955). They do not take hydrological fluxes into account, despite their importance in influencing soil moisture patterns at the toposcale. Indeed, many hydrological studies showed a redistribution of the soil moisture gradient along the hillslope gradient [START_REF] Brocca | Soil moisture spatial variability in experimental areas of central Italy[END_REF][START_REF] Ticehurst | Using a physically based model to conduct a sensitivity analysis of subsurface lateral flow in south-east Australia[END_REF], with large variations depending on the season and precipitation [START_REF] Weyman | Measurements of the downslope flow of water in a soil[END_REF]. The effect of topographical position has also been observed on vegetation. Several studies attributed changes in species composition [START_REF] Deblauwe | The global biogeography of semi-arid periodic vegetation patterns[END_REF][START_REF] Johnson | Ebrary Academic Complete Subscription[END_REF] and productivity [START_REF] Berges | Sessile oak (Quercus petraea) site index variations in relation to climate, topography and soil in even-aged high-forest stands in northern France[END_REF][START_REF] Curt | Alimentation en eau et production forestière, application d'indicateurs simples pour les résineux dans le massif central[END_REF][START_REF] Kobal | Influence of soil properties on silver fir (Abies alba Mill.) growth in the Dinaric Mountains[END_REF] along the topographical gradient to variations in soil moisture, suggesting that lateral fluxes could be an important consideration in the study of plant ecology. Lateral fluxes can be estimated by hydrological models such as TOPMODEL [START_REF] Beven | A physically based, variable contributing area model of basin hydrology[END_REF], TOPOG [START_REF] Oloughlin | Saturation regions in catchments and their relations to soil and topographic properties[END_REF], WET [START_REF] Moore | GIS and landsurface-subsurface modeling[END_REF] and SMR [START_REF] Frankenberger | A GIS-based variable source area hydrology model[END_REF], using soil hydraulic conductivity at saturation (Ksat: maximum rate at which a soil can transmit water) and the shape of the surface topography as data inputs. Most of them aim to determine where the runoff takes place in the catchment to reproduce river discharge at the basin outlet [START_REF] Beven | Spatially distributed modeling: conceptual approach to runoff prediction[END_REF][START_REF] Xu | A Review on Monthly Water Balance Models for Water Resources Investigations[END_REF]. They are not suitable to study plant distribution over large areas for different reasons:

-they do not describe PAW spatial variation. Moreover, many of them are semi-distributed, which means that hydrologically similar portions of the watershed are lumped together and are characterized by averaged ecological conditions, which do not provide precise estimation of soil moisture; -they are based on fine time step calculations, estimating PAW for long periods of time and over broad areas can be too time-consuming and data is not always available;

-some parameters should be calibrated at the catchment scale and cannot be extrapolated.

The aim of this study was to include lateral fluxes in a simple GIS-based soil water balance model at regional scale, to improve the estimation of available water for trees and evaluate the importance of lateral fluxes to explain tree species distribution. By building a program that could easily use available input data and that does not need calibration, we estimated the monthly time step average soil moisture for the 1961-1990 reference period, accounting or not for the lateral fluxes, in a 43,000-km² area in northeastern France. We used the model outputs to evaluate the influence of lateral fluxes in SWB, and determined their ability to explain the distribution of the most common tree species present in the study area.

Materials and Methods

To estimate available water for plants, we developed a fully-distributed water balance model coupled with a Geographic Information System (GIS) that requires easily available variables. To account for runoff that usually occurs at daily or shorter time scales, the model uses a daily subroutine whose soil water balance components are aggregated on a monthly basis to provide average monthly values of PAW that are representative of a long period of time in order to be related to tree species distribution. The routine was implemented and launched from R statistics software and executed in the environment of GRASS GIS through the R interface library for GRASS 6.4 spgrass6.

Input data

Climatic data

The study area (43,000 km²) covers a large climatic gradient in northeastern France, with altitudes ranging between 140 and 1424 m (Fig. 1), and mean annual temperature and precipitation ranging between 6 to 10.5° and 400 to 2400 mm, respectively. Average PAW values are required over long periods of time to understand plant distribution patterns, whereas runoff estimation requires data on a daily or shorter time scale. Since accounting for runoff at a daily time step for many decades is too time-consuming, we used P and PET daily values for a year that were representative of the 1961-1990 period.

Because averaging P over this time period will result in an unrealistic sequence, we disaggregated monthly average P for the 1961-1990 period into the most likely sequence of daily rainfall events (Supporting information S1). Daily PET were obtained using the Turc formula [START_REF] Turc | Estimation of irrigation water requirements, potential evapotranspiration: a simple climatic formula evolved up to date[END_REF]. This requires solar radiation obtained by dividing monthly values by the number of days in the month, and temperature obtained by averaging daily values over the 1961-1990 period for each Julian day of the year.

Fig.1: Elevation of the study area (m).

Mean temperatures were extracted from the SAFRAN model [START_REF] Quintana-Seguí | Analysis of near-surface atmospheric variables: Validation of the SAFRAN analysis over France[END_REF][START_REF] Vidal | A 50-year high-resolution atmospheric reanalysis over France with the Safran system[END_REF] and solar radiation from the Helios model [START_REF] Piedallu | Multiscale computation of solar radiation for predictive vegetation modelling[END_REF]. For P, we extracted averaged monthly values over the period 1961-1990 from a 1km resolution map [START_REF] Bertrand | Changes in plant community composition lag behind climate warming in lowland forests[END_REF] and occurrences of daily rainfall events used for disaggregation from SAFRAN data [START_REF] Quintana-Seguí | Analysis of near-surface atmospheric variables: Validation of the SAFRAN analysis over France[END_REF].

Soil properties

The model requires knowledge about the soil thickness (D), SWHC, hydraulic conductivity (K sat ) and volumetric water content at saturation (Θ sat ), field capacity (Θ fc ) and wilting point (Θ pwp ). All these data were gathered from plots monitored by the French National Forestry Inventory (IFN), part of the French National Geographic Institute (IGN). We used 15,369 plots distributed according to a quasi-systematic sampling over the forests, with a mean distance of approximately 1 km [START_REF] Drapier | La base de données écologiques de l'IFN[END_REF]. For each plot, a maximum of two horizons were determined as well as the texture class to which they belong, at a soil pit scale of a maximum depth of 1 m. For each texture horizon, Θ sat , Θ fc and Θ pwp were estimated with Al Majou et al. (2008b) class pedotransfer function (class-PTF) and K sat was estimated using Schaap et al. (2001) class PTFs (see Supporting information S6). SWHC was estimated according to Al Majou et al. (2008a) as the difference in potential water content between -100 hPa and -15000 hPa, and taking both rock outcrops at the plot scale and the stone content for each horizon into account [START_REF] Piedallu | Mapping soil water holding capacity over large areas to predict potential production of forest stands[END_REF]. D was estimated by summing the thickness of each horizon. For each plot, average values of hydraulic parameters for each horizon were weighted according to the thickness of soil horizons. They were assigned to the whole soil profile and mapped at a 50-m resolution. SWHC, Θ sat and Θ fc were mapped using linear models as a function of geology, altitude and topographic features derived from a digital elevation map [START_REF] Mcbratney | On digital soil mapping[END_REF]. In order to ensure consistency between databases, the amount of water at permanent wilting point was obtained by the difference between the amount of water at field capacity and the amount of water at SWHC. K sat maps were obtained through interpolation using ordinary kriging, which gave better results than linear models.

Modeling approach

Soil water balance calculation

Our calculations are based on Thornthwaite's model (Thornthwaite and Mather, 1955;[START_REF] Thornthwaite | Instructions and Tables for Computing Potential Evapotranspiration and the water balance[END_REF]) that provides monthly estimations of PAW [START_REF] Piedallu | Soil water balance performs better than climatic water variables in tree species distribution modelling[END_REF]

AET t = PAW t-1 + (P t -PAW t ) (2b) DE t = PET t -AET t (2c)
We simulated two reservoirs of water storage whose maximum capacities are given by the difference between Θ sat and Θ fc , respectively, corresponding the soil water drainage capacity (SWDC), and the difference between Θ fc and Θ pwp , corresponding to the SWHC. SWDC controls the water surplus, i.e., the runoff available soil water (RAW), whereas SWHC controls the PAW (Fig. 2).

Fig.2:

The double reservoir for keeping track of incoming and outgoing water fluxes in the soil (see Table 1 for abbreviations).

P falling from the atmosphere directly reaches the ground surface. The water can then either completely infiltrate the soil or be divided into two parts depending on the amount of P, initial soil moisture, land cover, slope and soil type that cannot infiltrate the soil generates water at the ground surface available for surface runoff. The part of P that infiltrates the soil can either be stored in the soil water column in the case where SWC is not recharged up to lateral subsurface runoff, corresponding to At the end of the time period t, the soil moisture storage of a given grid cell is then computed as follows:

SWC t = where Q surf = surface runoff and incoming/outgoing water fluxes. SWC is set to SWHC at the beginning of the first month.

PAW t is calculated using Formula 3, and is equal to SWC The daily subroutine first calculates the soil water balance at the daily scale in a loop over one month. At the end of the month, cumulative daily values of P, AET, Q Q sub are used to compute monthly values of SWC and PAW using Equations (3) monthly routine then uses these monthly SWC and PAW as starting values for the next month and so on.

The double reservoir for keeping track of incoming and outgoing water fluxes in the soil (see Table 1 for abbreviations).

P falling from the atmosphere directly reaches the ground surface. The water can then either completely infiltrate the soil or be divided into two parts depending on the amount of P, initial soil moisture, land cover, slope and soil type (Fig. 3). that cannot infiltrate the soil generates water at the ground surface available for surface runoff. The part of P that infiltrates the soil can either be stored in the soil water column in the case where SWC is not recharged up to field capacity (FC), or move laterally by , corresponding to interflow (Fig. 3).

At the end of the time period t, the soil moisture storage of a given grid cell is then

= SWC t-1 + P t -AET t ± Q surf t ± Q sub t
surface runoff and Q sub = lateral subsurface runoff. The sign ± stands for incoming/outgoing water fluxes. SWC is set to SWHC at the beginning of the first is calculated using Formula 3, and is equal to SWC t within the limit of The daily subroutine first calculates the soil water balance at the daily scale in a loop over one month. At the end of the month, cumulative daily values of P, AET, Q are used to compute monthly values of SWC and PAW using Equations (3) monthly routine then uses these monthly SWC and PAW as starting values for the next 8

The double reservoir for keeping track of incoming and outgoing water fluxes in the P falling from the atmosphere directly reaches the ground surface. The water can then either completely infiltrate the soil or be divided into two parts depending on the (Fig. 3). The part of P that cannot infiltrate the soil generates water at the ground surface available for surface runoff. The part of P that infiltrates the soil can either be stored in the soil water column or move laterally by At the end of the time period t, the soil moisture storage of a given grid cell is then

(3) = lateral subsurface runoff. The sign ± stands for incoming/outgoing water fluxes. SWC is set to SWHC at the beginning of the first within the limit of FC.

The daily subroutine first calculates the soil water balance at the daily scale in a loop over one month. At the end of the month, cumulative daily values of P, AET, Q surf and are used to compute monthly values of SWC and PAW using Equations (3). The monthly routine then uses these monthly SWC and PAW as starting values for the next

Fig.3: Schematic pathways of water movements simulated in the model (see Table 1 for abbreviations).

Surface runoff

Q surf (in mm) is the available surface runoff generated at the daily time scale by P in excess and computed using a modification of the Soil Conservation Service Curve Number method (SCS-CN, USDA, 1972;[START_REF] Usda | Urban hydrology for small watersheds[END_REF]:

Q surf = (P -I a )²/[(P -I a )+S] if P > I a (4) Q surf = 0 if P ≤ I a where:
-S is the potential retention in mm; -I a is the initial abstraction in mm (Supporting information S2).

The retention parameter S is related to soil types, land slope, antecedent soil moisture and land cover via the curve number CN according to the following relationship (USDA, 1986): S = 25.4 (1000/CN -10) (5) where CN ranges from 0 under totally dry conditions, to 100 under water-ponded conditions (Supporting information S2).

Once Q surf is calculated, it is discharged downslope using Manning's formula:

V = 1/n R 2/3 Slope 1/2 (6a) Q out,surf = V B Q surf (6b) R i = a (A i ) b (6c)
where: -V is flow velocity (m.s -1 ), limited between 0.02 m.s -1 to 2 m.s -1 , corresponding to a speed range that occurs naturally across a hillslope (Grimaldi et al., 2010;[START_REF] Maidment | Unit hydrograph derived from a spatially distributed velocity field[END_REF], -Q out,surf is the surface runoff discharge (m 3 .s -1 ), -n is the Manning roughness coefficient (s.m -1/3 ) (Supporting information S7), -R i is the hydraulic radius at cell i (m), -Slope is the slope at ground surface (m.m -1 ), -B is the flow width (cell width) (m), -A i is the upstream drainage area at cell i (m²), -a is a network constant (2.4 x 10 -4 ) and b a geometric scaling exponent (0.5), empirically estimated with respect with flow velocity boundaries imposed.

The quantity of surface water moving laterally as a sheet flow is limited by Q surf , and is equal to:

Q surf (mm) = minimum value between Q out,surf and Q surf (7)

Lateral subsurface runoff

The lateral subsurface runoff is calculated from a simple relationship for saturated flow based on Darcy's law [START_REF] Darcy | Les fontaines publiques de la ville de Dijon[END_REF] and a kinematic approximation of hydraulic gradient i.e., the hydraulic gradient is equal to the land slope at the site [START_REF] Beven | kinematic subsurface stormflow[END_REF][START_REF] Frankenberger | A GIS-based variable source area hydrology model[END_REF]. Under saturated flow conditions, the pores of the soil layer are full of water and flow discharge depends on hydraulic conductivity at saturation and the contributing saturated depth of the soil layer:

Q out,sub = K sat D s w tan β (10 -3 /Area of cell) (8)
where:

-Q out,sub is the subsurface lateral discharge (mm), -K sat is the (horizontal) hydraulic conductivity at saturation (m.d -1 ), -D s is the saturated depth area (m), -w is the flow width (dimension of the cell in m), -tan β is the tangent of the land slope (m/m or unitless) and 10 -3 /(Area of cell) is a factor of conversion to convert discharge from m 3 .d -1 to mm.

As the soil layer dries, a saturated area is assumed to build up at its bottom because of the tendency of water to accumulate downwards. Therefore, D s is assumed to be a function of the total soil layer depth D and relative soil water content storage:

D s = D (Θ -Θ fc / Θ sat -Θ fc ), Θ > Θ fc , (9) D s = 0, Θ ≤ Θ fc , where:
-D s is the saturated depth of soil layer, -Θ is the actual soil water content (cm 3 /cm 3 ),

-the term (Θ -Θ fc /Θ sat -Θ fc ) determines the fraction of the soil layer contributing to saturated flow.

The quantity of subsurface water draining laterally is limited by RAW and is equal to: Q out,sub = minimum value between Q out,sub and RAW (10)

Lateral redistribution of fluxes

The lateral redistribution of surface flow and lateral subsurface flow is processed by an automatic drainage path algorithm that divides the discharge of a particular cell among its downhill neighbors, based on the maximum downslope gradient (MFD-md) proposed by [START_REF] Qin | An adaptive approach to selecting a flow-partition exponent for a multiple-flow-direction algorithm[END_REF]:

d i = tanβ i f(e) L i / ∑ (tanβ j f(e) L j ) ( 11
)
Each cell is potentially surrounded by eight cells, and MFD-md allocates a fraction d i of the discharge from a particular cell into its i th neighboring cell; tanβ i is the directional slope with the i th neighbor (the difference in elevation divided by the horizontal distance between the centers of the cells); L i is the ''effective contour length'' of cell i (0.5 and 0.354 for downslope cells in cardinal directions and diagonal directions, respectively, and 0 for non-downslope [START_REF] Quinn | The prediction of hillslope flow paths for distributed hydrological modeling using digital terrain models[END_REF].

The function f(e) is used for modeling the flow [START_REF] Qin | An adaptive approach to selecting a flow-partition exponent for a multiple-flow-direction algorithm[END_REF]: f(e) = 8.9 min(e,1) + 1.1 ( 12) where e refers to the maximum downslope gradient.

The lateral inflow that a cell receives from its surrounding upslope neighbors is defined as:

Q surf in = ∑ d i Q surf out,i Or, Q sub in = ∑ d i Q sub,i
(13) where Q surf in /Q sub in is the quantity of runoff that enters the given cell and is equal to the cumulative contributions of water flow from neighbors that drain it. According to our flow velocity calculations, we allowed the subsurface water to flow from one cell to another in one day. In agreement with the USDA (1972), surface flow was restricted to a maximum of six cells per day.

Alluvial zone delineation

In the model, we accounted for fluvial processes that operate in the alluvial zone, subjected to variations in the water table. We delineated areas where Θ is always considered greater or equal to Θ fc using the depth-to-water (DTW, in m.) index at site i [START_REF] Murphy | Topographic modelling of soil moisture conditions: a comparison and verification of two models[END_REF]:

DTW i = [Σ (dz/dx)a]w c (14)
where [Σ (dz/dx)a] is the cumulative slope (sum of slope values) along the least cost path connecting any point of the landscape to a watercourse (m); a is a multiplier that is equal to 1 when the path is in the cardinal direction, and 1.414214 when it is diagonal; w c is the grid cell size (m).

To delineate the extent of the alluvial zone on either side of the watercourse, we assumed that:

-If DTW ≤ 1.6, then Θ ≥ Θ fc for main streams.

-If DTW ≤ 1, then Θ ≥ Θ fc for secondary watercourses.

-Every single grid cell located in a water body has an actual water content fixed to saturation i.e., Θ = Θ sat for watercourses, lakes, etc.

Model run and evaluation

The model was run over one complete year from January to December, assuming that SWC is at saturation for every single grid cell at the first day of the year (which also means that PAW is equal to SWHC).

Comparison between soil moisture estimated with and without lateral fluxes

In order to evaluate the importance of fluxes in SWB calculation, the model was run both with and without fluxes (a suffix f indicates that lateral fluxes are included). We calculated monthly maps for volumetric SWC, PAW and DE for the 1961-1990 period and for the growing season (March-September). The differences between SWC and SWC f indices were studied according to four geological strata (calcareous, metamorphic, sandstone and sedimentary; see Fig. 4), and three topographic strata (crest, mid-slope and valley; see Fig. 4), using a uniform mesh of sample points spaced on a 1-km grid over the study area (n=27311 plots when removing flat areas). The geological strata were identified through the aggregation of detailed soil information from the 1:50000 scale map provided by the French Geological Survey (BRGM), to determine areas of relatively homogenous soil textures that can lead to similar hydraulic behavior. The topographic strata allowed comparison for similar slope position and were calculated using Jenness Enterprises' Land Facet Corridor Designer (Jenness et al., 2013) and a 50-m digital terrain model input (Supporting information S3). The median, first quartile and third quartile values of the SWC f and SWC indices were plotted monthly and analyzed for each of the 12 strata using bootstrapping, with 100 random samples (n=350) with replacement.

Ability of lateral fluxes to improve tree species distribution models

Because validating such soil moisture maps would require a dataset of measurements collected over many decades at a regional scale, which is currently impossible to obtain, we have chosen to compare the ability of SWB and SWB f to explain tree species distribution patterns. A dataset of 2,419 georeferenced plots with the presence or absence of 19 common tree species was extracted from the French National Forest Inventory (IFN) database collected from 2005 to 2014 [START_REF] Drapier | La base de données écologiques de l'IFN[END_REF]. All the tree species that were present in the study area were studied except those that are commonly planted and those that had less than 50 presences in the dataset. Their distribution was modeled using generalized additive models (GAMs) with a logistic link function, binomial error distribution and smoothing splines (df=4) [START_REF] Hastie | Generalized additive models[END_REF]. Mean annual temperature for the 1961-1990 period [START_REF] Bertrand | Changes in plant community composition lag behind climate warming in lowland forests[END_REF] and soil pH were systematically included as explanatory variables and were supplemented by the water variables. Four combinations of water variables calculated for the growing season were evaluated per species: (i) including PAW or DE, respectively; and (ii) PAW plus the difference between PAW f and PAW (Diff PAW-PAWf ) or DE plus the difference between DE f and DE (Diff DE-DEf ). Each model was evaluated using the explained deviance (D²). The difference in D² between models including or not including Diff PAW-PAWf and Diff DE- DEf was evaluated, and its significance was estimated using an ANOVA with a Chi-square test.

Fig.4: Delineation of the four geological strata (on the left) and the three topographical strata (on the right) 3 Results

Comparison between SWB and SWB

SWB and SWB f indices representative for the 1961 monthly and mapped at a 50 calcareous soils in winter, to 25% for sandstone in PAW f shows considerable spatial variations both at the global and local scale, with lower values in mountains and in calcareous plateaus and near steep areas (Fig. 6). Lateral redistribution along the topography content near the crests than differences occur between SWC and SWC P (during the winter) since the effects of lateral fluxes appear when s between Θ s and Θ fc (Supporting information S5) globally decreases the amount of water available to plants ( for the metamorphic or sand strata where Ks is the h lateral fluxes were amplified by the high precipitation rates and steep slopes of the Vosges Mountains (Table 3). Calcareous areas show little difference between SWC and SWC to the similar values for Θ s Fig. 5: Monthly SWC (%) calculated with SWB and SWB according to geological and topographic strata. The median is represented by a solid line, while the transparent area illustrates the range from the 1 each strata.

Comparison between SWB and SWB F maps

indices representative for the 1961-1990 period were calculated monthly and mapped at a 50-m resolution. Average SWC values vary from 45% for calcareous soils in winter, to 25% for sandstone in summer (Fig. 5). The distribution of shows considerable spatial variations both at the global and local scale, with lower values in mountains and in calcareous plateaus and near steep areas (Fig. 6). Lateral redistribution along the topography logically leads to a more important decrease of than in the valleys (Table 3, Fig. 5, Supporting information S5 differences occur between SWC and SWC f when the soil is saturated and PET is lower than P (during the winter) since the effects of lateral fluxes appear when s Supporting information S5). Including lateral fluxes in the calculation globally decreases the amount of water available to plants (Table 3, Figs. 5 and7), mainly for the metamorphic or sand strata where Ks is the highest (Table 2). For these two strata, lateral fluxes were amplified by the high precipitation rates and steep slopes of the Vosges . Calcareous areas show little difference between SWC and SWC s and Θ fc and the low Ks that does not allow large fluxes (Table 2 Monthly SWC (%) calculated with SWB and SWB f models for the period 1961 according to geological and topographic strata. The median is represented by a solid line, while the transparent area illustrates the range from the 1 st to the 3 14 1990 period were calculated m resolution. Average SWC values vary from 45% for summer (Fig. 5). The distribution of shows considerable spatial variations both at the global and local scale, with lower values in mountains and in calcareous plateaus and near steep areas (Fig. 6). Lateral a more important decrease of water Supporting information S5). No when the soil is saturated and PET is lower than P (during the winter) since the effects of lateral fluxes appear when soil moisture ranges . Including lateral fluxes in the calculation Figs. 5 and7), mainly ighest (Table 2). For these two strata, lateral fluxes were amplified by the high precipitation rates and steep slopes of the Vosges . Calcareous areas show little difference between SWC and SWC f due that does not allow large fluxes (Table 2).

Effect of lateral fluxes to predict tree species distribution

Most of the species distribution models were improved by adding fluxes, regardless of whether using PAW or DE calculations (Table 4). Lateral fluxes have the greatest influence for two species present in wet environments, Alnus glutinosa, and Salix caprea, and five species among the most common in the study area, that represent half of the trees present in the database: Carpinus betulus, Fagus sylvatica, Pinus sylvestris, Quercus petraea and robur. Most of the species whose distribution is the least influenced by lateral fluxes are linked to rich soil and are mainly present in the calcareous plateaus where lateral fluxes are insignificant (Acer sp., Sorbus torminalis, Ulmus minor), except for Sorbus aucuparia present in the Vosges Mountains where P is always greater than PET and SWB f outputs are close to those of SWB.

Discussion

For the first time, a large-scale study combined fine resolution water balance and lateral fluxes to estimate water available to plants. The simulated soil moisture values seem to be in agreement compared to measurements provided by other studies. For example, [START_REF] Ali | Spatial relationships between soil moisture patterns and topographic variables at multiple scales in a humid temperate forested catchment[END_REF] found volumetric SWC values between 10 and 45% in Canada, and [START_REF] Brocca | Soil moisture spatial variability in experimental areas of central Italy[END_REF] between 13 and 52% in central Italy, while the majority of our simulations range between 19 and 45%. By accounting for soil moisture dynamics, we refined existing maps and improved the estimate of the spatial distribution of water available for plants. We showed that including lateral fluxes in SWB significantly changes the soil moisture patterns, with both a large area effect due to climate and soil properties, and a local effect due to the topography. We also showed that the intensity of fluxes differs according to the location and the period, with few changes when the soil water reserve is full (e.g., in the wet mountains), or when soil characteristics do not allow for large fluxes (e.g., on calcareous plateaus). This little lateral redistribution of soil moisture along hillslopes in chalky catchments has already been observed in previous studies [START_REF] Blyth | Can soil moisture be mapped onto the terrain?[END_REF]. Spatial redistribution is also meaningless when available water is low (< Θ fc ), mainly during dry periods in summer (Supporting information S4), in agreement with existing studies showing that the effects of topography become negligible when soil moisture decreases [START_REF] Ridolfi | Stochastic soil moisture dynamics along a hillslope[END_REF][START_REF] Western | Scaling of soil moisture: A hydrologic perspective[END_REF].

We logically observe a decrease in soil water content during summer, but with a lower intensity, as can be expected. Compared to winter, soil moisture is reduced by 10-40% for the driest month, compared to the 50-80% reduction reported in different studies in the same area [START_REF] Granier | Evidence for soil water control on carbon and water dynamics in European forests during the extremely dry year: 2003[END_REF][START_REF] Lebourgeois | Climate-tree-growth relationships of European beech(Fagus sylvatica L.) in the French Permanent Plot Network (RENECOFOR)[END_REF]. This difference can be due to the water that is stored between ϴ fc and ϴ s in our model, making it available for infiltration the following day, or because precipitation intercepted by vegetation and percolation are not accounted for. Improvements in the model's outputs could probably be obtained by further investigating the partitioning between vertical fluxes and lateral fluxes since the importance of vertical fluxes is frequently reported in the literature [START_REF] Grayson | Terrain and the distribution of soil moisture[END_REF][START_REF] Selby | Hillslope Materials and Processes[END_REF]. However, these values should be compared with caution because soil moisture in our study is averaged over 30 years, includes different soil layers up to 1 m, and is compared to plot measurements or estimations at a specific place, over a short time period and at a specific soil depth. As a follow-up to this work, the model should be validated with ground measurements that could be made at a more local scale and considering shorter time steps.

In contrast to SWC, SWC f shows a large gradient in water content along hillslopes, with between 5 to 25% less SWC f in the crests than in the valleys for metamorphic, sand and sedimentary bedrocks (Fig. 7). This soil moisture gradient was expected because it was already reported in numerous studies [START_REF] Famiglietti | Variability in surface moisture content along a hillslope transect: Rattlesnake Hill, Texas[END_REF][START_REF] Lookingbill | An empirical approach towards improved spatial estimates of soil moisture for vegetation analysis[END_REF]. The amount of water redistributed along the topographical gradient is difficult to compare between studies due to the considerable heterogeneity in climate, soils, topography, periods of the year, time step and measurement methods. Different authors reported that SWC is 15 to 90% lower in the crests compared to the valleys at the same period, with values higher than 40% for wet areas or wet periods [START_REF] Blyth | Can soil moisture be mapped onto the terrain?[END_REF][START_REF] Western | Spatial correlation of soil moisture in small catchments and its relationship to dominant spatial hydrological processes[END_REF], and ranging between 15 and 25% when the climate is dry [START_REF] Blyth | Can soil moisture be mapped onto the terrain?[END_REF][START_REF] Famiglietti | Variability in surface moisture content along a hillslope transect: Rattlesnake Hill, Texas[END_REF]. Compared to the literature, the differences we observed seem relatively small, despite the fact that they are probably a bit underestimated in Fig. 7 by the averaging of SWC f values from many plots with different environmental characteristics. Locally, differences exceeding 50% can be observed on steep slopes or sand strata, but they concern limited areas. These discrepancies may be caused by differences between our estimates based on the entire soil profile, and the shallow measurement depths (0-15 cm) used in some studies since the topographic control on SWC is greater near the surface than at greater depths [START_REF] Famiglietti | Variability in surface moisture content along a hillslope transect: Rattlesnake Hill, Texas[END_REF][START_REF] Grayson | Terrain and the distribution of soil moisture[END_REF][START_REF] Wilson | Identifying and quantifying sources of variability in temporal and spatial soil moisture observations[END_REF]. Additionally, our model provides monthly time step data, whereas studies of spatial patterns in soil moisture are often characterized by daily or hourly measurements of SWC during storm events or intensive rains. Finally, the disaggregation method used to generate daily data for P that is representative of the 1961-1990 period leads to obtaining a representative number of days with rainfall events of the same intensity. By reducing the likelihood of exceeding the thresholds that trigger lateral flow, we probably limit the intensity of runoff. Further work should evaluate the results obtained by using observed sequences of rainfall.

The largest magnitude of subsurface fluxes and the strongest effects of topographic gradients in SWC f are observed in the sand strata that have the highest transmissivity values. This result highlights the model outputs' sensitivity to soil textural assumptions and corresponding PTFs. The range of available PTFs is highly variable among similar soil textural classes, and different PTFs for any given soil texture could result in very different outcomes. For example, the K s value we used for sand is provided by Rosetta and is 3.739 m day -1 (Schaap et al., 2001), but other PTFs use K s values for sand of 5.315 m day -1 [START_REF] Carsel | Developing joint probability distributions of soil water retention characteristics[END_REF] and 0.690 m day -1 (Hypres, [START_REF] Wosten | Development and use of a database of hydraulic properties of European soils[END_REF]. Running the model using the Carsel and Parrish database should consequently increase the lateral fluxes in the sand strata, although they will be significantly decreased using the Hypres database. The limited fluxes observed in the calcareous strata are also explained by PTFs (mainly Ks, Θ s and Θ fc ), and using different PTFs will change their magnitude in this area. This comparison emphasizes the importance of the PTF selection and the high variability that can be generated according to the database used [START_REF] Loosvelt | Impact of soil hydraulic parameter uncertainty on soil moisture modeling[END_REF]. Using local PFT and improving the models' inputs by using continuous PTFs should be a solution for reducing these uncertainties [START_REF] Nemes | Functional evaluation of pedotransfer functions derived from different scales of data collection[END_REF]. Despite the fact that numerous studies have focused on the effects of hydrologic processes on plants [START_REF] Porporato | Plants in water-controlled ecosystems: active role in hydrologic processes and response to water stress -III. Vegetation water stress[END_REF], very few of them were interested in evaluating the effect of lateral fluxes on vegetation over broad areas [START_REF] Berges | Revisiting the use of soil water budget assessment to predict site productivity of sessile oak (Quercus petraea Liebl.) in the perspective of climate change[END_REF][START_REF] Hwang | Ecosystem processes at the watershed scale: Hydrologic vegetation gradient as an indicator for lateral hydrologic connectivity of headwater catchments[END_REF]. The studied tree species distributions were better explained when accounting for the lateral fluxes, demonstrating that soil moisture dynamics contribute to shaping species distribution patterns at a regional scale. The importance of topographical position on vegetation has been recognized for a long time by practitioners, and is commonly linked both to changes in soil depth, water or nutrient fluxes, and local climates [START_REF] Cajander | The theory of forest types[END_REF]. By controlling changes in soil depth, soil nutrition or temperature, we showed that the observed model improvement is not induced by these three environmental gradients, that can be linked with the topographical position [START_REF] Burke | Spatial variability of soil properties in the shortgrass steppe: The relative importance of topography, grazing, microsite and plant species in controlling spatial patterns[END_REF][START_REF] Dyer | Assessing topographic patterns in moisture use and stress using a water balance approach[END_REF]. The better performance of models accounting for the lateral fluxes can probably be attributed to the redistribution of water from the alluvial zone for the species that occur in wet areas, like Alnus glutinosa or Salix caprea. For the other species, soil moisture dynamics along hillslopes are probably the best explanation for these results. The species whose distributions are the least improved are mainly distributed on calcareous bedrocks where lateral fluxes are weak, or in mountains where water availability is not a limiting factor. We can expect that the importance of runoff in shaping species distribution would be underestimated in this study. A large proportion of the studied area is flat (35% is less than 3˚, and 92% is less than 10˚), and most of the steep slopes are located in the Vosges Mountains where little difference exists between SWB and SWBf. Moreover, water constraint is relatively low in this area, and we can expect lateral runoff to have a greater effect on plants in arid or Mediterranean ecosystems [START_REF] Porporato | Plants in water-controlled ecosystems: active role in hydrologic processes and response to water stress -III. Vegetation water stress[END_REF].

Identifying the driest and the wettest areas at the landscape scale is of crucial interest for plant ecologists and natural resources managers. We provided a new tool that makes it possible to study the redistribution of water along hillslopes over large areas. The program, that does not need calibration at the catchment scale, requires only readily available data. It can be used to characterize the spatial distribution in available water over long periods, or to simulate the evolution of soil moisture distribution among years. This information can be used for a wide variety of applications: to refine the delineation of areas with favorable ecological conditions for species, to better predict species composition, to improve our knowledge about species ecology, and to better understand tree growth dynamics, plant health and forest fire risks. It also has important implications for assessing ecosystem responses to climate change by determining topographic refugia where species could persist locally despite regional climatic conditions that may be unfavorable in the future or local areas where species could colonize [START_REF] Ashcroft | Moisture, thermal inertia, and the spatial distributions of near-surface soil and air temperatures: Understanding factors that promote microrefugia[END_REF].

CN III = (23 * CN II ) / (10 + 0.13 * CN II ) (S2.eq2) where CN is the curve number for moisture conditions I, II and III, referring to dry, normal and wet conditions, respectively. CN I and CN III values were used to calculate S in Equation (S2.eq4).

CN adjustment for different slope conditions

In the traditional method, CN II is assumed to be used for 5% slopes. To adjust CN for different slopes, we used the following equation [START_REF] Williams | The EPIC model, Computer models of watershed hydrology[END_REF]:

CN II slope = (CN III -CN II ) /3*(1-2*exp(-13.86*Slope)) + CN II (S2.eq3) where:

-CN II slope is the curve number for normal conditions adjusted for slope, -CN II is the curve number for normal conditions, -CN III is the curve number for wet conditions, -Slope is the slope at the ground surface (in %), CN adjusted for slopes were then recalculated by replacing CN II with CN II slope in Equations C.1 and C.2.

CN adjustment for soil moisture changes

Once CN are adjusted for different moisture conditions and slopes, we allowed the retention parameter S (Equation 5) to vary continuously with soil profile water content by applying the equation used in the hydrological model SWAT (Soil and Water Assessment Tool) [START_REF] Arnold | Large area hydrologic modeling and assessment part I: Model development1[END_REF]

: S = S max * (1 -SW/[SW + exp(w 1 -w 2 * SW)] ) (S2.eq4)
where: -S max is the maximum potential retention solved for dry conditions (i.e., CN I slope ) using equation 5, -SW is the soil water content of the entire profile excluding the water held at wilting point, -w 1 and w 2 are shape coefficients related to maximum potential retention S max , potential retention for wet conditions S 3 , soil water content at saturation SAT, and soil water content at field capacity FC: The first coefficient w 1 and the second w 2 are solved in Equation 5 by assuming that: the retention parameter for moisture conditions I corresponds to soil profile water content at the wilting point; and the retention parameter for moisture conditions III corresponds to soil profile water content at field capacity. These assumptions make the water retention parameter more dependent on soil water storage. Equation (S2.eq4) then replaces Equation 5 to estimate S.

I a /S ratio adjustment

I a represents all losses before runoff begins such as interception by vegetation, water retained in surface depression storage, evaporation and infiltration. Historically, the percentage of initial abstraction was derived from the study of many small, experimental U.S. watersheds, and may therefore not be appropriate in different situations. If needed, the SCS-CN user manual recommends using a relationship other than I a = 0.2S. [START_REF] Woodward | Runoff curve number method: examination of the initial abstraction ratio[END_REF] found that I a /S ratios of 0.05 would seem more appropriate. We used this 5% I a /S ratio to generate (more) surface runoff for smaller amounts of precipitation. Equation 4 becomes [START_REF] Woodward | Runoff curve number method: examination of the initial abstraction ratio[END_REF]:

Q surf = (P -0.05S 0.05 )²/(P -0.95S 0.05 ), if P > 0.05S (S2.eq7) Q surf = 0, if P ≤ 0.05S where S 0.05 = 1.33S 0.20

.15

Text S3: Topographic strata calculation Four topographic strata, including "crest", "mid-slope", "valley", and "flat", were calculated using Jenness Enterprises' Land Facet Corridor Designer (Jenness et al., 2013) for ArcGIS 10 Four-category Topographic Position Tool and a 50-m digital terrain model input. The tool uses a topographic position index (TPI), which is calculated as the difference between the elevation of a pixel and the average elevation of all the pixels in the surrounding neighborhood to assign pixels to different topographic positions (Jenness et al., 2013). Positive TPI values indicate that the pixel elevation is higher than the surrounding neighborhood, and negative values indicate that the elevation is lower than the surrounding neighborhood. If the TPI of a pixel is greater than a given crest threshold or lower than a given valley threshold, it will be classified as a crest or valley. If the TPI is near zero, the elevation of the pixel is similar to that of the surrounding cells and may be located in a flat area or mid-slope. The slope of the pixel is used to differentiate the two situations, and the slope threshold is the minimum slope required for a pixel to be considered mid-slope. Pixels with slopes lower than the slope threshold can be classified as flat, crest or valley.

Each pixel in our study area was assigned a topographic position based on its classification calculated using a 300-m circular neighborhood with crest and valley thresholds of +10 m and -10 m, respectively, and a 3° slope threshold (Class 1 ), and also using its classification based on a 1000-m neighborhood with thresholds of ±5 m and a 3° slope (Class 2 ). Pixels were assigned to crest and mid-slope positions based on their classification in Class 1 and to valley and flat positions based on their classification in Class 2 , given that crest and mid-slope classifications from Class 1 took priority, and Class 2 could only be used to classify remaining pixels as flat or valley. Class 2 was calculated with a larger neighborhood and more sensitive valley threshold so that pixels in broad valleys, which are in a water receiving position at the catchment scale but have a relatively similar elevation to their neighboring pixels within 300 m (and are thus classified as flat in Class 1 ), would be assigned to the valley classification (rather than flat). Areas classified as "flat" were not considered because the elevation in the surrounding neighborhood was so similar that flat pixels were not expected to generate or receive lateral fluxes.

Text S4: Lateral flux estimation

Estimation of lateral fluxes (in mm of water) according to the topographical position: inputs increase from the crest to the valleys, while outputs show the opposite pattern. The net value results from the difference between inputs and outputs. We showed that fluxes mainly occur during spring and autumn, and are very weak during summer due to soil water content < Θ FC .

  Fig.6: July PAW f (mm) estimated with lateral fluxes, for the period1961-1990.

Fig. 7 :

 7 Fig.7: Difference in July PAW estimated with and without fluxes (PAW f -PAW, mm), for the period 1961-1990.

Figure S4. 1 :

 1 Figure S4.1: Monthly subsurface lateral fluxes (mm) for the 1961-1990 period according to topographical strata (n = 1400 for each strata).

Figure S5 .

 S5 Figure S5. Difference between DE and DE f

Table 1 :

 1 List of abbreviations.

	Abbreviation	Description	Units
	SWB	soil water balance	mm
	P	precipitation	mm
	PET	potential evapotranspiration	mm
	AET	actual evapotranspiration	mm
	DE	deficit of evapotranspiration	mm
	RAW	runoff available soil water	mm
	PAW	plant available water	mm
	D	soil thickness	m
	SWDC	soil water draining capacity (Θ sat -Θ fc ) * D	mm
	SWHC	soil water holding capacity (Θ fc -Θ pwp ) * D	mm
	SWC	soil water content (PAW + RAW)	mm
	Volumetric SWC	volumetric SWC (SWC/D)	cm 3 /cm 3
	FC	soil water content at field capacity	mm
	SAT	soil water content at saturation	mm
	Ksat	hydraulic conductivity at saturation	m.d -1
	Θ	volumetric soil water content	cm 3 /cm 3
	Θ pwp	volumetric water contents at permanent wilting point	cm 3 /cm 3
	Θ fc	volumetric water contents at field capacity	cm 3 /cm 3
	Θ sat	volumetric water content at saturation	cm 3 /cm 3
	Q surf	surface runoff	mm
	Q sub	lateral subsurface runoff	mm
	I a	initial abstraction	mm
	S	potential retention	mm
	S max	maximum potential retention	mm
	CN	curve number	---
	V	flow velocity (0. 02 m.s -1 <V< 2 m.s -1 )	m.s -1
	Q out,surf	surface runoff discharge	m 3 .s -1 or mm
	n	Manning's roughness coefficient	s.m -1/3
	R	hydraulic radius at cell i	m
	Slope	slope at ground surface	m.m -1
	B	flow width (cell width)	m
	A i	upstream drainage area at cell i	m²
	a	network constant (2.4 10 -4 )	---
	b	geometric scaling exponent (0.5)	---
	D s	saturated depth area	m
	w	flow width (dimension of the cell)	m
	tan β	land slope	m/m
	d i	fraction of the discharge from a particular cell	---
	L i	effective contour length of cell i: 0.5 and 0.354 for
		downslope cells in cardinal directions and diagonal
		directions, respectively
	e	maximum downslope gradient
	DTW	depth-to-water index	m
	Σ (dz/dx)	the cumulative slope (sum of slope values) along the	m
		least cost path connecting any point of the landscape
		to a watercourse
	a	a is a multiplier equal to 1 when the path is in the	---
		cardinal direction, and 1.414214 when it is diagonal
	w c	grid cell size	m

Table 2 :

 2 Mean

values of model inputs linked to soil characteristics, by geological and topographic strata (n=27311). For abbreviations, see Table

1

.

Table 3 :

 3 Mean monthly values of climatic inputs, lateral fluxes, and PAW

	B. Ks (m day -1 )
	Geology	Crest	Mid-slope	Valley
	Calcareous	0.19	0.19	0.18
	Metamorphic	0.73	0.78	0.78
	Sand	2.08	1.28	1.97
	Sedimentary	0.57	0.46	0.48

f , by geological and topographic strata

(in mm, n=27311)

. For abbreviations, see Table

1

. The winter period includes December, January and February, and the summer period includes June, July and August.

Table 4 :

 4 Changes in performance for the 19 studied species distribution models when adding lateral fluxes: D² for the models without lateral fluxes, difference in D² when adding lateral fluxes and significance of this difference. All the models include T and pH in addition to PAW or DE calculated for

	PAWgs	DEgs

Acknowledgments

This work was supported by the French National Research Agency through the Laboratory of Excellence ARBRE , the GIP ECOFOR and the RMT AFORCE , the Regional Council of Lorraine , and the Direction Régionale de l'Alimentation, de l'Agriculture et de la Forêt of the Lorraine region. We also thanks François Ridremont for its proofreading and Bruno Ambroise for its help.

Supporting Information for

Including the lateral redistribution of soil moisture in a supra regional water balance model to better identify suitable areas for tree species

Text S1: Temporal disaggregation of monthly rainfall

The variation of rainfall event with statistics similar to the 1961-1990 period disaggregated into sequences of daily rainfall occurrences by proposed by [START_REF] Haan | A Markov chain model of daily rainfall[END_REF]. A that the climate alternates between two states S to the other is controlled by transitional probabilities estimated from historical data for each month of a year. probability of the occurrence of an event (or state) in a certain time step depends on the state o previous time step [START_REF] Hoel | Introduction to stochastic processes[END_REF] The state of the climate at time t S(t) is simulated by a Markov process by: S(t)|S(t-1) ~ Markov (Tr, where, S(t) is the state of the climate at day t, S(t 2x2 matrix of transition probabilities whose elements are defined by P ij = Pr(S(t) = i|S(t and where P 1 denotes the "first order" Markov process

Figure S1.1 Graphical model of states of the Markov process and transitional probabilities associated matrix of transition probabilities.

If the transition probability to the state; if not, it switches to the other state. Thus, if the transition probabilities to other states are sufficiently low, the climate persists in one state for a number of days,

Temporal disaggregation of monthly rainfall

events within a month is introduced by generating sequences of wet and dry days 1990 period. Each monthly rainfall data averaged for the 1961 disaggregated into sequences of daily rainfall occurrences by using a first-order two-state Markov Chain model as stochastic model is used for simulating daily rainfall occurrences by assuming that the climate alternates between two states S, either 'dry' or 'wet', and that the transition at time t from one state to the other is controlled by transitional probabilities estimated from historical data for each month of a year.

occurrence of an event (or state) in a certain time step depends on the state o , 1972). The state of the climate at time t S(t) is simulated by a Markov process by: 1) ~ Markov (Tr, P 1 ) here, S(t) is the state of the climate at day t, S(t-1) the state of the climate at the previous time step t 2x2 matrix of transition probabilities whose elements are defined by P ij such as:

= Pr(S(t) = i|S(t-1) = j) i, j = wet or dry, (S1.eq2) denotes the "first order" Markov process

Graphical model of states of the Markov process and transitional probabilities associated

If the transition probability to the other state is lower than the hazard, the climate remains in the current if not, it switches to the other state. Thus, if the transition probabilities to other states are sufficiently low, the climate persists in one state for a number of days, leading to a sequence of wet spells that 29 by generating sequences of wet and dry days Each monthly rainfall data averaged for the 1961-1990 period was state Markov Chain model as stochastic model is used for simulating daily rainfall occurrences by assuming the transition at time t from one state to the other is controlled by transitional probabilities estimated from historical data for each month of a year. The occurrence of an event (or state) in a certain time step depends on the state of the system in the (S1.eq1) 1) the state of the climate at the previous time step t-1, and Tr, a Graphical model of states of the Markov process and transitional probabilities associated with the other state is lower than the hazard, the climate remains in the current if not, it switches to the other state. Thus, if the transition probabilities to other states are sufficiently low, the that alternates with dry spells.

The amount of rainfall on rainy days is assumed to be constant and simply equal to monthly rainfall P divided by the number of simulated wet days within the month.

The Markov model requires knowledge of two monthly conditional probabilities to simulate sequences of rainy and no-rain days over a month: the probability P DW that a wet day occurs at time t, knowing that a dry day occurred at the precedent time t-1, and the probability P WD that a dry day occurs at time t, knowing that a wet day occurred at the precedent time t-1. Both can be estimated from the monthly probability density function (PDF) f M (t) of times of inter-arrivals, where the subscript M refers to the month considered, and t represents the time of interarrivals, i.e., the time elapsed between two successive events of interest. In our case, events are either rainy days P DW = f M (1) or dry days, P WD = f M (1). In order to estimate the PDF f M (t), we used available daily precipitation, and for each month M of the year over the period 1961-90, the within-month distribution of the time elapsed between two successive rainy days or dry days was fitted by testing three PDF: Exponential, Gamma and Weibull using a twosample Kolmogorov-Smirnov test. Whenever the exponential law proved to be in agreement with the data distribution, it was selected because of its ease of use. For each month, P DW and P WD were computed according to its corresponding PDF, and mapped at a 50-m resolution by kriging. Then, for a given month and for each unique couple of transition probabilities P DW and P WD , 5,000 first-order Markov chain models produced 5,000 sequences of states of the Markov model. Each sequence of observations consisting of successive weather, i.e., "wet, wet, dry" (two successive rainy days followed by a day with no rain), has a probability to be observed, which is calculated as:

The most likely sequence is defined as the sequence with the highest probability (to occur), and is used for generating daily maps of precipitation. The CN values for normal moisture condition cover are given in Table S2.1.

Table S2.1: Look-up table between curve number and land cover for the the USDA Soil Conservation Service, 1972

This classical SCS-CN calculation was modified according to the purpose of the study, to facilitate the calculation, to account for the high variability in topography and land use, and to be linked

CN adjustment for AMC

For changes in soil moisture conditions, CN Condition indices (AMC-index): AMC referring to normal and wet conditions depending on the quantity of precipitation calculation we used the following conver Table S2.1 (USDA Soil Conservation Service CN I = (4.2 * CN II

The Curve Number calculation.

The curve number CN used in Equation 5(USDA, 1972) is estimated by: i) Dividing the soils into four classes or hydrologic groups according to their infiltration capacity: low runoff potential] (K sat > 0.18 m -d ); medium capacity [moderate runoff potential] (0.09 m -d < K sat < 0.18 m -d ); apacity [high runoff potential] (0.02 m -d < K sat < 0.09 m -d ); and very low capacity [very high runoff potential] (K sat < 0.02 m -d ).

ii) Assigning specific values to each of the previous classes with respect to land cover. Land cover data was CORINE project (European Environment Agency, 2000).

The CN values for normal moisture conditions (CN II ) used for the four soil classes and the different land up table between curve number and land cover for the four hydrologic soil groups Conservation Service, 1972).

CN calculation was modified according to the purpose of the study, to facilitate the variability in topography and land use, and to be linked For changes in soil moisture conditions, CN II values are adjusted depending on AMC-I referring to dry conditions; and AMC-II (Table referring to normal and wet conditions, respectively. In the traditional method, CN is adjusted on the quantity of precipitation that fell in the last five days and the season. In order to simplify this calculation we used the following conversion formulas adjusted depending on Antecedent Moisture (Table S2.1) and AMC-III is adjusted for a particular AMC days and the season. In order to simplify this to use only the CN values displayed in Θ sat , and Θ fc were estimated with the (Al Majou et al., 2008) class pedotransfer function (class-PTF), and K sat was estimated using Schaap et al. (2001) class PTFs.

Table S7: Look-up table between land cover and roughness coefficient n used in

Manning's roughness coefficient (n, s.m tables reported in the literature [START_REF] Chow | CORINE (Coordination of Information on Environment) Database, a key database for European integrated environmental assessment. Programme of the European Commission[END_REF] up table between land cover and roughness coefficient n used in Equation 6a.

Manning's roughness coefficient (n, s.m -1/3 ) obtained for Corine Land cover codes, by combining multiple value [START_REF] Chow | CORINE (Coordination of Information on Environment) Database, a key database for European integrated environmental assessment. Programme of the European Commission[END_REF]Grimaldi et al., 2010).

quation 6a.

) obtained for Corine Land cover codes, by combining multiple value