Mael Le Corre

Maryline Pellerin

David Pinaud

Guy Van Laere

Hervé Fritz Ac

Dr Sonia Saïd
email: sonia.said@oncfs.gouv.fr

A multi-patch use of the habitat: testing the First-Passage Time analysis on roe deer paths

Keywords: area-restricted search, first-passage time analysis, roe deer, spatial heterogeneity

A heterogeneous environment includes several levels of resource aggregation. Individuals do not respond in the same way to this heterogeneity depending on the scale at which they perceive it, and develop different foraging tactics accordingly. The development of methods to analyse animal movements has enabled the study of foraging tactics at several scales.

, to detect a multi-patch use of the habitat. First-Passage Time analysis results showed that most of the female roe deer exploited their home range as a patchy resource, ranging from 1 to 5 areas of intensive use in their home range. These areas were identified as the most attractive sites within the roe deer female home range. Moreover, this method allowed us to rank the attractive areas according to the time spent in each area.

Coupled with habitat selection analysis to identify what makes these areas attractive, the First-Passage Time analysis should offer a suitable tool for landscape ecology and management.

Introduction

Animals live in a patchy and hierarchical structured environment. Resources are aggregated and habitat complexity is generally reflected by this aggregation at different scales, lower levels nesting into upper levels [START_REF] Kotliar | Multiple scales of patchiness and patch structure: a hierarchical framework for study of heterogeneity[END_REF]). According to the scale, resources differ in density, turnover and predictability [START_REF] Wu | From balance of nature to hierarchical patch dynamics: a paradigm shift in ecology[END_REF][START_REF] Fauchald | Foraging in a hierarchical patch system[END_REF][START_REF] Fauchald | Scale-dependent predator-prey interactions: the hierarchical spatial distribution of seabirds and prey[END_REF]. The spatial and temporal distribution of a resource has a dominant influence on the searching efficiency and behaviour of a predator [START_REF] Bell | Searching behavior patterns in insects[END_REF]). Hierarchical organization of resource aggregation thus induces both a differential perception of the habitat by predators according to the scale, and a modification of their search behaviour [START_REF] Fauchald | Foraging in a hierarchical patch system[END_REF]). Animals will change their movement patterns in relation to changes in the spatial distribution of resources [START_REF] Fritz | Scale-dependent hierarchical adjustments of movement patterns in long-range foraging seabird[END_REF]. According to the Optimal Foraging Theory, when resources are distributed in patches, a consumer should increase its intake rate in a high resource density patch [START_REF] Stephens | Foraging Theory[END_REF], slowing down its speed and increasing its turning rate. [START_REF] Kareiva | Swarms of predators exhibit "preytaxis" if individual predators use area-restricted search[END_REF] defined this concentration of consumer activity within an area of interest as an "Area-Restricted Search" behaviour (ARS). In fact, animal movement trajectories should result from the interaction between animal decisions and landscape properties [START_REF] Morales | Adaptive models for large herbivore movements in heterogeneous landscapes[END_REF].

Movement patterns of animals have mainly been analyzed using point distribution and path trajectories [START_REF] Turchin | Quantitative analysis of movements: measuring and modelling population redistribution in animals and plants[END_REF]. Two main but different ways have been used to study vertebrate movements: studies based on random walks [START_REF] Mårell | Foraging and movement paths of female reindeer: insights from fractal analysis, correlated random walks, and Lévy flights[END_REF] or on fractal dimension [START_REF] Nams | Using animal movement paths to measure response to spatial scale[END_REF]. Correlated random walks have been used to model random movements with a certain degree of directionality [START_REF] Kareiva | Analyzing insect movement as a correlated random walk[END_REF][START_REF] Bergman | Caribou movement as a correlated random walk[END_REF], while the fractal dimension measures a path tortuosity. The tortuosity of foraging paths represents animal reactions to landscape heterogeneity in which animals translate environmental stimuli into movements [START_REF] Crist | Animal movement in heterogeneous landscapes: An experiment with Eleodes beetles in shortgrass prairie[END_REF][START_REF] With | Using fractal analysis to assess how species perceive landscape structure[END_REF]. Our understanding of 4 these relationships and processes has heavily relied on studies of small organisms or insects for which the monitoring of movement and the manipulation of the microlandscapes are feasible (e.g. [START_REF] Crist | Animal movement in heterogeneous landscapes: An experiment with Eleodes beetles in shortgrass prairie[END_REF][START_REF] With | Using fractal analysis to assess how species perceive landscape structure[END_REF], With et al. 1999). However, the monitoring of foraging paths of large animals has recently been developed (e.g. [START_REF] Gross | Movement rules for herbivores in spatially heterogeneous environments: responses to small scale pattern[END_REF][START_REF] Etzenhouser | Foraging behavior of browsing ruminants in a heterogeneous landscape[END_REF][START_REF] Viswanathan | Optimizing the success of random searches[END_REF][START_REF] Johnson | Movement parameters of ungulates and scale specific responses to the environment[END_REF][START_REF] Mårell | Foraging and movement paths of female reindeer: insights from fractal analysis, correlated random walks, and Lévy flights[END_REF][START_REF] Fauchald | Using First-Passage Time in the analysis of area-restricted search and habitat selection[END_REF][START_REF] Fritz | Scale-dependent hierarchical adjustments of movement patterns in long-range foraging seabird[END_REF]. The methods developed in these studies are mainly based on the shift between displacement patterns to define spatial scale changes. The study of trajectories allows the understanding of how the habitat heterogeneity affects foraging efficiency of the predators and how they react to changes in the spatial structure of resources.

In a recent study, [START_REF] Fauchald | Using First-Passage Time in the analysis of area-restricted search and habitat selection[END_REF] used variation in the search effort to assess at which scale consumers interact with their environment. The "First-Passage Time" (FPT) analysis, developed by [START_REF] Fauchald | Using First-Passage Time in the analysis of area-restricted search and habitat selection[END_REF], estimates the search effort, measuring the time spent by a predator to cross a circle of a given radius and so, highlights the part of the path where it adopts an ARS behaviour. The FPT analysis provides an estimation of the scale at which a consumer perceives resource by determining the exploited patch size [START_REF] Fauchald | Using First-Passage Time in the analysis of area-restricted search and habitat selection[END_REF].

Determination of the perception scale of resource by an animal is important in understanding exploitation tactics, and this perception often differs from our own [START_REF] Levin | The problem of pattern and scale in ecology[END_REF]). As most of the studies on large animal movements (e.g. [START_REF] Bascompte | Fractals and search paths in mammals[END_REF][START_REF] Viswanathan | Optimizing the success of random searches[END_REF][START_REF] Bergman | Caribou movement as a correlated random walk[END_REF][START_REF] Johnson | Movement parameters of ungulates and scale specific responses to the environment[END_REF][START_REF] Mårell | Foraging and movement paths of female reindeer: insights from fractal analysis, correlated random walks, and Lévy flights[END_REF][START_REF] Fritz | Scale-dependent hierarchical adjustments of movement patterns in long-range foraging seabird[END_REF][START_REF] Austin | Intraspecific variation in movement patterns: modeling individual behaviour in a large marine predator[END_REF], FPT analysis was performed on seabirds with high travelling abilities, making long foraging trips in low predictable resource habitats [START_REF] Fauchald | Using First-Passage Time in the analysis of area-restricted search and habitat selection[END_REF], 2006[START_REF] Pinaud | Scale-dependant habitat use in a long-ranging central place predator[END_REF][START_REF] Pinaud | At-sea distribution and scale-dependent foraging behaviour of petrels and albatrosses: a comparative study[END_REF]. However [START_REF] Frair | Scales movements by elk (Cervus elaphus) in response to heterogeneity in forage resources and predation risk[END_REF] and [START_REF] Bailey | Quantitative analysis of bottlenose dolphin movement patterns and their relationship with foraging[END_REF] also performed FPT at a smaller scale, respectively on a large ungulate 5 (elk, Cervus elaphus) living in a more restricted home range and on groups of bottlenose dolphins (Tursiops truncatus Montagu) in a small surface area. In this study, we also attempted to apply the FPT analysis at a finer scale than in [START_REF] Fauchald | Using First-Passage Time in the analysis of area-restricted search and habitat selection[END_REF] and [START_REF] Pinaud | Scale-dependant habitat use in a long-ranging central place predator[END_REF] in order to investigate the structure of roe deer (Capreolus capreolus) daily paths.

It has been assumed that herbivores generally feed at sites with high food quality and/or quantity, following qualitative predictions from Optimal Foraging Theory (e.g. [START_REF] Stephens | Foraging Theory[END_REF]. These predictions are supported by feeding-patch selection by deer which has been shown to be related to nutrient content and plant biomass [START_REF] Wilmshurst | Patch selection by red deer in relation to energy and protein intake: a re-evaluation of Langvatn and Hanley's (1993) results[END_REF].

Effects of habitat quality on roe deer population dynamics are quite well understood [START_REF] Pettorelli | Population density and small-scale variation in habitat quality affect phenotypic quality in roe deer[END_REF][START_REF] Pettorelli | Variations in adult body mass in roe deer: the effects of population density at birth and of habitat quality[END_REF](Pettorelli et al. , 2003a and b) and b), but less so on the individual use of space [START_REF] Mysterud | Sacle-dependent trade-offs in foraging by European roe deer (Capreolus capreolus) during winter[END_REF]. Though many studies on roe deer have been performed on the relationship between environmental factors and home range variation [START_REF] Tufto | Habitat use and ecological correlates of home range size in a small cervid: the roe deer[END_REF][START_REF] Mysterud | Sacle-dependent trade-offs in foraging by European roe deer (Capreolus capreolus) during winter[END_REF], Saïd et al. 2005, Saïd & Servanty 2005), very few studies have investigated the heterogeneous use of the habitat within the home range.

In this paper, we used data from roe deer females equipped with GPS collars to investigate the structure of their movements throughout a day, using FPT analysis. Roe deer are known to have quite fixed and restricted home ranges [START_REF] Strandgaard | The roe deer (Capreolus capreolus) population at Kalø and the factors regulating its size[END_REF][START_REF] Liberg | The function of male territoriality in roe deer[END_REF]. It is thus particularly challenging to determine how roe deer may perceive or define the spatial organization of their resource in their daily environment and the impact of environmental structure on their movements.

The aim of our study was to see if roe deer perceive their habitat as heterogeneous. We hypothesized 1) that in a heterogeneous environment roe deer will use areas where they will increase their research of palatable plants and slow down their speed to increase their intake rate [START_REF] Fauchald | Using First-Passage Time in the analysis of area-restricted search and habitat selection[END_REF]; 2) roe deer will use several ARS throughout the day as this 6 species presents several activity phases with two main peaks at dawn and dusk [START_REF] Bubenik | Rythme nycthéméral et le régime journalier des ongulés sauvages, Problèmes théoriques, Rythme d'activité du chevreuil[END_REF], Maublanc et al. 1991).

Methods

Study area

This study was carried out in the Réserve Nationale de Chasse et de Faune Sauvage of Chizé (see Figure 1). The Chizé reserve is an enclosed forest of 2614 ha in western France (46°05'N, 0°25'W). The elevation varies between 47 and 101 meters and the oceanic climate is characterized by mild winters and hot and dry summers. The Chizé forest includes three habitats contrasting in quality: an oak (Quercus spp.) forest with resource-rich coppices dominated by hornbeam (Carpinus betulus) in the north-eastern part, an oak forest with coppices of medium quality dominated by Montpellier maple (Acer monspessulanum) in the north-western part, and a poor beech (Fagus sylvatica) forest in the south part (Pettorelli et al. 2003a, Saïd & Servanty 2005; see Figure 1). The roe deer population at Chizé was estimated from Capture-Mark-Recapture methods to be approximately 400 adults in 2003and 450 in 2004(e.g. Gaillard et al. 1992, 1993[START_REF] Pettorelli | Variations in adult body mass in roe deer: the effects of population density at birth and of habitat quality[END_REF], 2003b, Gaillard unpubl. data).

Data collection

Does were equipped with Lotek's GPS_3300 radio collars (Lotek Wireless, Fish & Wildlife Monitoring). These collars, weighing 285 g, provided information on GPS positioning in differential mode (i.e., latitude, longitude, date and time) at pre-programmed intervals, fixed quality (DOP = dilution of precision) and ambient temperature. We scheduled collars programming one location every 5 minutes during 24 hours several days (1 to 3 days) per month. In September to December 2003, 2 females were equipped with GPS collars. 7

Then, in January-February 2004, 7 does were collared, in addition to the 2 previously equipped ones. Data location and data activity were recovered during the capture campaign in January-February 2004 and 2005. We obtained 15 daily paths with one point every 5 minutes from September 2003 to January 2004 and 15 daily paths from May to November 2004 (see Figure 2).

First-Passage Time analysis The First-Passage Time method

To detect one or several ARS, we perform an FPT analysis on our data, following [START_REF] Fauchald | Using First-Passage Time in the analysis of area-restricted search and habitat selection[END_REF] License. The First-Passage Time (FPT) corresponds to the time needed by an animal to cross a circle of a given radius r. In order to perform the analysis, we assumed linear path with a constant speed between locations (e.g., [START_REF] Kareiva | Analyzing insect movement as a correlated random walk[END_REF] and completed interlocation path by points spaced every 5 meters as 90% of the path segments were longer than 5 meters. FPT was calculated at each point of individual paths for a given radius r, the location corresponding to the centre of the circle. We determined this measure of time for r ranging from 5 to 400 meters, as locations were spaced every 5 meters maximum and as the mean radius of the roe deer home range in the Chizé Forest is about 300 meters. The radius where the highest peak of variance in FPT is observed corresponds to the spatial scale (radius) of the area where the individual concentrates its search effort. We thus calculated the variance in FPT given by Var[log t(r)], S(r), for each radius, and plotted it against the radius to obtain the spatial scale corresponding to the peak in S(r), of the globality of individual path. Spatial scale was confirmed when the plot of the FPT against time for the radius for which S(r) was maximum showed intensive search pattern, that is a high and constant FPT (see Fauchald & 8 Tveraa 2003, see Figure 3). Then, the intensive search area which corresponds to an ARS was identified on a path as the circle with the longest FPT for the determined spatial scale.

Multi-patch detection

Using the FPT analysis we detected the area where the search effort was the highest.

Moreover, the circle radius gave us an estimation of the size of the exploited patch. However, the majority of our paths were clearly composed of more than one ARS (see Figure 2 and3).

The scale used to estimate the size of the main exploited patch can not be applied to the other ARS. Indeed, these ARS do not have necessarily the same radius than the one with the highest FPT. In order to assess the scale of these ARS, we need to remove the effect of the main ARS.

We therefore modified the original process described by [START_REF] Fauchald | Using First-Passage Time in the analysis of area-restricted search and habitat selection[END_REF], and programmed by [START_REF] Pinaud | Scale-dependant habitat use in a long-ranging central place predator[END_REF] to detect other potential ARS (program available contacting the authors). Part of the path included in the first determined FPT circle was removed and replaced by a segment that individuals cover in 5 minutes, as if the point before and after the ARS were really two successive points in the path. We then corrected the time of the locations following the removed path part. FPT analysis was finally performed on the modified path and we repeated the analysis until we obtained no more potential ARS (see Figure 4).

ARS validation

All circles for which the FPT plotted against time did not show any intensive search pattern were removed. Remaining circles were re-sorted: circles overlapping previously detected ARS and circles with radius less than 25 meters (GPS error calculated for Chizé forest = 25.56 ± 33.99 meters) were removed, especially in order to minimize the bias induced by the positioning error of the GPS. Finally, circles with a mean speed displacement higher 9 than the mean speed in the whole path were eliminated. Indeed, when applying the method in order to detect other ARS, the substitution of previously detected circles by a straight line covered in 5 minutes sometimes produced a bias in the detection. The identified parts of the path correspond to the parts where individuals spend more time and slow down their speed according to the ARS behaviour. As these parts of the path were substituted by a 5 minutes covered segment, the mean speed displacement along the path increased. A new ARS was detected where individual moved slower than the whole new generated path. However, the speed of the new generated path was higher than the one of the original path. The difference between the mean speed inside the ARS and the mean speed of the original path was smaller than with the generated path where the FPT method was applied. In extreme cases, speed inside the ARS was even higher than the speed of the original path and animals move faster in the detected circle than in the whole original path, that is not congruent with the definition of the ARS. Consequently, such circles were not taken into account. In order to identify these circles, we estimated the individual's mean speed in the ARS and compared it with the speed that the deer could have during the same time period. Given that the part of the path included in the circle is composed of n segments with associated speed, we randomly drew n segments, with replacement, of the whole path and estimated the mean speed. We performed this 100 000 times to obtain a mean speed distribution over the path. If the mean speed in the ARS was higher than the median value of the distribution, we did not take into account the corresponding circle. At last, the selected ARS were divided into 3 classes to assess whether the ARS differed according to the detection order: class 1, the first detected ARS; class 2, the second detected ARS; and class 3, the last ARS.

Statistical analysis 1

To test the effect of female identity, year, month, class of circle and circle radius as a covariable on the FPT in the circle, we fitted mixed linear models ("lme" function in the "nlme" R package, [START_REF] Laird | Random-Effects Models for Longitudinal data[END_REF][START_REF] Pinheiro | Mixed-effects models in S and S-Plus[END_REF][START_REF] Pinheiro | nlme: Linear and nonlinear mixed effects models[END_REF] with the FPT in the circle, the year (2003 or 2004) as a two-modality fixed factor, the month (May, August, September, November and December) as a five-modality fixed factor, the class of circle (1, 2 or 3) as a three-modality fixed factor, and the females identity as random factor.

We did the same with the log-transformed FPT divided by the circle surface (FPT/m²).

This measurement represents the utilization intensity of an area. We fitted mixed linear models with log(FPT/m²), year, month, class of circle (1, 2 or 3) as fixed factors, and only the females identity as random factor. We tested the random effect with a Restricted Maximum Likelihood (REML) procedure and the fixed effect with a Maximum Likelihood (ML) procedure [START_REF] Vaida | Conditional Akaike information for mixed-effects models[END_REF]. To select the best model, we used the Akaike Information Criterion (AIC) [START_REF] Burnham | Model selection and inference: a practical information-theoretic approach[END_REF] and we retained the model with the lowest AIC value. When the difference between two models was less than 2, we retained the simplest model according to the parsimony rules [START_REF] Burnham | Model selection and inference: a practical information-theoretic approach[END_REF]. We checked the normality of residuals of our selected models for FPT (Shapiro-Wilk test: W = 0.994, pvalue = 0.94) and log(FPT/m²) (Shapiro-Wilk test: W= 0.991, p-value = 0.77). Standard deviation of the random factor Female Id and the residuals were respectively for each model: Female Id SD = 0.852, Residual SD = 0.994, and Female Id SD = 0.517, Residual SD = 0.330.

Results

Over the 30 analyzed paths, 24 (13 in autumn/winter, 11 in spring/summer) presented a peak of the relative variance S(r) (See Figure 5) with a radius that seemed to correspond to an ARS. There was a total of 55 ARS detected (Class 1 = 23, Class 2 = 20, Class 3 = 12), 1 ranging from 1 to 5 ARS per path (mean = 2.2, SE = 0.96) (See Figure3). Radius of the detected circles ranged from 25 to 105 meters (mean = 50.1, SE = 20.47m), and the FPT ranged from 1h 38 min to 10h 05min (mean = 4h 13min, SE = 1h 59min). The best mixed linear model for the First-Passage Time included ARS radius (See Figure 6), circle class (See Figure 7) and month (See Figure 8), as well as the female ID (See Table1 and Table 2, AIC = 189.96). We observed a positive relationship between the FPT and the ARS radius (slope = 0.052, See Figure 6).

For log(FPT/m²), the best mixed linear model included circle class (See Figure 7), month (See Figure 8) and female ID (See Table 1 andTable 2, AIC = 75.36). Moreover, the result showed that log(FPT/m²) was negatively correlated with circle class: the use of the patches by the females was more intensive in the first detected ARS than in the second ones and the third ones (See Figure 7).

Discussion

Using the First-Passage Time analysis [START_REF] Fauchald | Using First-Passage Time in the analysis of area-restricted search and habitat selection[END_REF], we detected at least one Area-Restricted Search behaviour (sensu [START_REF] Kareiva | Swarms of predators exhibit "preytaxis" if individual predators use area-restricted search[END_REF] on most of the paths. Our study highlighted the existence of a second scale of landscape use which was finer than the home range scale for the roe deer (0.8 ha versus 24.5 ha for home range in Chizé forest (Saïd & Servanty 2005)). Many authors showed that certain portions within the home range are more frequently used than others [START_REF] Adams | The internal anatomy of home range[END_REF][START_REF] Dixon | Harmonic mean measurement of animal activity areas[END_REF][START_REF] Springer | Movement patterns of coyotes in south central Washington[END_REF][START_REF] Macdonald | Enduring social relationships in a population of crab-eating fox zorro, Cerdocyon thous, in Amazonian Brazil[END_REF]. These areas of concentrated use by resident animals, loosely termed core areas, commonly include nest sites, daytime roost sites, refuges, and regions with the most dependable food sources [START_REF] Burt | Territoriality and home range concepts as applied to mammals[END_REF][START_REF] Kaufmann | Ecology and social behavior of the coati, Nasua nirica on Barro Colorado Island Panama[END_REF][START_REF] Ford | Home range in a patchy environment: optimal foraging predictions[END_REF]). 1 Methods such as harmonic mean [START_REF] Dixon | Harmonic mean measurement of animal activity areas[END_REF] or bivariate models [START_REF] Koeppl | A bivariate home range model with possible application to ethological data analysis[END_REF] were used to estimate core areas based on the cumulative proportion of an animal's location and the related increase of the area estimated by these locations within the total home range. We decided to use FPT method, because it allows identifying smaller sectors where animals spend more time using ARS than the core area method. Although the limited accuracy of the GPS collars did not allow us to explore nested levels of resource aggregation, our study highlights a patchy exploitation of their home range by the roe deer females. A large proportion of their activity is concentrated in restricted parts of their environment. In our study, we focused on the third scale of habitat selection defined by [START_REF] Johnson | The comparison of usage and availability measurements for evaluating resource preference[END_REF] See Figure 9) but more accurate data are needed to highlight and define the finest scale (i.e., selection of food item inside feeding site).

The modification of the original method in order to detect further ARS along the path allowed us to highlight a multi-patch use of the home range at a daily level in roe deer females for 20 out the 24 paths. One major methodological questions associated with the use of our modified version of the FPT is obviously when to stop running the iterative analysis, i.e. how many patches make biological sense or can truly be distinguished with some statistical power. We chose to rely on a very simple and repeatable methodological rule to set our limit: we stopped when the size of the ARS circles fell below the accuracy of the GPS measures (i.e. below 25 m). There certainly are other decision rules, for instance for animals showing a very high selection for one or two patches, and for which the addition of another three minute circles, although significant, may not make much biological sense. We never had to question ourselves about this on our roe deer data, probably because our decision rule fell nicely within the scales at which roe deer make decisions about patch choices. However, we encourage users of this modified method to adapt the decisions rules to their equipment but more so to the biology of their studied animals. 1

In our study, we used the First-Passage Time value (FPT) and the FPT by surface unit.

The first one represents the time an individual stayed in the patch and gives us the time he spent foraging. The second one represents the exploitation effort of an individual in a patch. It indicates the intensity with which an individual foraged in the patch. We found a continuum in the FPT between the circle classes: animals spent more time in circles of first class than in circles of second and third classes. However, this result can be directly linked to the way our extension of the FPT method works. Indeed, as the method focuses on the ARS where the FPT is the maximum, the first detected ARS is systematically the one with the highest FPT, the second detected ARS is the one with the second highest FPT and so on. Nevertheless, the exploitation intensity (log(FPT/m²)) of the patch really differs between circle classes. The time spent in an ARS of a given radius is higher in circles of the first class than in circles of the two other classes. Thus, the detection order seems to reflect an order of preference in the use of patches within the home range and roe deer should spend more time in the most valuable patch.

There could be another way to interpret differences between circle classes. They could reflect a spatial discrimination of different behaviours (e.g., resting and feeding activities) that imply different constraints and resources, or environmental features. In [START_REF] Fauchald | Using First-Passage Time in the analysis of area-restricted search and habitat selection[END_REF], 2006) and [START_REF] Pinaud | Scale-dependant habitat use in a long-ranging central place predator[END_REF], studied paths were foraging trips. They performed respectively the FPT analysis on the Antarctic petrel (Thalassoica antarctica) and the yellow-nosed albatross (Thalassarche carteri) during the breeding period. In [START_REF] Pinaud | At-sea distribution and scale-dependent foraging behaviour of petrels and albatrosses: a comparative study[END_REF], the analysis was also performed on seven breeding species of Procellariforms. During this period, the two species display a central-place forager behaviour [START_REF] Ashmole | Seabird ecology and the marine environment[END_REF]): they make a long foraging trip of several days after they left the colony, and then come back. Birds thus have to maximize their energy gain, spending most of their trip in the food patch so that their activity within the ARS is mainly related to foraging. In our study, 1 like in Frair et al.'s (2005) study, roe deer path did not seem to be typically displaying centralplace foraging as starting point hours were chosen arbitrarily and the path remained all the day including feeding and resting activity for which roe deer look for different structures of vegetation (Mysterud & Øtsbye 1995). [START_REF] Frair | Scales movements by elk (Cervus elaphus) in response to heterogeneity in forage resources and predation risk[END_REF] worked at a large temporal scale using elk (Cervus elaphus) paths taken over several months with a 2 hour inter-location interval. They used the FPT analysis to define 3 scales of movements (inter-patch relocations, foraging displacements and resting movements). In our study we used daily paths, working at a finer temporal scale with a 5 minutes inter-location interval. Performing the FPT analysis on paths with locations every 5 minutes gives us a sharp information of the animal movement and allows us to highlight a succession of used patches at the scale of the day. These patches can also be characterized by activity data, as GPS collars can be used in combination with activity sensors (e.g. [START_REF] Adrados | GPS approach to study fine-scale site use by wild red deer during active and inactive behaviors[END_REF], Frair et al. 2005). The use of such tools to identify activity within FPT circles will provide means to assess whether animals display a multipatch use of the environment because of the resource patch heterogeneity, versus a spatial discrimination between different behaviours (e.g. resting vs feeding sites).

The patch use behaviour reported in this study corresponds, in the case of a food search strategy, to the animal's perceptions of the resource aggregation, and it is probably related to the habitat heterogeneity. This habitat heterogeneity, natural or due to human disturbance, drives ecological processes [START_REF] Fortin | Landscape ecology comes of age[END_REF]. Understanding the link between these processes and the heterogeneity of the environment is a key question for management and landscape ecology [START_REF] Wiens | Ecological mechanisms and landscape ecology[END_REF]. Coupling a Geographic Information System with the FPT analysis would allow us to characterize areas of interest (eg. ARS for feeding, resting site) in terms of environmental features (topography, vegetal cover, thermic cover, vegetation quality) according to animal activity. FPT analyses should offer a suitable tool to be incorporated in habitat selection analysis in order to improve identification of the 1 environmental features that make a habitat attractive to animals. Fine approaches as FPT analysis enable to highlight the response of species to habitat heterogeneity at very small scales. The study of this behavioural response allows assessing in a finest way the consequences of management decisions on population dynamics and then to orientate these decisions to improve the management of these populations and their habitats.

# # # # # # # # # # # # # ## # # # # ### ## # # # # # # # # # # # #### # # # ## # # # # # # # # ## # # ## # # # # # # # ## # # # # # # # # # # # # # # # # # # # # # # ## # # # # # # # # # # # # # ## # # # ## # # # # # # # # # # # # # # # # # # # # # ## ## # # # # # # # # # # # ## # # # ## # # # # # # # # # # # # # # # # # # # # # # # # # # ## # # # # # # # # ### # # ## # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ## 0 1 2 3 4
Oak forest with coppices of hornbeam Oak forest with coppices of Montpellier maple Beech forest

  method, and using the software R (version 2.1.0; R Development Core Team 2005; Ihaka & Gentleman 1996) distributed under the GNU General Public

Figure 1 :

 1 Figure legends:

Figure 2 :

 2 Figure 2: Example of a path where several areas were detected (female 35, the 24th of

Figure 3 :

 3 Figure 3: Example of First-Passage Time all along the path for the scale assessed by the

Figure 4 :

 4 Figure 4: Example of detection of a second area: (a) a first circle was detected after the FPT

Figure 5 :

 5 Figure 5: Mean variance in First-Passage Time S(r) against the radius r of all paths showing

Figure 6 :

 6 Figure 6: First-Passage Time (Hours) as a function of the corresponding radius (meters) of

Figure 7 :

 7 Figure 7: First-Passage Time (Hours) (empty box whisker plots) and Mean log(First-Passage

Figure 8 :

 8 Figure 8: First-Passage Time (Hours) (empty box whisker plots) and Mean log(First-Passage

Figure 9 :

 9 Figure 9: The four scales of habitat selection defined by Johnson (1980), with radius in

  Fig.2 1

  Fig.8 1

Table 1 .

 1 Results of the mixed linear models for the First-Passage Time and for the logtransformed First-Passage Time divided by the circle area of female roe deer paths in Chizé forest, with the female identity as random factor.

	Tested models
	Response variable	Random	Fixed effect	Number of	∆ AIC
		effect		Parameters
			Month		52.65
			Year		48.58
			Constant		46.54
			Circle Class + Month		39.13
	First-Passage Time	Female	Circle Class + Year Circle Class		34.56 32.57
			Radius		29.61
			Radius + Circle Class + Year		7.34
			Radius + Circle Class		6.02
			Radius + Circle Class + Month		0
			Year		41.61
	log(First-Passage Time/ Surface unit)	Female	Year + Circle Class Circle Class Month		31.67 29.99 15.28
			Month + Circle Class		0

Table 2 .

 2 Estimates for the mixed linear models selected in table 1 for the First-Passage Time and for the log-transformed First-Passage Time divided by the circle area of female roe deer paths in Chizé forest, with the female identity as random factor. Standard deviation of random

	factors are respectively for each model: Female Identify SD = 0.852, Residual SD = 0.994,
	and: Female Identify SD = 0.517, Residual SD = 0.330.
	Response	Tested	Explanatory	Value	S.E	D.F	t-value	p-value
	variable	models	Variables
			(Intercept)	+1.42	0.762		+1.864	0.070
			Radius	+0.08	0.010		+7.539	<0.001
			Circle Class 2	-1.66	0.339		-4.895	<0.001
	First-Passage Time	Radius + Circle Class + Month	Circle Class 3 May August	-2.24 -0.42 +0.06	0.444 0.622 0.719		-5.032 -0.674 +0.081	<0.001 0.505 0.936
			October	+0.24	0.784		+0.303	0.764
			November	+0.50	0.804		+0.621	0.539
			December	-1.67	0.732		-2.274	0.029
			(Intercept)	-7.18	0.251		-28.626	<0.001
			Circle Class 2	-0.46	0.111		-4.178	<0.001
	log(First-Passage Time/ Surface unit)	Circle Class + Month	Circle Class 3 May August October	-0.50 -0.37 +0.42 +0.15	0.146 0.238 0.248 0.279		-3.418 -1.562 +1.706 +0.523	0.002 0.126 0.096 0.604
			November	+0.03	0.287		+0.121	0.905
			December	-0.78	0.259		-3.017	0.005

Acknowledgements

This work was made possible by financial support form GIPEcoFor ("Groupement d'Intérêt Public Ecosystèmes Forestiers"). We are grateful to all the students and volunteers for help to collecting data on the study site. We would like to thank Jean-Michel Gaillard, Clément Calenge, Gilles Yoccoz, Philippe Aubry, Bram Van Moorter, Olive Said, François Klein and two anonymous reviewers for valuable discussions and review of earlier versions of this manuscript and Timothee Cook for English revision.


Simone Di Prima
email: sdiprima@uniss.it

Ryan D Stewart

Mirko Castellini

Majdi Vincenzo Bagarello

Abou Najm

Mario Pirastru

Filippo Giadrossich

Massimo Iovino

Rafaël Angulo-Jaramillo

Laurent Lassabatère

Majdi R Abou Najm

Rafael Angulo-Jaramillo

Laurent Lassabatere

Estimating the macroscopic capillary length from Beerkan infiltration experiments and its impact on saturated soil hydraulic conductivity predictions

Keywords: infiltration, macroscopic capillary length, Beerkan, hydraulic conductivity

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

 [START_REF] Gardner | Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a water table[END_REF]proposed an exponential relation between the soil hydraulic conductivity, K (L T - 1 ) and the water pressure head, h (L):

𝐾(ℎ) = 𝐾 𝑠 𝑒𝑥𝑝(αℎ) (1)
where Ks (L T -1 ) is the saturated soil hydraulic conductivity and α (L -1 ) is a fitting parameter. Many authors pointed out the importance of the α parameter. For instance, [START_REF] Philip | Reply To "Comments on Steady Infiltration from Spherical Cavities[END_REF] introduced a parameter related to α called sorptive length that expresses the role of sorption in multidimensional

(2-3D) steady state flows. The reciprocal of α, i.e., the macroscopic capillary length, λc (L), is a more general length scale parameter that was firstly defined as an integral parameter by [START_REF] Bouwer | Unsaturated Flow in Ground-Water Hydraulics[END_REF] (called capillary fringe in his work) and later discussed in detail by [START_REF] White | Macroscopic and microscopic capillary length and time scales from field infiltration[END_REF] and [START_REF] Philip | Reply To "Comments on Steady Infiltration from Spherical Cavities[END_REF]. This parameter expresses the relative importance of gravity and capillary forces during water movement in unsaturated soils [START_REF] Raats | Analytical Solutions of a Simplified Flow Equation[END_REF]. More specifically, low λc values indicate a dominance of gravity over capillarity, which occurs in coarse-textured and/or highly structured porous media.

On the contrary, high λc values indicate dominance of capillarity over gravity, which occurs in finetextured and/or unstructured porous media (Reynolds et al., 2002).

The λc parameter depends on the soil hydraulic characteristics and the initial condition [START_REF] Smith | Infiltration theory for hydrologic applications[END_REF], and is defined as follows [START_REF] Philip | Reply To "Comments on Steady Infiltration from Spherical Cavities[END_REF]:

λ 𝑐 = ϕ 𝑚 Δ𝐾 (2)
where ΔK (=Ks -Ki) stands for the difference between Ks and the initial soil hydraulic conductivity, Ki (L T -1 ), and ϕm (L 2 T -1 ) is the matric flux potential, defined by [START_REF] Gardner | Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a water table[END_REF]:

ϕ 𝑚 = ∫ 𝐾(ℎ)𝑑ℎ 0 ℎ 𝑖 ℎ 𝑖 ≤ ℎ ≤ 0 (3)
where hi (L) is the initial soil water pressure head.

Eq. ( 1) has been used in numerous studies to solve the flow problem (e.g., [START_REF] Philip | Steady Infiltration From Buried Point Sources and Spherical Cavities[END_REF][START_REF] Wooding | Steady Infiltration from a Shallow Circular Pond[END_REF]. Therefore, the knowledge λc has a noticeable practical interest since these parameters were included in many infiltrometer methods for calculating the saturated soil hydraulic conductivity [START_REF] Bagarello | A Simplified Falling-Head Technique for Rapid Determination of Field-Saturated Hydraulic Conductivity[END_REF](Bagarello et al., , 2014d(Bagarello et al., , e.g., 2017;;[START_REF] Elrick | Methods for analyzing constant-head well permeameter data[END_REF][START_REF] Nimmo | Rapid Measurement of Field-Saturated Hydraulic Conductivity for Areal Characterization[END_REF][START_REF] Reynolds | Ponded Infiltration From a Single Ring: I. Analysis of Steady Flow[END_REF]Stewart and Abou Najm, 2018a;[START_REF] Wu | Measuring Saturated Hydraulic Conductivity using a Generalized Solution for Single-Ring Infiltrometers[END_REF]. Sometimes, its reciprocal, α, is considered, instead.

Many infiltrometer techniques have been developed overtime to estimate λc (e.g., [START_REF] Bagarello | Testing a Simplified Approach to Determine Field Saturated Soil Hydraulic Conductivity[END_REF][START_REF] Reynolds | Ponded Infiltration From a Single Ring: I. Analysis of Steady Flow[END_REF]Stewart and Abou Najm, 2018a;[START_REF] Wu | Measuring Saturated Hydraulic Conductivity using a Generalized Solution for Single-Ring Infiltrometers[END_REF]. However, all these methods present some limitation. For instance, the two-ponding depths method by Reynolds and [START_REF] Reynolds | Ponded Infiltration From a Single Ring: I. Analysis of Steady Flow[END_REF] requires the execution of complicate experiments in the field for measuring the steady state flow rates corresponding to two ponding depths of water. [START_REF] Bagarello | Testing a Simplified Approach to Determine Field Saturated Soil Hydraulic Conductivity[END_REF] proposed some empirical equations, that however need to be calibrated on experimental infiltration data.

Approach 3 by Stewart and Abou Najm (2018b) uses water retention data to constrain λc, thus this method requires the estimation of the θ(h) relationship. Other methods are based on the analysis of transient state data, as for the case of the method 1 by [START_REF] Wu | Measuring Saturated Hydraulic Conductivity using a Generalized Solution for Single-Ring Infiltrometers[END_REF] and Approach 2 by Stewart and Abou Najm (2018a), that need an accurate and often not easy to obtain description of the transient state [START_REF] Di Prima | Experimental assessment of a new comprehensive model for single ring infiltration data[END_REF]. Therefore, alternative methods for estimating λc from simple and easily replicable infiltration experiments could substantially reduce the amount of work necessary to obtain accurate estimates of soil hydraulic properties in the field (Bagarello et al., 2014d).

The objective of this investigation was to validate a simple field method to estimate the macroscopic capillary length by only using a single-ring infiltration experiment of the Beerkan type [START_REF] Lassabatere | Beerkan estimation of soil transfer parameters through infiltration experiments-BEST[END_REF] and a measurement of the initial and saturated soil water contents. The theoretical analysis was initially developed. The proposed method was then validated using analytically generated data and a Beerkan infiltration database including 433 field experiments.

Theory

According to [START_REF] White | Macroscopic and microscopic capillary length and time scales from field infiltration[END_REF], Eq. ( 2) could be rewritten as:

λ 𝑐 = 𝑏 𝑆 2 (∆θ ∆𝐾) (4)
where b is a dimensionless constant dependent on the shape of the soil water diffusivity function, S (L T -0.5 ) is the soil sorptivity [START_REF] Philip | The theory of infiltration: 4. Sorptivity and algebraic infiltration equations[END_REF], Δθ (=θs -θi) stands for the difference between the saturated, θs (L 3 L -3 ), and initial, θi (L 3 L -3 ), volumetric soil water content. For field soils, b is commonly set equal to 0.55 and Ki is negligible [START_REF] White | On the variability and use of the hydraulic conductivity alpha parameter in stochastic treatments of unsaturated flow[END_REF].

The estimation of the macroscopic capillary length with equation (4) requires the previous determination of sorptivity and hydraulic conductivity. These quantities can be estimated thanks to water infiltration experiments and fitting to the quasi-exact implicit (QEI) model developed by [START_REF] Haverkamp | Three-dimensional analysis of infiltration from the disc infiltrometer: 2. Physically based infiltration equation[END_REF] or its related approximate expansions (see Lassabatere et al., 2009 for more details). [START_REF] Haverkamp | Three-dimensional analysis of infiltration from the disc infiltrometer: 2. Physically based infiltration equation[END_REF] proposed the following approximate expansion for the description of the steady-state for water infiltrated through a disc source while maintaining a zero water pressure head at surface:

𝐼 3𝐷 +∞ (𝑡) = (𝐾 𝑠 + γ𝑆 2 𝑟∆θ ) 𝑡 + 𝑆 2 2(1-β)∆𝐾 𝑙𝑛 ( 1 β ) (5)
Where r (L) is the radius of the source, γ and β are two infiltration constants, often fixed at 0.75 and 0.6, respectively [START_REF] Haverkamp | Three-dimensional analysis of infiltration from the disc infiltrometer: 2. Physically based infiltration equation[END_REF]. 𝑏 𝑠 (L) and 𝑖 𝑠 (L T -1 ) that are the intercept and slope of the steady state approximate expansion are defined as a function of hydraulic conductivity and sorptivity as follows:

𝑖 𝑠 = 𝐾 𝑠 + 𝛾𝑆 2 𝑟 𝑑 ∆𝜃 (6a) 𝑏 𝑠 = 𝑆 2 2(1-𝛽)∆𝐾 𝑙𝑛 ( 1 𝛽 ) (6b)
In this study, we use the equation defining the intercept, 𝑏 𝑠 , to quantify the ratio between sorptivity and the difference in hydraulic conductivity, as already suggested by Castellini et al. (2018):

𝑆 2 ∆𝐾 = 𝑏 𝑠 𝐶 (7a) 𝐶 = 1 2(1-β) 𝑙𝑛 ( 1 β ) (7b)
The implementation of such expression into equation (4) leads to the following expression for the macroscopic capillary length:

λ 𝑐 = 𝑏 ∆𝜃 𝑏 𝑠 𝐶 (8)
Eq. ( 8) constitutes a considerable simplification, allowing to estimate the macroscopic capillary length, λc, by only using the infiltration experiment and a measurement of difference between the initial and saturated soil water contents, ∆𝜃. Indeed, once the intercept of the steady state, bs, is estimated, using for instance the linear regression of the last points of the cumulative infiltration, and once the parameter b is determined, λ 𝑐 can be easily determined. Note that [START_REF] White | Macroscopic and microscopic capillary length and time scales from field infiltration[END_REF] fixed b at 0.55 for the case when the initial hydraulic conductivity is far lower than the saturated hydraulic conductivity, Ki << Ks. In addition, β is fixed at 0.6, leading to a value of 0.639.

The previous estimation of the macroscopic capillary length may allow the estimation of the saturated hydraulic conductivity. Many methods for calculating the saturated soil hydraulic conductivity from single-ring data require an estimation of λc or its inverse α. Thus, in this investigation, we also assessed the impact of the proposed approach on Ks calculation. Among the λcdependent methods, we considered the following methods for estimating the saturated hydraulic conductivity:

i) the One-Ponding Depth (OPD) method by [START_REF] Reynolds | Ponded Infiltration From a Single Ring: I. Analysis of Steady Flow[END_REF];

ii) the method 2 by [START_REF] Wu | Measuring Saturated Hydraulic Conductivity using a Generalized Solution for Single-Ring Infiltrometers[END_REF]

Material and methods

Analytically generated data

We assessed the accuracy of the proposed calculation approach to yield reliable λc and Ks predictions by using the same six soils (i.e., sand, loamy sand, sandy loam, loam, silt loam, silty clay loam) considered by [START_REF] Hinnell | Explicit infiltration function for boreholes under constant head conditions: INFILTRATION FUNCTION FOR BOREHOLES[END_REF] and [START_REF] Bagarello | Estimating saturated soil hydraulic conductivity by the near steady-state phase of a Beerkan infiltration test[END_REF] to cover a wide range of hydraulic responses. We modelled the infiltration experiments for these synthetic soils using the parameters listed by [START_REF] Carsel | Developing joint probability distributions of soil water retention characteristics[END_REF] to describe the water retention curve and the hydraulic conductivity function according to the van Genuchten-Mualem model [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF][START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF]) (Table C):

{ θ(ℎ) = θ 𝑟 + (θ 𝑠 -θ 𝑟 )[1 + α 𝑉𝐺 |ℎ| 𝑛 ] -𝑚 ℎ < 0 θ(ℎ) = θ 𝑠 ℎ ≥ 0 (9a) 𝑚 = 1 - 1 𝑛 (9b) 𝐾(ℎ) = 𝐾 ̂𝑠𝑆𝑒 𝑙 [1 -(1 -𝑆𝑒 1 𝑚 ⁄ ) 𝑚 ] 2 (9c)
where 𝐾 ̂𝑠 (L T -1 ) is the actual value of the saturated soil hydraulic conductivity, αVG is an empirical parameter related to the inverse of the air-entry pressure head, n is the pore size distribution index, Se = (θ -θr)/(θs -θr) is the effective saturation degree, θr (L 3 L -3 ) is the residual water content, l is the pore connectivity parameter, assumed to be 0.5 by [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF]. To also test the effect of the initial soil water content on parameters predictions, initial values of Se ranging from 0.1 to 0.8 were considered when modeling infiltration. For a given soil and an initial Se value, the actual value of the macroscopic capillary length, λ ̂𝑐, was estimated by Eq. (2).

The cumulative infiltration curves were generated using the infiltration model proposed by [START_REF] Haverkamp | Three-dimensional analysis of infiltration from the disc infiltrometer: 2. Physically based infiltration equation[END_REF]. More specifically, for any zero or negative value of the pressure head at the soil surface, 1D cumulative infiltration, I1D (L), into an initially uniform unsaturated soil profile can be modelled by the following 1D implicit equation [START_REF] Haverkamp | Infiltration under ponded conditions: 3. A predictive equation based on physical parameters[END_REF]:

2Δ𝐾 2 𝑆 2 𝑡 = 1 1-β [ 2Δ𝐾 𝑆 2 (𝐼 1𝐷 (𝑡) -𝐾 𝑖 𝑡) -𝑙𝑛 ( 𝑒𝑥𝑝(2β Δ𝐾 𝑆 2 (𝐼 1𝐷 (𝑡)-𝐾 𝑖 𝑡))+β-1 β )] (10)
with ΔK = 𝐾 ̂𝑠 -Ki. Eq. ( 10) was extended to three-dimensional (3D) infiltration from a surface disk source by [START_REF] Smettem | Three-dimensional analysis of infiltration from the disc infiltrometer: 1. A capillary-based theory[END_REF], adding a term representing 3D geometrical effects [see Eq. ( 23)

of [START_REF] Smettem | Three-dimensional analysis of infiltration from the disc infiltrometer: 1. A capillary-based theory[END_REF]:

𝐼(𝑡) = 𝐼 1𝐷 (𝑡) + 𝛾𝑆 2 𝑟 𝑑 Δθ ( 11
)
where I is 3D cumulative infiltration, γ is a shape parameter for geometrical correction of the infiltration front shape, commonly set at 0.75 for θi < 0.25 θs [START_REF] Haverkamp | Three-dimensional analysis of infiltration from the disc infiltrometer: 2. Physically based infiltration equation[END_REF]. For a given soil and an initial Se value, sorptivity was estimated from soil hydraulic properties in the interval of soil water pressure head between hi and 0 as follows [START_REF] Parlange | Convergence and validity of time expansion solutions: a comparison to exact and approximate solutions[END_REF]:

𝑆 = √∫ (θ 𝑠 + θ -2θ 𝑖 )𝐾(ℎ)𝑑ℎ 0 ℎ 𝑖 (12)
Then, Eqs. ( 10) and ( 11) were used to generate cumulative infiltration curves. Shape parameters β and γ were assumed equal to the default values, i.e., 0.6 and 0.75, respectively (Angulo-Jaramillo et al., 2019). The duration of the infiltration process was fixed at three times the maximum time, tmax (T), for which the transient expression of the 3D Haverkamp equation (see Eq. (3a) in [START_REF] Lassabatere | Beerkan estimation of soil transfer parameters through infiltration experiments-BEST[END_REF]) can be considered valid, with tmax calculated as follows [START_REF] Lassabatere | Beerkan estimation of soil transfer parameters through infiltration experiments-BEST[END_REF]:

𝑡 𝑚𝑎𝑥 = 1 4(1-𝐵) 2 ( 𝑆 𝐾 ̂𝑠) 2 (13a) 𝐵 = 2-𝛽 3 (1 - 𝐾 𝑖 𝐾 ̂𝑠) + 𝐾 𝑖 𝐾 ̂𝑠 (13b)
These analytical data were used to estimate the intercept, bs, and the slope, is, by linear regression analysis of the last three data points of the cumulative infiltration versus time plot, which were assumed to provide a reliable description of steady state conditions. Then, we estimated λc and Ks by Eqs. ( 8) and (t1-4), respectively. These estimates were compared with the corresponding actual values, i.e., the λ ̂𝑐 values estimated with Eq. ( 2) and the 𝐾 ̂𝑠 values by [START_REF] Carsel | Developing joint probability distributions of soil water retention characteristics[END_REF] listed in Table C, using the relative error, Er, defined as follows:

𝐸𝑟(𝑥) = 100 × 𝑥-𝑥 x ̂ (14)
Where 𝑥 ̂ is the actual value for λ ̂𝑐 or 𝐾 ̂𝑠 and x is the corresponding estimated value λc or Ks. Note that positive Er(x) values indicate overestimations whereas negative values indicate underestimations. Small deviations, i.e., Er(x) ~ 0, suggest that the estimated soil hydraulic parameters are close to actual values. According to the accuracy criterion by [START_REF] Reynolds | An assessment of borehole infiltration analyses for measuring field-saturated hydraulic conductivity in the vadose zone[END_REF], the estimates were deemed accurate when they fell within the range 0.75 < x/ 𝑥 ̂ < 1.25 (i.e. < 25% error).

Such a stringent accuracy criterion was used because the parameters were estimated by analytically generated data, and therefore they were free of the perturbations embedded in field and laboratory measurements (e.g., measurement error, random noise and natural variability).

The Beerkan infiltration database

In this investigation we considered a large database of Beerkan infiltration experiments carried out in four different countries, i.e., Italy, Burundi, France and Spain, during the period 2010-2017 (Table S). However, nearly half of the runs, i.e., the 46.7% were carried out in Sicily (Italy) (202 out of 433).

While, the 35,1% of the runs (152 out of 433) were carried out in Burundi in the African great lakes region. The runs were carried out at 433 sampling points on soils having contrasting textural characteristics, i.e., from sandy to clayey (Figure T).

All Beerkan experiments were carried out according to the methodology described by [START_REF] Lassabatere | Beerkan estimation of soil transfer parameters through infiltration experiments-BEST[END_REF]. As prescribed by the Beerkan method, a stainless steel ring was inserted shallowly into the soil (~1 cm). For each run, different volumes of water were successively poured on the confined soil surface. The number of poured volumes varied depending on time needed to reach steady state, as required by the Beerkan method [START_REF] Angulo-Jaramillo | Beerkan Estimation of Soil Transfer parameters (BEST) across soils and scales[END_REF]. The Beerkan experiments considered in this investigation were performed using rings with inner diameters from 5 to 30 cm.

The volume of water for each pouring was chosen in order to establish ponding condition on the entire infiltration surface; the volumes ranged from 17 to 800 mL, depending on the used ring (Table S).

The energy of the falling water was dissipated with fingers to minimize the soil disturbance owing to water pouring, as commonly suggested (e.g., Di Prima et al., 2019). For each poured volume, the time needed for the water to infiltrate was recorded, and the cumulative infiltration, I (mm), was plotted against time, t (h).

Similarly to the procedure previously described to analyse analytical generated curves, for each Beerkan experiment we estimated the intercept, bs (mm), and the slope, is (mm h -1 ), of the regression line fitted to the data points describing steady state conditions. Then, we estimated λc and Ks by Eqs.

(8) and (t1-4), respectively.

For the field data, we chose not to use a reference λc value, as an independent datum that can be used for assessing simplified procedures or validating new developed methods does not currently exist. In addition, methods for estimating λc from single-ring infiltration experiments are based on cumbersome field procedures.

According to [START_REF] Elrick | Methods for analyzing constant-head well permeameter data[END_REF], we classified the λc estimates in five soil capillarity categories, i.e., from negligible (λc ≤ 10 mm) to very strong (λc ≥ 1000 mm) (Table L). Note that these categories were originally proposed to select a representative value of the sorptive number, α (L -1 ), which is equal to λ -1 , on the basis of five soil texture-structure categories when calculating Ks by the OPD method [START_REF] Angulo-Jaramillo | Infiltration Measurements for Soil Hydraulic Characterization[END_REF]. In this investigation, we also propose range values for each category as detailed in Table L.

For Ks, we used as a benchmark the values obtained with the BEST-steady method proposed by Bagarello et al. (2014c), which estimates Ks as follows:

𝐾 𝑠,𝐵𝐸𝑆𝑇 = 𝐶 𝑖 𝑠 𝐴 𝑏 𝑠 +𝐶 (15a) 𝐴 = 𝛾 𝑟(𝜃 𝑠 -𝜃 𝑖 ) (15b)
Thus, the Ks values estimated by Eqs. (t1-4) were compared with those obtained by Eq. ( 15) in terms of relative errors using Eq. ( 14). We chose this method, as it involves the same experimental information of the λc-dependent methods considered in this investigation for estimating Ks, but it does not require an estimation of λc.

We avoided to use laboratory measurements as benchmark, as they can induce experimental artifacts [START_REF] Haverkamp | Soil properties and moisture movement in the unsaturated zone[END_REF], such as soil compaction and samples biased by pores, that may limit their comparability with in-situ measurements (e.g., [START_REF] Bagarello | Testing laboratory methods to determine the anisotropy of saturated hydraulic conductivity in a sandy-loam soil[END_REF]Di Prima et al., 2018c). Discrepancies are also expected when different measurement techniques are applied in the field or even when the same dataset is analyzed by alternative calculation approaches [START_REF] Mertens | Characterisation of the field-saturated hydraulic conductivity on a hillslope: in situ single ring pressure infiltrometer measurements[END_REF], though in the latter case the results can still be compared to one another [START_REF] Wu | Measuring Saturated Hydraulic Conductivity using a Generalized Solution for Single-Ring Infiltrometers[END_REF].

For comparisons between paired observations, the paired differences, i.e., Ks -Ks,BEST for given method, were calculated and the hypothesis of normality of these differences was checked by the Kolmogorov-Smirnov test. For non-normally distributed data the Wilcoxon signed rank test was used to test the median difference between paired observations at P < 0.05.

Results

Analytical validation

Actual values of the macroscopic capillary length

As expected, Eq. ( 2) yielded higher λ ̂𝑐 values for fine-textured soils and lower Se values (Table P). This is logical, since for fine soils the capillary contribution to water flow was higher than for coarser soils. More specifically, high λc values are associated to initially flat K(h) relationships, i.e., when a decrease in pressure head determines a moderate pore emptying [START_REF] Angulo-Jaramillo | Infiltration Measurements for Soil Hydraulic Characterization[END_REF][START_REF] Reynolds | Letter to the editor. Comment on "comparison of three methods for assessing soil hydraulic properties" by G. B. Paige and D. Hillel[END_REF]. As an example, this behavior can be observed on the hydraulic conductivity

Estimating λc from analytically generated data

Cumulative infiltration could be calculated for all cases trough Eqs. ( 10) and ( 11) and using the parameters listed in Table C. All curves exhibited usual shapes, with a concave shape as a function of time, due to a decreasing hydraulic gradient as the wetting front moves away from the source and the influence of capillarity decreases [START_REF] Xu | Analysis of single-ring infiltrometer data for soil hydraulic properties estimation: Comparison of BEST and Wu methods[END_REF]. As the process approaches steady state, cumulative infiltration curves become approximately linear with time (Angulo-Jaramillo et al., 2016)

(Figure H). Note that the duration of the infiltration process decreased for higher values of the initial effective saturation degree. Indeed, for wetter initial conditions steady state was reached more quickly and with less water [START_REF] Di Prima | Testing a new automated single ring infiltrometer for Beerkan infiltration experiments[END_REF]. In this case, the influence of capillary forces on infiltration was very short, and the cumulative infiltration curves exhibited a more linear shape [START_REF] Angulo-Jaramillo | Beerkan Estimation of Soil Transfer parameters (BEST) across soils and scales[END_REF].

The values of the intercept, bs, estimated from these curves, were used in conjunction with the θi and θs values to calculate λc by Eq. ( 8) (Table P) G). Under these conditions, the ratio Ki / Ks departed from the assumption made to approximate Eq. (7a) (i.e., Ki ≪ Ks). More specifically, near the saturation, the ratio Ki / Ks was higher for the coarser soils than for the fine-textured soils because the capillary forces were weaker and the increase of the hydraulic conductivity near the saturation was less abrupt (Figure K). As a consequence, for the coarser soils, neglecting the Ki / Ks ratio in Eq. (7a) negatively affected λc estimates, yielding higher errors. Nevertheless, λc estimates were sufficiently accurate also in these cases, with errors always lower than 25%.

Estimating Ks from analytically generated data

The values of the slope, is, estimated from the analytically generated curve, were used to calculate

Ks by the four λc-dependent methods (i.e., OPD, WU2, SSBI, A4) using Eqs. (t1-4) (

Field testing

Peculiarities and experiments exclusion

For six Beerkan experiments infiltration rates increased with time, such that the cumulative infiltration curves exhibited convex shapes. The cumulative infiltration data in that instance yielded a negative value for the intercept, bs. This circumstance led to negative λc data. In addition, in case of negative intercept values, the BEST-steady algorithm is also expected to estimate meaningless Ks data [START_REF] Angulo-Jaramillo | Beerkan Estimation of Soil Transfer parameters (BEST) across soils and scales[END_REF]. For this reason these Beerkan experiments were excluded from subsequent analyses. These cases occurred at the Kinyami site (two instances), Palermo -SAAF site (one instance) and Crépieux-Charmy site (three instances), and reflect that the early wetting phase was impeded due to soil water repellency [START_REF] Concialdi | An open-source instrumentation package for intensive soil hydraulic characterization[END_REF]. Hydrophobia may be attributed to locally high OC content [START_REF] Goebel | Soil water repellency and its implications for organic matter decomposition -is there a link to extreme climatic events?[END_REF], exudates produced by the plant root systems or living organisms like arbuscular mycorrhizal fungi [START_REF] Rillig | Mycelium of arbuscular mycorrhizal fungi increases soil water repellency and is sufficient to maintain water-stable soil aggregates[END_REF], and organic pollutants such as hydrocarbons [START_REF] Durand | Characterisation of complex organic matter present in contaminated sediments from water retention ponds[END_REF].

Estimating λc from the Beerkan infiltration database

The Beerkan experiments that were not affected by flow impedance (427 out of 433) could be analyzed to determine λc. Eq. ( 8) yielded λc values ranging from 1.5 to 737.7 mm (Table Q).

According to the soil capillarity categories suggested by [START_REF] Elrick | Methods for analyzing constant-head well permeameter data[END_REF], the λc values were classified as representative of strong capillarity conditions in 29.7% of the cases (127 out of 427), moderate in the 44.3% (189), weak in the 25.1% (107), and negligible in the 0.9% (4) (Figure D). The categories weak and negligible were associated with the higher sand contents, i.e., 35% and 37%, respectively. Thus, for the field experiments, Eq. ( 3) estimated the lowest λc values for the coarsest soils, i.e., when gravity was expected to prevail over capillarity. This result increases our confidence on the proposed calculation approach. 2019), these cases correspond to three typical responses when applying the Beerkan protocol respectively on i) regular soils, ii) very fine soils and iii) very coarse soils (see their Figure 13).

In the first case (Figure Fb), cumulative infiltration exhibits an usual concave shape as a function of time. For this run, we estimated an intercept value of 30.9 mm, and a λc value of 49.2 mm falling within the range of 42 ≤ λc ≤ 125 mm (Figure Fa). According to the capillarity categories suggested by [START_REF] Elrick | Methods for analyzing constant-head well permeameter data[END_REF], the estimated λc value is representative of moderate capillarity conditions (Table L).

In the second case (Figure Fc), the cumulative infiltration curve exhibits a strong concave shape with an intercept value of 209.6 mm, yielding a λc value of 598.7 mm, falling within the strong capillary category (Table L andFigure Fa). This behavior is typical for very fine soils, with low permeability. For this run, capillary forces predominates for almost the duration of the experiment and the attainment of steady state required more time [START_REF] Di Prima | Testing a new automated single ring infiltrometer for Beerkan infiltration experiments[END_REF].

In the third case (Figure Fd), the cumulative infiltration curve exhibits an almost linear shape with an intercept value of 1.0 mm, yielding a λc value of 1.5 mm, i.e., lower than the considered threshold of 10 mm for negligible capillarity forces (Table L andFigure Fa). This behaviour is typical for coarse-textured soils and occurs when the infiltration process is mainly driven by gravity. Negligible values of λc were also reported in other investigations (e.g., [START_REF] Khaleel | Variability of Gardner's α for coarse-textured sediments[END_REF][START_REF] Russo | On the spatial variability of parameters of the unsaturated hydraulic conductivity[END_REF][START_REF] White | On the variability and use of the hydraulic conductivity alpha parameter in stochastic treatments of unsaturated flow[END_REF]. In this investigation we obtained negligible λc values in only four cases out of 427. These cases corresponded to the lowest measured intercept values (Figure Fa).

Estimating Ks from the Beerkan infiltration database

The differences between Ks and Ks,BEST were always non-normally distributed according to the Kolmogorov-Smirnov test. The Wilcoxon signed rank test showed that the four methods (i.e., OPD, WU2, SSBI, A4) yielded Ks estimates significantly different from the BEST-steady values. However, the discrepancies between methods were always lower than < 25%, with the exception of the SSBI method (Figure Z). More specifically, the method that yielded the lowest discrepancies was the WU2 method. This method yielded Ks values ranging from 1.4 to 3493.8 mm h -1 (Table Q), with 76 and 24% of the runs yielding respectively lower and higher Ks estimates than the BEST-estimated values (Figure Z). The median Ks values differed by only a factor of 1.002. The relative error, Er(Ks), ranged from -12.1 to 22.1%. Therefore, according to selected criterion (i.e., discrepancies < 25%), the WU2 and BEST-steady methods yielded similar Ks estimates. The OPD and A4 methods also gave |Er(Ks)| values always lower than 25%, suggesting that also these methods yielded similar Ks estimates to those values obtained by BEST-steady. Larger discrepancies were only obtained when using the SSBI method. This method always yielded Ks values higher than Ks,BEST, with only one exception. The Er(Ks) values ranged from -1.0 to 32.8%, with 21% of the runs (88 out of 427) yielding higher values than the considered threshold of 25%.

Discussion

In this investigation, we proposed to estimate the macroscopic capillary length using a simple

Beerkan infiltration experiment and a measurement of the initial and saturated soil water contents.

Previous investigations also suggested that the measured infiltration curve contains the necessary information to estimate λc [START_REF] Bagarello | Testing a Simplified Approach to Determine Field Saturated Soil Hydraulic Conductivity[END_REF](Bagarello et al., , e.g., 2014d;;Stewart and Abou Najm, 2018a;[START_REF] Wu | Measuring Saturated Hydraulic Conductivity using a Generalized Solution for Single-Ring Infiltrometers[END_REF]. However, these methods are based on the analysis of the transient infiltration process.

Indeed, they all require an estimation of the infiltration coefficients by the fitting of transient models to early infiltration data. It is not uncommon that these methods perform poorly, especially in case of uncertainties on the transient phase (e.g., Di Prima et al., 2019[START_REF] Di Prima | Comparing Transient and Steady-State Analysis of Single-Ring Infiltrometer Data for an Abandoned Field Affected by Fire in Eastern Spain[END_REF]. While these methods may alleviate the experimental efforts needed to determine λc and Ks in the field, for instance, reducing the required measurement time for low permeability soils (Bagarello et al., 2014d), they may introduce uncertainties on the hydraulic characterization and increase parameter variability [START_REF] Di Prima | Comparing Transient and Steady-State Analysis of Single-Ring Infiltrometer Data for an Abandoned Field Affected by Fire in Eastern Spain[END_REF]. This limitation can ultimately complicate an accurate assessment of spatial-temporal variability of soil hydraulic properties [START_REF] Angulo-Jaramillo | Infiltration Measurements for Soil Hydraulic Characterization[END_REF].

On the other hand, the proposed method uses the experimental information collected during the final stage of the infiltration process, thus, it is expected to be less sensitive to experimental data accuracy [START_REF] Bagarello | Testing a Simplified Approach to Determine Field Saturated Soil Hydraulic Conductivity[END_REF]. Estimating λc requires the intercept (bs) that can be easily calculated by linear regression analysis of the final portion of the cumulative infiltration. While the final infiltration rate (is) is totally independent from the early phase of the infiltration process, the intercept depends on the shape of the cumulative infiltration curve, thus, it restitutes a certain information on the transient phase (Bagarello et al., 2014c). The shape of the cumulative infiltration, in turn, reveals the relative importance of capillary and gravity forces during a ponding infiltration process, being more concave or linear respectively in case of prevalence of capillarity over gravity or vice versa [START_REF] Angulo-Jaramillo | Beerkan Estimation of Soil Transfer parameters (BEST) across soils and scales[END_REF]. More specifically, small values of the intercept are indicative of a linear infiltration curve, i.e., when gravity prevails over capillarity, which occurs primarily in coarsetextured and/or highly structured porous media. On the contrary, high intercept values are indicative of a concave infiltration curve, i.e., when capillarity prevails over gravity, which occurs primarily in fine-textured soils. Therefore, bs is expected to be a reliable predictor of the macroscopic capillary length.

In this investigation, this hypothesis was validated by assessing the proposed calculation approach using both analytical and field data. The analytical validation demonstrated that Eq. ( 8) was able to provide reliable λc estimates over a wide range of soil textural characteristics and initial soil water contents. We also tested the new method on a large database including 433 Beerkan infiltration experiments carried out on different soils having contrasting textural characteristics, i.e., from sandy to clayey. The 99.1% of the experiments yielded λc values falling within the realistic range 10 ≤ λc ≤ 1000 mm (W. D. Reynolds and Elrick, 2002), with the only expectation of four cases yielding negligible λc values, i.e., lower than 10 mm.

An expected advantage of the proposed method is to improve Ks prediction since it avoids to choose a λc value on the basis of the soil capillarity categories. Indeed, using Eqs. (t1-4) needs knowledge of λc, that, according to [START_REF] Elrick | Methods for analyzing constant-head well permeameter data[END_REF], can be estimated on the basis of a general description of soil textural and structural characteristics (Table L). The effect of an erroneous choice of the capillarity category on Ks predictions was investigated by Bagarello et al. (2014d). These authors assessed the sensitivity of Ks to the right choice of the capillarity category, reporting that an erroneous choice of the capillarity category, for instance differing by plus or minus one category from the correct one, may also lead to discrepancies between Ks estimates of more than a factor of three.

A sensitivity analysis to determine the effect of errors in constraining λc on Ks estimations was also performed by Stewart and Abou Najm (2018b). These authors reported that Ks estimates are more sensitive to underpredictions of λc compared with overpredictions. Therefore, choosing an erroneous λc value from the previous capillarity category may have a more relevant effect than choosing the λc value from the successive category. In any case, the subjectivity on the choice of an appropriate λc value constitutes a clear limitation when using λc-dependent methods for estimating Ks, such as those considered in this investigation. On the contrary, the proposed method avoids any subjectivity on the estimation of λc.

A possible limitation when using this calculation approach, especially on very fine soils, is related to the need of infiltrating enough water during the experiment in order to approach steady state conditions. Indeed, fine-textured and unstructured soils often exhibit very low permeability, thus, the attainment of steady state may require excessively long experiments [START_REF] Di Prima | Testing a new automated single ring infiltrometer for Beerkan infiltration experiments[END_REF], which would complicate reaching steady state conditions with the regular Beerkan protocol (Alagna et al., 2016a). This problem can be easily solved using automated devises, such as the one developed by Di Prima (2015) and [START_REF] Concialdi | An open-source instrumentation package for intensive soil hydraulic characterization[END_REF], that allows accurate and automated infiltration rates measurements. The first obvious advantage in using automatic infiltrometer in case of low permeability soils is that data collection is automated and it does not require any further manipulation by the user once the device is installed. Therefore, tests can be simultaneously performed alleviating the experimental efforts in case of long experiments durations [START_REF] Alagna | Alternative analysis of transient infiltration experiment to estimate soil water repellency[END_REF].

Summary and conclusions

In this investigation, we assessed a simple field method for estimating the macroscopic capillary length, λc, by only using a single-ring infiltration experiment of the Beerkan type and a measurement of the initial and saturated soil water contents. We validated the proposed method by using both analytical generated data and a large database of 433 Beerkan infiltration experiments carried out in four different countries, i.e., Italy, Burundi, France and Spain, during the period 2010-2017. While the analytical validation supported our hypothesis that the intercept, bs, is a reliable predictor of the macroscopic capillary length, the testing carried out using the Beerkan database increased our confidence on the proposed calculation approach. This method constitutes an easy solution to estimate λc, and it can be applied in conjunction with a number of λc-dependent methods for estimating the saturated soil hydraulic conductivity, allowing to improve Ks prediction in the field. The proposed method may also avoid any subjectivity caused by the selection of a representative λc value on the only base of textural and structural characteristics. is (L T -1 ) = steady-state infiltration rate estimated from the cumulative infiltration versus time plot. λc (L) = macroscopic capillary length. r (L) = ring radius. γ = dimensionless constant set equal to 0.75 and related to the effect of disk curvature and gravity [START_REF] Smettem | Three-dimensional analysis of infiltration from the disc infiltrometer: 1. A capillary-based theory[END_REF]. γw = dimensionless constant related to the shape of the wetting front [START_REF] White | Macroscopic and microscopic capillary length and time scales from field infiltration[END_REF]. A value of γw = 1.818 was suggested by Reynolds and Elrick (2002).

  function of the silty clay loam soil (Figure K), where the functional relationship between K and Se is very flat for low Se values, i.e. high pressure head. On the other hand, low λc values are associated to steep relationships, i.e., when a small decrease in pressure head determines a more relevant pore emptying (e.g., sand, Figure K).

.

  Then, the λc values were compared with the corresponding actual values, λ ̂𝑐, in terms of relative errors (Figure G). All λc values were accurate according to the chosen criterion, since Er(λc) ranged from -23.2 to 7.9%. The higher |Er(λc)| values were obtained for the coarse-textured soils, namely the sandy and loamy sand soils, and only for initial wet conditions (Figure

  Figure Fd). According to Angulo-Jaramillo et al. (2019), these cases correspond to three typical

  = ring insertion depth into the soil. 𝐺 * = d + r/2. H (L) = ponding depth of water.

Figure H .

 H Figure H. Cumulative infiltration curves for different soils and initial effective saturation degrees, Se. The

Figure G .

 G Figure G. Values of the relative error of the macroscopic capillary length, Er(λc), predictions (expressed as a percentage of the actual values, 𝛌 ̂𝒄) estimated for the six synthetic soils listed in Table C.

Figure W .

 W Figure W. Values of the relative error of the saturated soil hydraulic conductivity, Er(Ks), predictions (expressed as a percentage of the actual values, 𝑲 ̂𝒔) estimated for the six synthetic soils by the four

Figure

  Figure F. (a) Comparison between the estimated macroscopic capillary length, λc, and the intercept, bs,

Figure D .

 D Figure D. Sample size, N, of the macroscopic capillary length predictions and mean sand content (%) for

Figure Z .

 Z Figure Z. Cumulative empirical frequency distribution of the saturated soil hydraulic conductivity, Er(Ks),

Table P )

 P . Then, the Ks values were compared with the corresponding actual values, 𝐾 ̂𝑠 (TableC), in terms of relative errors (FigureW). The Er(Ks) values ranged from -3.1 to 9.5%, from -7.4 to 5.4%, from -24.7 to -2.2%, and from 2.4 to 14%, for the OPD, WU2, SSBI and A4 methods, respectively. Mean |Er(Ks)| values were ordered as OPD < A4 < WU2 < SSBI. Therefore, the OPD method was the most accurate overall. However, the four methods always yielded accurate estimates, since |Er(Ks)| values were always lower than 25% (FigureW).

Table E .

 E Relationships for calculating the saturated soil hydraulic conductivity, Ks, by the One-Ponding Depth (OPD) method, the method 2 byWu et al. (1999) (WU2), the Steady version of the Simplified method based on a Beerkan Infiltration run (SSBI), and the Approach 4 (A4).

				Relationships
	OPD		WU2	SSBI	A4
	𝐾 𝑠 =	𝑖 𝑠 𝜋𝑟 2 𝜆 𝑐 𝜆 𝑐 𝑟( 𝐻 +1)+(0.316 (0.316 𝑑 𝑟 +0.184) 𝑑 𝑟 +0.184) 𝜋𝑟 2 𝜆 𝑐	𝐾 𝑠 =	𝑖 𝑠 0.9084( 𝐻+λ 𝑐 𝐺 * +1)
	(t1)

Table C .

 C Soil hydraulic parameters for the five studied soils used to model the infiltration experiments. From[START_REF] Carsel | Developing joint probability distributions of soil water retention characteristics[END_REF].

	Soil texture	Sand	Loamy Sand	Sandy Loam	Loam	Silt Loam	Silty Clay Loam
	θr	0.045	0.057	0.065	0.078	0.067	0.089
	θs	0.43	0.41	0.41	0.43	0.45	0.43
	αVG (mm -1 )	0.0145	0.0124	0.0075	0.0036	0.002	0.001
	n	2.68	2.28	1.89	1.56	1.41	1.23
	𝐾 ̂𝑠 (mm h -1 )	297.0	145.9	44.2	10.44	4.5	0.7
	l	0.5	0.5	0.5	0.5	0.5	0.5

Table L .

 L Soil capillarity categories suggested by[START_REF] Elrick | Methods for analyzing constant-head well permeameter data[END_REF], and representative α (mm -1 ) and λc (mm) values. The suggested range values for the two parameters are also reported.

	Soil capillarity category	Representative α (mm -1 )	Suggested α range values	Representative λc (mm)	λc range values
	Very strong	≤0.001		≥1000
	Strong	0.004	0.001<α<0.008	250	125<λc<1000
	Moderate	0.012	0.008≤α≤0.024	83	42≤λc≤125
	Weak	0.036	0.024<α<0.1	28	10<λc<42
	Negligible	≥0.1		≤10

Table S .

 S Summary of the Beerkan infiltration database. Total number of Beerkan infiltration

	607
	experiments, Ntot = 433.

Table P .

 P Summary of the soil hydraulic properties for the six synthetic soils.

	Soil	Se	hi	θi	S	ϕ	λ ̂𝑐	λc	is	bs		Ks (mm h -1 )
			(mm)	(m 3 m -3 ) (mm h -0.5 ) (mm 2 h -1 )	(mm)	(mm) (mm h -1 ) (mm) OPD	WU2	SSBI	A4
					(Eq. 12)	(Eq. 3)	(Eq. 2)	(Eq. 8)			(Eq. t1) (Eq. t2) (Eq. t3) (Eq. t4)
	Sand	0.1	269.0	0.083	86.5	11308.6	38.1	37.2	516.7 15.0	304.3	319.1	341.3	289.8
		0.2	174.5	0.122	81.3	11298.6	38.0	37.0	515.3 13.2	304.2	319.0	341.0	289.7
		0.3	133.1	0.160	75.7	11268.5	37.9	36.8	513.3 11.5	303.7	318.5	340.3	289.3
		0.4	107.8	0.199	69.6	11202.4	37.7	36.7	510.6	9.8	302.4	317.2	338.8	288.1
		0.5	89.7	0.237	62.9	11076.7	37.3	36.7	506.5	8.2	299.8	314.4	335.9	285.6
		0.6	75.2	0.276	55.5	10854.0	36.5	37.1	500.4	6.6	294.9	309.2	330.7	280.9
		0.7	62.5	0.314	46.9	10467.3	35.2	38.2	491.1	5.1	286.1	299.8	321.6	272.3
		0.8	50.2	0.353	36.8	9776.1	32.9	40.5	476.0	3.6	270.3	282.7	305.4	256.8
	Loamy Sand	0.1	483.8	0.092	58.2	5602.2	38.4	37.4	254.5 13.8	149.5	156.7	167.8	142.4
		0.2	276.4	0.128	54.7	5598.0	38.4	37.2	253.8 12.2	149.5	156.7	167.6	142.4
		0.3	195.6	0.163	50.9	5584.2	38.3	36.9	252.8 10.6	149.3	156.6	167.4	142.3
		0.4	150.0	0.198	46.8	5552.5	38.1	36.6	251.4	9.0	149.0	156.2	166.9	141.9
		0.5	119.2	0.233	42.3	5490.2	37.6	36.4	249.3	7.5	148.1	155.4	165.8	141.2
		0.6	95.9	0.269	37.3	5376.9	36.9	36.3	246.1	6.0	146.4	153.6	163.9	139.5
		0.7	76.5	0.304	31.5	5176.0	35.5	36.5	241.3	4.5	143.2	150.2	160.4	136.5
		0.8	58.9	0.339	24.6	4811.2	33.0	37.4	233.4	3.1	137.2	143.8	153.9	130.6
	Sandy Loam	0.1	1765.2	0.100	36.0	2145.5	49.7	48.3	86.5 17.4	45.4	47.2	52.1	42.9
		0.2	799.2	0.134	33.8	2145.5	49.7	48.0	86.3 15.4	45.4	47.2	52.0	42.9
		0.3	494.2	0.169	31.5	2140.3	49.5	47.6	85.9 13.4	45.4	47.2	52.0	42.9
		0.4	344.1	0.203	29.0	2129.8	49.3	47.1	85.3 11.3	45.3	47.2	51.8	42.9
		0.5	253.1	0.238	26.2	2107.9	48.8	46.5	84.5	9.3	45.2	47.0	51.6	42.7
		0.6	190.3	0.272	23.1	2066.0	47.8	45.7	83.3	7.3	44.8	46.7	51.2	42.5
		0.7	142.4	0.307	19.5	1988.8	46.0	44.9	81.4	5.4	44.2	46.1	50.3	41.8
		0.8	102.3	0.341	15.2	1843.8	42.7	44.0	78.1	3.5	42.8	44.6	48.7	40.6
	Loam	0.1 16941.8	0.113	20.9	719.7	69.2	67.6	24.3 24.9	10.7	11.0	12.6	10.0
		0.2	4883.1	0.148	19.6	719.6	69.2	67.2	24.2 22.0	10.7	11.0	12.6	10.0
		0.3	2330.8	0.184	18.3	718.8	69.1	66.6	24.1 19.1	10.7	11.0	12.6	10.0
		0.4	1354.4	0.219	16.8	716.5	68.9	65.9	23.9 16.2	10.7	11.0	12.6	10.0
		0.5	866.2	0.254	15.2	711.0	68.4	64.8	23.7 13.2	10.7	11.0	12.6	10.0
		0.6	579.6	0.289	13.4	699.3	67.2	63.2	23.3 10.3	10.7	11.0	12.5	10.0
		0.7	390.5	0.324	11.3	676.0	65.0	61.0	22.7	7.5	10.6	10.9	12.4	9.9
		0.8	252.5	0.360	8.8	629.2	60.5	57.7	21.6	4.7	10.4	10.8	12.1	9.8
	Silt Loam	0.1	1.4×10 5	0.105	16.3	402.8	89.5	87.8	12.3 35.1	4.6	4.7	5.6	4.3
		0.2 25267.3	0.144	15.3	402.8	89.5	87.3	12.2 31.1	4.6	4.7	5.6	4.3
		0.3	9318.5	0.182	14.3	402.6	89.5	86.6	12.2 27.0	4.6	4.7	5.6	4.3
		0.4	4529.8	0.220	13.2	401.8	89.3	85.7	12.1 22.9	4.6	4.7	5.6	4.3
		0.5	2531.6	0.258	11.9	399.6	88.8	84.4	12.0 18.8	4.6	4.7	5.6	4.3
		0.6	1519.5	0.297	10.5	394.4	87.6	82.4	11.8 14.7	4.6	4.7	5.6	4.3
		0.7	933.0	0.335	8.9	383.2	85.2	79.3	11.5 10.6	4.6	4.7	5.5	4.3
		0.8	553.5	0.373	7.0	359.1	79.8	74.2	10.9	6.6	4.6	4.7	5.4	4.3
	Silty Clay Loam 0.1	2.2×10 7	0.123	6.0	60.0	85.7	84.9	1.9 30.3	0.7	0.7	0.9	0.7
		0.2	1.1×10 6	0.157	5.6	60.0	85.7	84.5	1.9 26.8	0.7	0.7	0.9	0.7
		0.3 1.9E×10 5	0.191	5.2	60.0	85.7	84.1	1.9 23.3	0.7	0.7	0.9	0.7
		0.4 53399.3	0.225	4.8	60.0	85.7	83.4	1.9 19.8	0.7	0.7	0.9	0.7
		0.5 19954.8	0.259	4.4	59.9	85.5	82.5	1.8 16.3	0.7	0.7	0.9	0.7
		0.6	8725.5	0.294	3.9	59.5	85.0	81.0	1.8 12.8	0.7	0.7	0.9	0.7
		0.7	4137.5	0.328	3.3	58.5	83.6	78.5	1.8	9.3	0.7	0.7	0.9	0.7
		0.8	1966.9	0.362	2.6	55.9	79.9	73.6	1.7	5.8	0.7	0.7	0.9	0.7

Table Q .

 Q Summary of the soil hydraulic properties estimated from the Beerkan database.

	Statistic	λc		Ks (mm h -1 )
		(mm)	BEST-steady	OPD	WU2	SSBI	A4
		(Eq. 8)	(Eq. 15)	(Eq. t1)	(Eq. t2)	(Eq. t3)	(Eq. t4)
	N	427	427	427	427	427	427
	min	1.5	1.3	1.4	1.4	1.6	1.2
	max	737.7	3550.9	3294.1	3493.8	3758.7	3173.7
	mean	112.6	270.7	258.5	269.3	305.2	244.6
	median	68.0	156.8	153.3	156.4	182.4	142.0
	CV	101.3	142.0	140.1	141.3	137.6	141.3

† Department of Agricultural, Food and Forest Sciences (SAAF = Scienze Agrarie, Alimentari e Forestali).

N = Number of Beerkan infiltration experiments; D (cm) = ring diameter; V (mL) = water volume applied with each pouring.

Funding: This work was supported through the INFILTRON Project (ANR-17-CE04-0010, Package for assessing infiltration & filtration functions of urban soils in stormwater management; https://infiltron.org/) funded by the French National Research Agency (ANR), and the "Programma Operativo Nazionale (PON) Ricerca e Innovazione 2014-2020 (Linea 1 -Mobilità dei ricercatori, AIM1853149, CUP: J54I18000120001) funded by the European Regional Development Fund (ERDF) and the Italian Ministry of Education, University and Research (MIUR).

Author Contributions: S. Di Prima outlined the investigation, developed the theory and analyzed the data. All authors contributed to discussing the results and writing the manuscript.

Conflicts of Interest:

The authors declare no conflict of interest.