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

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-km2 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.


Abstract
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  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.

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 (Porporato et al., 2004). 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 (Barbour and Billings, 2000;Botkin and Keller, 1995;Chabot and Mooney, 1985). This knowledge is particularly important in the current climate change context, with an expected decrease in water availability in large parts of the world (Bates et al., 2008). 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 (Choisnel, 1992), 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 (Dyck, 1985;Saxton, 1985). 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 (Thornthwaite, 1948). 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 (Granier et al., 1999), 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 (Bruand et al., 2003). 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 (Thornthwaite, 1948).
Numerous water balance models have been developed at various time scales (e.g., hourly, daily, monthly and yearly) and to varying degrees of complexity (Xu andSingh, 1998, Schwärzel et al 2011). Most existing PAW maps available over broad areas have been comprised using simplified equations (Van der Schrier et al., 2006;Zierl, 2001). They are often based on models at the monthly scale pioneered in the middle of the last century by Thornthwaite and Palmer (Palmer, 1965;. 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 (Brocca et al., 2007;Ticehurst et al., 2003), with large variations depending on the season and precipitation (Weyman, 1973). The effect of topographical position has also been observed on vegetation. Several studies attributed changes in species composition (Deblauwe et al., 2008;Johnson et al., 2007) and productivity (Berges et al., 2005;Curt et al., 1996;Kobal et al., 2015) 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 (Beven and Kirkby, 1979), TOPOG (Oloughlin, 1981), WET (Moore et al., 1993) and SMR (Frankenberger et al., 1999), 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 (Beven, 1991;Xu and Singh, 1998). 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.
2.1 Input data 2.1.1 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 (Turc, 1961). 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 (Quintana-Seguí et al., 2008;Vidal et al., 2010) and solar radiation from the Helios model (Piedallu and Gégout, 2007). For P, we extracted averaged monthly values over the period 1961-1990 from a 1km resolution map (Bertrand et al., 2011) and occurrences of daily rainfall events used for disaggregation from SAFRAN data (Quintana-Seguí et al., 2008).

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 (Drapier and Cluzeau, 2001). 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  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 . 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 (McBratney et al., 2003). 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.

Soil water balance calculation
Our calculations are based on Thornthwaite's model Thornthwaite and Mather, 1957) that provides monthly estimations of PAW (Piedallu et al., 2013): then PAW(t) = minimum value between PAW(t-1) + P(t) -PET(t) and SWHC (1) and AET equals PET (1b) and PET t -AET t , representing the deficit of evapotranspiration DE t , is equal to zero.
If P(t) < PET(t): then PAW(t) = maximum value between PAW(t-1) 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  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;USDA, 1986): -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) 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: 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 Maidment et al., 1996), -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 (Darcy, 1856) and a kinematic approximation of hydraulic gradient i.e., the hydraulic gradient is equal to the land slope at the site (Beven, 1981;Frankenberger et al., 1999). 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: 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 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 Qin et al. (2007): 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 (Quinn et al., 1991).
The lateral inflow that a cell receives from its surrounding upslope neighbors is defined as: 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 (Murphy et al., 2009): 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.
-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  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 (Drapier and Cluzeau, 2001). 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) (Hastie and Tibshirani, 1990). Mean annual temperature for the 1961-1990 period (Bertrand et al., 2011) 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.

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 3.1 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 (  (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 and 7), 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). models for the period 1961-1990, according to geological and topographic strata. The median is represented by a solid line, to the 3 rd quartile; n=350 for

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, Ali et al. (2010) found volumetric SWC values between 10 and 45% in Canada, and Brocca et al. (2007) 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 (Blyth et al., 2004). 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 (Ridolfi et al., 2003;Western et al., 2002).
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 (Granier et al., 2007;Lebourgeois et al., 2005). 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 (Grayson and Western, 2001;Selby, 1993). 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 (Famiglietti et al., 1998;Lookingbill and Urban, 2004). The amount of water redistributed along the topographical gradient is difficult to 20 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 (Blyth et al., 2004;Western et al., 2004), and ranging between 15 and 25% when the climate is dry (Blyth et al., 2004;Famiglietti et al., 1998). 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 (Famiglietti et al., 1998;Grayson and Western, 2001;Wilson et al., 2004). 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 , but other PTFs use K s values for sand of 5.315 m day -1 (Carsel and Parrish, 1988) and 0.690 m day -1 (Hypres, Wosten et al., 1999). 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 (Loosvelt et al., 2011). Using local PFT and improving the models' inputs by using continuous PTFs should be a solution for reducing these uncertainties (Nemes et al., 2003).
Despite the fact that numerous studies have focused on the effects of hydrologic processes on plants (Porporato et al., 2001), very few of them were interested in evaluating the effect of lateral fluxes on vegetation over broad areas (Berges and Balandier, 2010;Hwang et al., 2012). 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 21 (Cajander, 1926). 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 (Burke et al., 1999;Dyer, 2009). 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 (Porporato et al., 2001).
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 (Ashcroft and Gollan, 2013).

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

Contents of this file
Text S1: Temporal disaggregation of monthly rainfall   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 Haan et al. (1976). 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 (Hoel et al., 1972) 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 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: P(WWD..) = P(W (t) ) * P(W (t) |W (t-1) ) * P(D (t) |W (t-1) ) … (S1.eq3) The most likely sequence is defined as the sequence with the highest probability (to occur), and is used for generating daily maps of precipitation.

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 (Williams and Singh, 1995): 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) (Arnold et al., 1998): 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. Woodward et al. (2003) 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 (Woodward et al., 2003): 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 1 .15 33 Text S3: Topographic strata calculation Four topographic strata, including "crest", "mid-slope", "valley", and "flat", were calculated using Jenness Enterprises' Land Facet Corridor Designer  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 . 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 . Figure S4.1: Monthly subsurface lateral fluxes (mm) for the  period according to topographical strata (n = 1400 for each strata).

Figure S5. Difference between DE and DE f
Difference in July DE estimated with and without fluxes (DE f -DE, mm), for the period 1961-1990. 36 (Chow et al., 1988; up table between land cover and roughness coefficient n used in Equation 6a.
) obtained for Corine Land cover codes, by combining multiple value