Duy-Huy Nguyen
email: duyhuy.nguyen@telecom-sudparis.eu

Hang Nguyen
email: hang.nguyen@telecom-sudparis.eu

Renault Éric
email: eric.renault@telecom-sudparis.eu

Samovar

E-MQS -A new Downlink scheduler for Real-time flows in LTE network

Keywords: Scheduling algorithm, Resource allocation, Realtime services, LTE, QoS, E-model, QoE, MQS

This paper proposes a new scheduling scheme which based on the extended E-model, Channel-and QoS-Aware (known as E-MQS scheduler) for real-time traffics in LTE downlink direction. The real-time services (VoIP, Video, etc.) are very sensitive to network impairments such as delay, packet loss, jitter, etc. The proposed scheduling scheme is based on the extension of the E-model and the consideration of Maximum Queue Size (MQS) as a factor for the metric. Since this scheduling scheme considers Mean Opinion Score (MOS) values, thus, it gets higher user perception. The simulation results show that the proposed scheme has the performance which not only satisfies QoS requirements of real-time services but also outperforms the Frame Level Scheduler (FLS), Modified Largest Weighted Delay First (M-LWDF) and Exponential/Proportional Fair (EXP/PF) schedulers in terms of delay, cell throughput, Fairness Index (FI) and Spectral Efficiency (SE), especially for Video flow. Our proposed scheduler also significantly improves the Packet Loss Rate (PLR) in comparison with the M-LWDF and EXP/PF schedulers for both VoIP and Video flows. The performance evaluation is compared in terms of Delay, PLR, Throughput, FI and SE for FLS, M-LWDF, EXP/PF schedulers and our proposed one.

I. INTRODUCTION

The Long Term Evolution (LTE) is developed by the Third Generation Partnership Project (3GPP) [1]. It is a mobile network which has high data rate, low delay and fully packetbased. It means to improve the capability of legacy system by increasing data rates and extending superior Quality of Service (QoS) for various multimedia applications. Basic components of LTE network include a powerful eNodeB (eNB) station and several User equipments (UEs) in addition to a gateway [START_REF] Ali | A utility based resource allocation scheme with delay scheduler for lte service-class support[END_REF]. The eNB station combines with core network through several standard complicated protocols. Basic packet scheduling is carried out by the network operator in both UE and eNB station for both uplink as well as downlink. However, according to the 3GPP, there are no firm specifications for scheduling technique in LTE network. One of the most important modules of packet scheduling is Radio Resource Management (RRM) which decides users that would transmit their data on the air interface. The packet scheduling should integrate fairness in terms of throughput as well as the service policies to which users subscribe [START_REF] Alasti | Quality of service in wimax and lte networks [topics in wireless communications[END_REF].

In this paper, we propose a new downlink scheduling scheme for real-time communications in LTE network with the consideration of user perception. There are several papers which mentioned it. Specifically, authors in [START_REF] Alfayly | Qoe-based performance evaluation of scheduling algorithms over lte[END_REF], [START_REF] Alfayly | Qoe-driven lte downlink scheduling for voip application[END_REF] proposed a new Quality of Experience (QoE)-driven LTE Downlink Scheduling for VoIP Application that is based on QoE min (i.e. MOS score at least equals 3.5 for VoIP Application) and they optimize number of users to access a cell. In [START_REF] Nasimi | Qoe-oriented cross-layer downlink scheduling for heterogeneous traffics in lte networks[END_REF], the authors a cross-layer design scheme that jointly optimises three different layers of wireless protocol stack, namely Application, Medium Access Control (MAC) and Physical layers. The purpose of this paper is to maximize network resource utilization and user-perceived quality of service (also called Quality of Experience -QoE). These papers proposed a new LTE Downlink Scheduler but not based on the E-model. The paper which is closest to our paper described in [START_REF] Mushtaq | Qoe-based lte downlink scheduler for voip[END_REF]. In this paper, the authors proposed a new scheduling scheme for VoIP service in LTE networks by using the user satisfaction as a metric for their scheduler. The authors used the E-model to predict user perception via MOS score, and then this factor used in the metric for scheduling decision. However, in this paper, authors did not consider the impact of network jitter for E-model. In addition, in the metric of the scheduler, there is no the presence of the MQS, and the authors evaluated only for VoIP traffic. In our paper, we extend the idea in [START_REF] Mushtaq | Qoe-based lte downlink scheduler for voip[END_REF] by extending the E-model and propose to consider the MQS as an essential and effective factor for the metric. We used the extended E-model to predict the MOS score and use this score as a main factor in the metric. Besides, we see that, the MQS factor has significant effects on the system performance. In the LTE-Sim [START_REF] Piro | Simulating lte cellular systems: an open-source framework[END_REF], this factor is fixed equal to 0. This means the MQS is infinite. So that, in the scheduling process, the MQS is not considered. However, in fact, the MQS should be a finite value because if the MQS value is infinite then the delay will increase and the congestion could be increased. Therefore, the MQS needs to be considered as a essential factor in the metric of the scheduling algorithms.

The rest of this paper is organized as follows: Overview of the system model is described in section II. In section III, we present the proposed scheduling scheme. The simulation results and performance evaluation of the proposed scheduler are analysed in section IV. The conclusion and future work is represented in section V.

II. THE SYSTEM MODEL

A. The E-model E-model is a computational model developed and standardized by ITU-T [START_REF]Itu-t recommendation g.107: The e-model, a computational model for use in transmission planning[END_REF]. It is used to estimate the MOS for narrow band audio quality. The output of the model is R-factor. The values of this R-factor in range of 0-100 where 100 is the best and 0 is the worst quality. And then, it is mapped to the corresponding MOS value. The standard R-factor in the E-model is defined as follows:

R = R 0 -I s -I d -I ef + A (1)
In which: R 0 : The basic signal-to-noise ratio which consists of noise sources such as circuit and room noise. In this model, its value is set to 94.2. I s : The simultaneous impairment factor, it is the sum of all impairments which may occur more or less simultaneously with the voice transmission. In this model, the default value is set to 0. I d : The delay impairment factor, representing all impairments due to delay of voice signals. I ef : The equipment impairment factor, capturing the effect of signal distortion due to low bit rates of the codec and packet losses of random distribution. A: The advantage factor, capturing the fact that some users can accept a reduction of quality due to the mobility of cellular networks. In this model, this factor is set to 0.

In above factors, I d and I ef are affected by end-to-end delay and packet loss, respectively, while R 0 and I s do not depend on network performance. The R-factor is then translated into the MOS as follows [START_REF]Itu-t recommendation g.107: The e-model, a computational model for use in transmission planning[END_REF]:

M OS =          1, if R < 0 1 + 0.035 × R + 7 × 10 -6 × R × (R -60)× (100 -R), if 0 ≤ R ≤ 100 4.5, otherwise (2)
After setting the default values for the E-model, Equation (1) can be rewritten as follows:

R = 94.2 -I d -I ef (3)
We see that, when voice packet transmitted over an IP network, is is affected by many network impairments such as PLR, delay, jitter, etc. In the E-model, there is no presence of network jitter. In order to improve use satisfaction, we propose to add the I j factor to the E-model. Hence, the E-model can be described as the following formula:

R = 94.2 -I d -I ef -I j (4)
Equation ( 4) shows that the R-factor depends on end-to-end delay (I d ), total loss probability (I ef ), and network jitter (I j ). Hence, in order to compute the R-factor, we must to count these factors. The I d is a factor which is affected by end-toend delay and is calculated as follows [START_REF] Olariu | Voip quality monitoring in lte femtocells[END_REF]:

I d = 0.024 × d + 0.11 × (d -177.3) × H(d -177.3) (5)
In which: H(x) is the Heavyside function:

H(x) = 0, if x < 0 1, otherwise (6)
In equation ( 5), d represents the total end-to-end delay (or mouth-to-ear delay) of speech packet. It can be calculated via some functions in LTE-Sim tool. The I ef is determined according to packet loss. In order to compute this factor, we use the equation in [START_REF] Mushtaq | Qoe-based lte downlink scheduler for voip[END_REF] as follows:

I ef = λ 1 + λ 2 × ln(1 + λ 3 × e l ) (7)
Where: The λ 1 represents the voice quality impairment factor caused by the encoder, λ 2 and λ 3 represent the effect of loss on voice quality for a given codec. Such that, these factors depend on the voice codec used. In this study, we use LTE-Sim [START_REF] Piro | Simulating lte cellular systems: an open-source framework[END_REF] to simulate. This simulation tool supports only G.729 codec, thus, for this codec, the factors above has values as follows: λ 1 = 11, λ 1 = 40, λ 3 = 10. While e l is the total loss probability (consisting of network and buffer layout) which has the value in range of 0..1. This factor is computed directly via some functions in the LTE-Sim tool.

The I j represents the impacts of network jitter to voice quality. It also depends on the voice codec. In this paper, we use the method proposed in [START_REF] Ding | Speech quality prediction in voip using the extended e-model[END_REF] as follows:

I j = C 1 × H 2 + C 2 × H + C 3 + C4 × e -T /K (8)
In which: C 1 , C 2 , C 3 , C 4 are coefficients, K is time instant. These factors depend on the voice codec, for the G.729 codec, these factors have the values as follows: C 1 = -15.5, C 2 = 33.5, C 3 = 4.4, C 4 = 13.6, K = 30. The factor of T is the fixed buffer size of the voice codec. For the G.729 codec, the packet size is 20ms, thus, we select T = 40ms. The H is a factor of Pareto distribution and in range of 0.55 to 0.9. According to [START_REF] Ding | Speech quality prediction in voip using the extended e-model[END_REF], the MOS drops when H increases, thus, in this study, we select H = 0.6 for the simulation.

The final expression of the R-factor when utilizing the G.729 codec is described in Equation [START_REF]Itu-t recommendation g.107: The e-model, a computational model for use in transmission planning[END_REF]. (9) The R-factor is then mapped to the MOS via equations [START_REF] Ali | A utility based resource allocation scheme with delay scheduler for lte service-class support[END_REF] . MOS is one of the important factors for the metric in our scheduler.

B. The correlated Scheduling algorithms

We assume that the metric assigned to the stream i on j-th sub-channel is noted by w i,j . In order obtain the metric, the scheduler usually need to know the average transmission rate ( Ri ) of flow i, and the flow rate available to the UE on the j-th sub-channel. In particular, at each TTI, the estimate Ri is given by [START_REF] Piro | Simulating lte cellular systems: an open-source framework[END_REF]:

Ri (k) = 0.8 × Ri (k -1) + 0.2 × ri (k) (10)
Where: Ri (k -1): The average transmission data rate estimating at the (k-1)-th TTI. ri (k): The rate allocated to i-th flow during the k-th TTI.

In the following subsections, we will describe the metric of several well-known scheduling algorithms which are related to our proposed scheduling scheme including: FLS, M-LWDF, and EXP/PF schedulers. We select these scheduling because they perform well and are suitable for real-time services.

1) The FLS scheduler: FLS is a two-level scheduling algorithm which are called upper level and lower level. These levels are distinct and communicate with each other to dynamically allocate the RBs to the users. At the upper level, a resource allocation method (called FLS) which uses a Discrete-Time (D-T) linear control theory is performed. FLS defines the amount of data that each real-time source should transmit within a single frame to meet its delay constraint. At the lower level, the algorithm uses Proportional Fair (PF) method to allocate RBs to the users at each TTI with considering the bandwidth requirements of FLS to ensure a good level of fairness among multimedia flows. Also at this layer, the scheduler determines number of TTIs/RBs via that each Real-time source will send its packets. In order to calculate the amount of data transmitted, the FLS scheduler uses the following formula:

V i (k) = h i (k) * q i (k) (11)
In which: V i (k) is the amount of the data transmitted by the flow i in LTE frame k, '*' operator is the discrete time convolution, q i (k) is the queue level. It can be said that, V i (k) is computed by filtering the signal q i (k) via a time-invariant linear filter with pulse response h i (k).

2) The M-LWDF scheduler: M-LWDF scheduling algorithm is used to support multiple real-time services in CDMA-HDR systems [START_REF] Ameigeiras | Performance of the mlwdf scheduling algorithm for streaming services in hsdpa[END_REF]. For each real-time flow, by considering the maximum time τ i , the probability is defined as the maximum probability δ i which is the time of the first packet of the queue exceeds the fixed maximum time D HOL,i . In this algorithm, the metrics for real-time and non real-time services are different. In order to offer priority to real-time flows, the metric was given as follows:

w i,j = α i × D HOL,i × ri,j Ri (12)
Where:

• r i,j : The rate assigned to i-th flow during the k-th TTI • Ri : The average transmission data rate estimating • α i : A factor and is given by: α i = -log(δi) τi

3) The EXP/PF scheduler: EXP/PF is a scheduling algorithm which supports multimedia applications in an adaptive modulation and coding and time division multiplexing (AMC/TDM) system [START_REF] Rhee | Scheduling of real/nonreal time services: adaptive exp/pf algorithm[END_REF]. The main purpose of this scheduling algorithm is to enhance the priority for the real-time flows by adding the average fixed maximum time of all active real-time flows. For the real-time services, they receive the increased priorities when their HOL packet delays are approaching the delay deadline. The metric of EXP/PF is calculated as follows:

w i,j = exp( αi×D HOL,i -X a+ √ X ) × ri,j Ri ( 13
)
Where X is given by: X = 1 Nrt ×α i × D HOL,i , with N rt is the number of active real-time flows in the downlink direction.

The remaining parameters are similar to the descriptions above.

III. THE PROPOSED SCHEDULING SCHEME

In our proposed scheduling scheme, we consider the characteristics of the real-time services such as VoIP, Video, etc. These services are sensitive to packet loss and delay, thus, scheduling process should consider various factors. In the related scheduling algorithms above, the authors almost focused on Head of Line packet delay, virtual token length besides other factors such as α i , r i,j and Ri . MOS is a parameter which represents user perception, thus, it should appear in the metric of scheduling algorithms. The higher MOS, the higher user satisfaction. MOS needs be automatically calculated at the receiver and is sent to the eNodeB via feedback technique. For the MQS, according to our knowledge, there are no articles which mention about it. We think that, this factor has strong effects on the system performance. In the LTE-Sim [START_REF] Piro | Simulating lte cellular systems: an open-source framework[END_REF], this factor is fixed equal to 0. This means the MQS is infinite. Hence the MQS is not considered in the scheduling process. However, in fact, the MQS should be finite. If the MQS value is infinite then the delay will increase and the congestion could be thus increased. Therefore, the MQS should be considered as a necessary factor in the metric of the scheduling algorithms.

The main idea of our scheduling algorithm is the consideration of user satisfaction (MOS) and the MQS factor (called also Q i,max ) included into the metric of the scheduling algorithm. This means the higher MOS and the lower (Q i,max -Q i ) values, the higher priority for the UE. The fixed maximum time D HOL,i and the maximum probability δ i are included in the Equation ( 9) to calculate the factors of I d , I ef , respectively. The metric in our scheduling scheme for the real-time services is defined as follows:

w i,j = M OSi×(Qi,max-Qi) τi × ri,j Ri (14)
Where:

• Q i , τ i , r i,j and Ri have the same significances in the previous formulas. • Q i,max : The MQS of the user i. This value can be obtained in bytes via some functions in LTE-Sim [START_REF] Piro | Simulating lte cellular systems: an open-source framework[END_REF]. For the non real-time services, we propose to use the method in the PF scheduler [START_REF] Jalali | Data throughput of cdma-hdr a high efficiency-high data rate personal communication wireless system[END_REF]. The w i,j is a priority matrix for each RB j is assigned to U E i . It is calculated based on the MOS, the remaining queue size (Q i,max -Q i ), the maximum time τ i and the channel condition. MOS is computed at the receiver and is feedbacked to the eNodeB in order to make the scheduling decision of UE. MOS included in the metric will fully exploit the user perception.

IV. SIMULATION ENVIRONMENT AND PERFORMANCE EVALUATION

A. Simulation parameters

The basic parameters used in the simulation are represented in the Table I.

B. Performance evaluation

In order to evaluate the performance of our scheduling scheme, we use the LTE-Sim [START_REF] Piro | Simulating lte cellular systems: an open-source framework[END_REF] to simulate the proposed scheduler with the other schedulers including the FLS, M-LWDF and EXP/PF. The performance evaluation is compared in terms of delay, PLR, cell throughput, FI and SE. The analysis of the simulation results are represented in the following subsection.

1) Delay: End-to-end delay (called one-way delay) is the time required for a packet to be transmitted from source to destination in the network. Figure 1(a) illustrates the delay of VoIP flow. It is clear that, the M-LWDF and EXP/PF scheduler has the same lowest delay. Our proposed scheduler (E-MQS) has the delay which is not significantly higher than the M-LWDF and EXP/PF schedulers and it outperforms the FLS scheduler. For the video flow, all schedulers have the delay which slightly increase when the UEs increase as shown in Figure 1(b). Our proposed scheduler, the M-LWDF and EXP/PF schedulers have nearly the same delay. The FLS has the delay which is higher than the other schedulers. However, all the schedulers have the good end-to-end delay when the number of UE increases equal to 25. It can be included that all the schedulers are also very suitable for video flow.

2) Packet Loss Rate: PLR shows the failure of one or more transmitted packets to reach their destination across a network. For the Video flow, as shown in the Figure 2(b), also as for VoIP flow, the PLR decreases when the number of UE increases. However, for all the schedulers, the PLR under 1% when the number of UE more than 12. For the FLS and E-MQS schedulers, the PLR less than 1% for all number of UE. For the number of UE from 5 to 12, the proposed scheduler has the lowest PLR, when the number of UE greater than 12, the FLS scheduler has the lowest PLR and the E-MQS keeps the second position. In general, the FLS scheduler has the lowest PLR and the M-LWDF has the highest PLR.

3) Cell throughput: As shown in Figure 3(a), for the VoIP flow, the cell throughput of all the schedulers increases when the number of UE increases. The proposed scheduler always has the cell throughput in the top of 2 highest schedulers for all cases of the number of the UE. This means the proposed scheduler is very suitable for VoIP flow. shows the FI of the schedulers. The FI changes when the number of UE increases. It is clear that, the FI of all the schedulers is nearly the same. This is also a very interesting case in the proposed scenario. This not usually happens. This demonstrates that all the schedulers consider the FI is very important in the metric of them.

For the FI, it can be conclusive that the FLS scheduler has the best performance. This is due to at the lower layer of this scheduler uses PF algorithm, thus, it ensures the good grade of fairness among multimedia flows. The FI is high which demonstrates that the cell-edge users have been guaranteed the minimum performance. Hence these users can be served when they move to the edge of the cell.

5) Spectral efficiency: The successful usage of radio resources is a basic purpose of scheduling algorithms. The SE is seen as the performance measurements for the entire cell. As shown on Figure 5, the SE increases when the number of UE increases. In almost case of the number of UE, the proposed and the FLS schedulers have the same highest SE. However, there is no the difference among schedulers. In this paper, we propose a new Channel-, QoS-and QoE-Aware scheduling scheme for the downlink direction in LTE network. The main idea in the proposed scheduler is the consideration of the MOS and the MQS factors into the metric in the proposed scheduler. The metric is based on the user perception (MOS), the remaining queue size, the fixed maximum time, and the channel condition. The simulation results show that the proposed scheduler not only meets QoS requirements for real-time services but also outperforms the FLS, M-LWDF, EXP/PF schedulers in terms of delay, cell throughput, FI and SE, especially for Video flow. For the PLR in the case of VoIP flow, the proposed scheduler has the performance not as well as the FLS scheduler. For the FI in the case of VoIP flow, the proposed scheduler does not have the stable performance it obtains or the highest either the lowest FI values. The proposed scheduler also significantly improves the delay in comparison with the M-LWDF and EXP/PF when the number of UE more than than 12 for Video flow. It can be said that when considering the MOS and the MQS as factors for the metric in the proposed scheduler, the system performance has been improved significantly, specially for Video flow. Through all simulation results, it can be said that, the FLS has the best performance, and the proposed scheduler keeps the second position in scheduling for both VoIP and Video flows. So that, the proposed scheduler is very suitable for real-time services such as VoIP, Video, etc. for the downlink direction in LTE system.

In this study, we have not yet considered the presence of best-effort flow. They might bring new directions to evaluate the proposed scheduler in heterogeneous traffic in LTE network.

R

  = 64.28 -[0.024 × d + 0.11 × (d -177.3)× H(d -177.3)] -40 × ln(1 + 10 × e l ) -13.6 × e -4/3

  Fig. 1. Effects of Delay vs VoIP and Video users

Figure 2 (

 2 a) represents the PLR of VoIP flow. When we set the MQS equal to 10 5 bytes (this factor is set equal to 0 as default in the LTE-Sim), for the VoIP flow, all schedulers have the decreased PLR when the number of UE increases. Normally, the PLR increases when the number of UE increases, thus, this case is quite special and it represents the unstableness of a real system. As shown in the Figure2(a), all the schedulers have the PLR which are under 1% while the FLS has the lowest PLR. The E-MQS scheduler has the second position while the M-LWDF has the highest PLR.

Fig. 2 .

 2 Fig. 2. Effects of PLR vs VoIP and Video users

Fig. 3 .

 3 Fig. 3. Effects of Throughput vs VoIP and Video usersFor the Video flow, in the proposed scenario, all schedulers have the near values of cell throughput for all cases of the number of UE as shown on Figure3(b). This is very interesting case. This concludes that all the schedulers are very conformable to Video flow. Also as VoIP flow, the cell throughput increases when the number of UE increases. This is very good for real-time flows such as VoIP, Video, etc.4) Fairness index: For the VoIP flow as shown on Figure4(a), the FIs of all schedulers are not stable when the number of UE increases. However, all these schedulers have the high FIs. The proposed scheduler has the highest FI when the number of UE equals 5, 15, and 20 and obtains the lowest FIs in the case of the number of UE equals 10, and 20. In general, the FLS has the best FI for VoIP flow.

Fig. 4 .

 4 Fig. 4. Effects of Fairness Index vs VoIP and Video users

Fig. 5 .

 5 Fig. 5. Spectral efficiency vs number of user V. CONCLUSION


Leroy Boris
email: leroy.boris@gmail.com

Gallon Régis

Feunteun Eric

Robuchon Marine

Ysnel Frédéric

Boris Leroy

Cross-taxon congruence in the rarity of subtidal rocky marine assemblages: no taxonomic shortcut for conservation monitoring

Keywords: Occurrence-based rarity, subtidal conservation, rare species, Brittany

published or not. The documents may come

Introduction

To mitigate the ongoing loss of marine biodiversity, the European Union has successively adopted several directives aiming at protecting and conserving marine habitats, ecosystems and biological diversity: the OSPAR convention (1998), the habitat directive (HD, 92/43/EEC, 1992), the Natura 2000 network, and more recently the Marine Strategy Framework Directive (MSFD; 2008/56/EC) (2010). Effectively preserving marine biodiversity requires the evaluation and monitoring of the diversity of different marine taxonomic groups and habitats. With respect to benthic subtidal habitats, much of the literature has focused on subtidal soft bottom habitats for which a plethora of tools is available [START_REF] Van Rein | A review of sublittoral monitoring methods in temperate waters: a focus on scale[END_REF]. On the other hand, rocky subtidal habitats remain poorly evaluated because of the difficulties and costs of sampling these diversified habitats. This lack of consideration impedes their inclusion in the implementation of European directives, in spite of their high taxonomic diversity and functional role. Consequently, researchers need to both optimise sampling procedures [START_REF] Gallon | Optimization of an "in situ" subtidal rocky-shore sampling strategy for monitoring purposes[END_REF] and identify potential surrogates that could be used to describe biodiversity based on a reduced set of data. In this study, we seek to identify biodiversity surrogates as one facet of the importance of subtidal rocky habitats for conservation, hereby not focusing on other equally important facets such as ecosystem functioning.

The concept of biodiversity surrogates has been extensively explored in the conservation literature and needs to be refined to be relevant to rocky subtidal habitats. Biodiversity surrogates fall in two main categories [START_REF] Grantham | Effectiveness of biodiversity surrogates for conservation planning: different measures of effectiveness generate a kaleidoscope of variation[END_REF]: environmental surrogates, which use a combination of physical and biological data to estimate or predict biodiversity, and taxonomic surrogates, predominantly based on the use of one or several taxonomic groups to estimate biodiversity. We focus here on the second category, i.e. the identification of surrogate phyla by cross-taxon congruence, since it has been advocated as appropriate to predict patterns of marine biodiversity for conservation purposes when data are scarce [START_REF] Mellin | Effectiveness of biological surrogates for predicting patterns of marine biodiversity: a global meta-analysis[END_REF]. The identification of surrogate taxa in marine habitats has mostly been based on species richness [START_REF] Mellin | Effectiveness of biological surrogates for predicting patterns of marine biodiversity: a global meta-analysis[END_REF] and, to a lesser extent, on multivariate patterns of assemblage structure (e.g., [START_REF] Hirst | Surrogate measures for assessing cryptic faunal biodiversity on macroalgaldominated subtidal reefs[END_REF][START_REF] Smale | Monitoring marine macroalgae: The influence of spatial scale on the usefulness of biodiversity surrogates[END_REF][START_REF] Sutcliffe | Biological surrogacy in tropical seabed assemblages fails[END_REF]). However, it has been suggested that species richness is not appropriate because of high spatial and taxonomic variability [START_REF] Su | Beyond Species Richness: Community Similarity as a Measure of Cross-Taxon Congruence for Coarse-Filter Conservation[END_REF]. More importantly, species richness is not a good enough indicator for biodiversity conservation as it does not take the identity of species into account or their varying degrees of vulnerability to extinction (e.g., [START_REF] Orme | Global hotspots of species richness are not congruent with endemism or threat[END_REF]. Alternatively, surrogacy across taxa can be assessed by studying congruency of rarity across assemblages of species, an appropriate approach for taxa with limited data availability [START_REF] Leroy | Integrating multiple scales in rarity assessments of invertebrate taxa[END_REF][START_REF] Leroy | Improving occurrence-based rarity metrics in conservation studies by including multiple rarity cut-off points[END_REF]. The choice of rarity as a criterion is based on the greater extinction risk of rare species relative to ecologically similar common species [START_REF] Flather | Species rarity: definition, causes and classification[END_REF][START_REF] Gaston | Rarity[END_REF][START_REF] Roberts | Extinction risk in the sea[END_REF], because they are appropriate indicators for other species of conservation concern [START_REF] Larsen | Improving the performance of indicator groups for the identification of important areas for species conservation[END_REF][START_REF] Lawler | Rare Species and the Use of Indicator Groups for Conservation Planning[END_REF], and because they have been shown to sometime support unique ecosystem functions unsupported by other species [START_REF] Mouillot | Rare species support vulnerable functions in high-diversity ecosystems[END_REF].

It has been established repeatedly that a large proportion of intertidal and subtidal marine benthic species exhibit narrow geographic ranges (see for instance [START_REF] Sanderson | Rarity of marine benthic species in Great Britain: development and application of assessment criteria[END_REF][START_REF] Chapman | Are there adequate data to assess how well theories of rarity apply to marine invertebrates?[END_REF]. Therefore the frequency distributions of benthic marine species range size are typically strongly right-skewed with a large number of low occurrence species at local or regional scales [START_REF] Ellingsen | Rare species, habitat diversity and functional redundancy in marine benthos[END_REF] similar to terrestrial taxa [START_REF] Gaston | Rarity[END_REF][START_REF] Gaston | Pattern and Process in Macroecology[END_REF]. Hence, testing the ability of the rarity of surrogate taxa to capture the rarity of other taxa may help to substantially reduce the cost of sampling protocols to evaluate the rarity of subtidal assemblages.

The methods used to measure rarity in species assemblages must be chosen carefully, because several methods have been proven to provide inappropriate results under particular conditions [START_REF] Leroy | Improving occurrence-based rarity metrics in conservation studies by including multiple rarity cut-off points[END_REF]. Hence, [START_REF] Leroy | Improving occurrence-based rarity metrics in conservation studies by including multiple rarity cut-off points[END_REF] proposed a new, flexible method to assess the rarity of species assemblages (the Index of Relative Rarity) with respect to the considered phyla. The rationale of this method is that rarity should be defined according to the taxon considered; therefore an inflexible method may bias the analysis toward a particular phylum. Indeed, the threshold of geographic range size below which species are considered rare is generally defined specifically for each particular phylum because of the large differences in range size among phyla [START_REF] Grenyer | Global distribution and conservation of rare and threatened vertebrates[END_REF]. This rationale is especially important for subtidal rocky habitats because the diversity of their phyla may result in very different rarity patterns among phyla.

We aimed to investigate patterns of rarity and cross-taxon congruency among sessile and lowmobility animal phyla of subtidal rocky areas, based on the Index of Relative Rarity developed by [START_REF] Leroy | Improving occurrence-based rarity metrics in conservation studies by including multiple rarity cut-off points[END_REF]. The opportunity for such an approach was offered by the compilation of a database from 137 inventories of assemblages of subtidal rocky habitats around the Brittany (western France) coast. These inventories were sampled between 1993 and 1998 by biologist scuba divers of the "Association pour la Découverte du Monde Marin" [START_REF] Girard-Descatoire | Inventaire de la faune et de la flore sur les fonds rocheux de la Baie de Lannion[END_REF][START_REF] Chapman | Are there adequate data to assess how well theories of rarity apply to marine invertebrates?[END_REF][START_REF] Girard-Descatoire | Inventaire de la faune et de la flore sur les fonds rocheux de l'archipel des Sept-Iles[END_REF][START_REF] Girard-Descatoire | Inventaire de la faune et de la flore sur les fonds rocheux de Saint-Malo et de Dinard. Rapport Diren[END_REF], 1996a, 1996b, 1996c[START_REF] Girard-Descatoire | Inventaire de la faune et de la flore sur les fonds rocheux autour de l'île d'Ouessant[END_REF][START_REF] Girard-Descatoire | Inventaire des fonds rocheux de grand intérêt biologique en milieu marin dans la zone subtidale (Côte de Granite Rose)[END_REF]L'Hardy-Halos et al., 2001;L'Hardy-Halos and Castric-Fey, 2000a, 2000b, 2001) in addition to samples by our diving team from 2005 to 2009. Because this database was initially not designed for such analyses, it may contain biases that could have a negative impact on the outcomes of our study [START_REF] Pearman | Identifying potential indicators of conservation value using natural heritage occurrence data[END_REF]. Hence, we applied a completeness metric to assess sampling quality across sites and phyla [START_REF] Soberón | Assessing completeness of biodiversity databases at different spatial scales[END_REF], and improved our database by applying appropriate corrections.

The objectives of this study are to [START_REF] Van Rein | A review of sublittoral monitoring methods in temperate waters: a focus on scale[END_REF] investigate the distribution of rarity among the sampled phyla, in order to calculate appropriate rarity metrics for each phylum;

(2) search for potential surrogate phyla with a cross-taxon congruence approach on the rarity of assemblages of species;

(3) propose a multi-phyla indicator and discuss its potential use to evaluate the importance of monitored subtidal rocky habitats for conservation.

Material and Methods

Database compilation

We compiled a database on the biodiversity of benthic assemblages of species of subtidal rocky habitats around Brittany based on 122 inventories sampled between 1993 and 1998 by biologist scuba divers [START_REF] Girard-Descatoire | Inventaire de la faune et de la flore sur les fonds rocheux de la Baie de Lannion[END_REF][START_REF] Chapman | Are there adequate data to assess how well theories of rarity apply to marine invertebrates?[END_REF][START_REF] Girard-Descatoire | Inventaire de la faune et de la flore sur les fonds rocheux de l'archipel des Sept-Iles[END_REF][START_REF] Girard-Descatoire | Inventaire de la faune et de la flore sur les fonds rocheux de Saint-Malo et de Dinard. Rapport Diren[END_REF], 1996a, 1996b, 1996c[START_REF] Girard-Descatoire | Inventaire de la faune et de la flore sur les fonds rocheux autour de l'île d'Ouessant[END_REF][START_REF] Girard-Descatoire | Inventaire des fonds rocheux de grand intérêt biologique en milieu marin dans la zone subtidale (Côte de Granite Rose)[END_REF]L'Hardy-Halos et al., 2001;L'Hardy-Halos and Castric-Fey, 2000a, 2000b, 2001), and 15 inventories sampled by our diving team from 2005 to 2009.

The inventories were grouped into 13 major sites around the Brittany coast, with 4 to 20 inventories per site. The 13 sites span the range of conditions around Brittany: estuaries and inland seas (Morbihan gulf, Etel ria, Brest bay, Rance estuary), open bays (Lannion bay, St Malo bay, Morlaix bay, Iroise sea), coastal sites (Granite rose coast, Crozon peninsula, Cape Sizun) and islands (Sept-Îles archipelago, Ushant island). These inventories of species presence-absence focused on benthic rocky communities between the infralittoral fringe and the nearest circalittoral, and were limited to depths of 30 m.

In this substantial inventory, eight animal taxa (conveniently called "phyla" in this paper) were kept for our analyses: Porifera, Cnidaria, Annelida, Mollusca, Arthropoda, Bryozoa, Echinodermata and Chordata, for a total of 548 species. This choice focused on only sessile and low-mobility species to minimise bias in scuba-diving surveys throughout the 13 investigated sites. The systematic nomenclature of the database was checked following the World Register of Marine Species (WoRMS Editorial Board, 2016) to avoid both orthographic mistakes and synonymies [START_REF] Costello | European register of marine species. A check-list of the marine species in Europe and a bibliography of guides to their identification[END_REF].

Database robustness and occurrence estimation

We first analysed the robustness of the database to assess and mitigate potential biases in sampling efforts, which could lead to an uneven representation of species diversity and occurrence among sites. Our analysis was divided into a three step process based on a completeness index [START_REF] Soberón | Assessing completeness of biodiversity databases at different spatial scales[END_REF]. The completeness index assesses the completeness of a set of samples by dividing the observed richness by the total estimated species richness on the basis of a richness estimator [START_REF] Soberón | Assessing completeness of biodiversity databases at different spatial scales[END_REF]. We calculated three richness estimators (Chao2, ICE and Jack1) [START_REF] Hortal | Evaluating the performance of species richness estimators: Sensitivity to sample grain size[END_REF][START_REF] Soberón | Assessing completeness of biodiversity databases at different spatial scales[END_REF] which yielded similar results; we included average values across the three estimators in the main text, and values for all three indicators in Appendix A.

Firstly, we analysed the completeness of the whole database. Secondly, all phyla pooled together, we analysed the completeness of each site in order to identify and remove under-sampled sites with a completeness threshold (see below). Thirdly, we analysed the completeness of each phylum in each site. Then we calculated the average completeness of each phylum across all sites in order to remove phyla that were on average insufficiently sampled across all sites. We applied a conservative rule to remove phyla and sites based on a threshold of completeness that we defined at 75%, i.e. when the observed richness did not reach 75% of estimated richness for a particular site or phylum and for at least two richness estimators, this site or phylum was removed.

The number of sampled stations differed between sites, from 4 to 20 sampled stations (Table 1). This difference introduced a bias if the occurrence was calculated from the number of sampled stations: a species occurring in all stations of a particular site would receive a higher occurrence if the site had 20 sampled stations than if the site had 4 sampled stations. Hence, this would result in an artificial "commonness" for better-sampled sites, versus an artificial rarity for lesssampled sites. Therefore, to remove this bias, we estimated the frequency of occurrence of each species in each site with the following formula:

j ij ij N q Q 
where: Qij is the frequency of occurrence of species i in site j; Nj is the number of stations sampled in the site j, and qij is the number of stations of site j in which species i has been found. The total occurrence (Qi) of each species i is then calculated as the sum of its occurrence in all sites:

  ij i Q Q .

Species rarity and calculation of rarity weights

We first analysed the frequency distribution of species occurrences to demonstrate differences in the distribution of rarity among phyla. Then, we calculated rarity weights for each species on the basis of a method that could be adjusted according to a user-chosen rarity cut-off point [START_REF] Leroy | Integrating multiple scales in rarity assessments of invertebrate taxa[END_REF][START_REF] Leroy | Improving occurrence-based rarity metrics in conservation studies by including multiple rarity cut-off points[END_REF]. With this method, rare species receive rarity weights that increase exponentially when their occurrence falls below a rarity cut-off point. Thus, weights of rare species (with occurrence lower than the cut-off) are amplified, whereas weights of common species (with occurrence higher than the rarity cut-off) tend to zero. At the cut-off point, species weight is always equal to 5% of the weight of the rarest species [START_REF] Leroy | Integrating multiple scales in rarity assessments of invertebrate taxa[END_REF]. The rarity cut-off point should be defined specifically for each phylum rather than choosing a general cutoff for all phyla, especially when phyla have different life histories [START_REF] Flather | Species rarity: definition, causes and classification[END_REF].

The rarity weight wi of each species i was calculated with the following formula [START_REF] Leroy | Integrating multiple scales in rarity assessments of invertebrate taxa[END_REF]:

                       2 min max min 05 . 1 97 . 0 exp Q Q r Q Q w i i (1)
Where all parameters were defined for the considered phylum: Qi is the occurrence of species i; Qmin and Qmax, minimum and maximum occurrences, respectively, of the considered phylum; and r, chosen rarity cut-off point (as a percentage of maximum occurrence) for the considered phylum.

To select a rarity cut-off point for each phylum, we followed Gaston's recommendations (1994): the rarity cut-off was the first quartile of the frequency distribution of species occurrences (i.e. rare species are the 25% species with the lowest occurrence). We analysed the effect of this cutoff choice on our analyses with a sensitivity analysis (Appendix B).

Index of Relative Rarity

The IRR of an assemblage of species is calculated as the average weight of rarity of all the species of the assemblage. The IRR is subsequently normalised between 0 and 1:

min max min 1 w w w S w I S i i RR      (2)
Where: wi is the weight of the i th species of the assemblage; S, species richness; and wmin and wmax, minimum and maximum possible weights, respectively. The IRR ranges from 0 (all species of the assemblage have the minimum weight, i.e. ubiquitous species) to 1 (all species of the assemblage have the maximum weight, i.e. very rare species).

Analyses

For each of the 133 inventories of our database, we calculated phylum-specific IRR for each of the 8 phyla. To search for potential surrogate phyla, we analysed the congruence between each pair of phyla with the phylum-specific IRR. Given the nested nature of the sampling design, we could not directly analyse the correlation among the 133 inventories. We therefore analysed the correlations at two levels: between sites and within sites. For the between sites level, we tested for cross-phyla correlations using mean IRR values for each site. For the within sites level, we tested for cross-phyla correlations within each site and reported the average within-site correlation, as well as the percentage of sites for which a significant correlation was detected. Our expectations were that a good surrogate phylum should be correlated to other phyla within each site and between sites as well. The congruence was measured with the rank-based correlation coefficient of Spearman. We finally constructed and analysed a multi-phyla index of rarity with phylum-specific rarity cut-offs (see Appendix C).

All data analyses were performed using R (R Core Team, 2016); completeness analyses were done with the "fossil" package [START_REF] Vavrek | fossil: palaeoecological and palaeogeographical analysis tools[END_REF], rarity analyses with the "Rarity" package [START_REF] Leroy | Rarity: Calculation of Rarity Indices for Species and Assemblages of Species[END_REF], correlation analyses with the "psych" package [START_REF] Revelle | psych: Procedures for Personality and Psychological Research[END_REF] and graphics with "ggplot2" [START_REF] Wickham | ggplot2: Elegant Graphics for Data Analysis[END_REF], all available on the Comprehensive R Archive Network.

Results

Database analysis

The completeness of the whole database was very high with 93% of the estimated total species richness sampled (Table 1). On average, 265±54 (mean±sd) species were sampled per site, with an average estimated species richness of 314±56. The completeness of almost all sites was relatively high with an average completeness value of 0.84±0.08). The only under-sampled site was the Rance estuary, with an estimated completeness of 0.60 and 4 sampled sites. Hence, the Rance estuary was removed for the calculation of species occurrence and cross-taxon congruence analyses. Interestingly, the heterogeneity in sampling intensities across sites did not lead to significant differences in completeness indices (observed richness vs. sampling effort: Pearson's r = -0.09, df = 10, p = 0.78; completeness vs. sampling effort: Pearson's r = 0.53, df = 10, p = 0.07). The average completeness per site of phyla ranged from 0.81 to 0.91, indicating that all these phyla were sampled with similar intensities for the 13 sites (Table 2). Out of the remaining phyla, the average number of species per inventory and per phylum ranged from 4.7±3.3 for Annelida to 21.1±5.5 for Mollusca. The average number of species per site and per phylum ranged from 11.4±6.3 for Annelida to 46.3±21.6 for Cnidaria. The 8 phyla were found in all of the sites with similar richness proportions among sites (See Appendix D), with three major phyla being Bryozoa, Cnidaria and Porifera.

Rarity distribution among phyla

Corrected occurrences (see methods) of species ranged from 0.05 to 12, and their distribution varied between phyla (Figure 1A). As expected, for the majority of phyla, the median was low, indicating that most species had low occurrences. Four phyla had ubiquitous species (i.e., occurrences equal to the maximum: 12): Arthropoda, Bryozoa, Cnidaria, Mollusca (species list available in Appendix E). The rarity cut-off point, defined as the first (lowest) quartile of the frequency distribution of occurrences, varied substantially between phyla, from an occurrence of 0.24 (Mollusca) to an occurrence of 0.80 (Bryozoa) (Figure 1A). Hence, rare species of Mollusca had occurrences ranging from 0.05 to 0.24, while rare species of Bryozoa had occurrences ranging from 0.05 to 0.80. On the other hand, the global rarity cut-off for all phyla implied a definition of rare species for occurrences ranging from 0.05 to 0.44.

Species rarity weights were calculated for each phylum with respect to phylum-specific rarity cut-off points. Rare species of each phylum received weights ranging from 0.05 (weight at the cut-off) to 1 (weight of the rarest species) (Figure 1B; full list of species with rarity scores in Appendix E). Although the range of weights was identical among phyla, several phyla had a greater variety of weights because of a higher dispersion of occurrence values among rare species (e.g., Bryozoa, Chordata and Cnidaria).

Rarity scores of species assemblages

Indices of relative rarity of inventories ranged from 0 to a maximum of 0.25 for an inventory of Annelida at the Ushant Island (Figure 2). A substantial number of inventories had indices of 0 for all phyla, but not necessarily for the same sites. Interestingly, all sites had at least one inventory with a high value for a phylum. Distributions of rarity indices across sites were clearly different between phyla. Some phyla had high rarity values concentrated in only a few sites (Annelida, Arthropoda, Echinodermata, Mollusca); other phyla had high rarity values in many sites (Bryozoa, Chordata, Cnidaria, Porifera).

Cross-taxon congruence of rarity

Between sites, the correlation among phyla of rarity indices was not significant for all comparisons (Figure 3). Within-site correlations were generally low, with the highest values obtained between Arthropoda and Porifera for which 30% of the station-level correlations were significant. This absence of congruence among phyla was reflected in the graphical comparison of inventory indices (Figure 3). For any given phylum, inventories with no rare species could correspond to the best as well as the worst ranked inventories of most other phyla. The diagonal shows the different phyla and the number of inventories available for each.

Multi-phyla indices of rarity

Given that no phylum could be identified as a clear indicator of another phylum rarity, we included all of them in the multi-phyla index of rarity. The multi-phyla IRR had values ranging from 0 to 0.041 (Figure 4). The most important result was that high values of IRR were found in all sites (Figure 4 and Figure 5). In addition, the top 25% assemblages were distributed in almost all of the sites, spread throughout Brittany, with the notable exception of the Morbihan Gulf. Nevertheless, higher values were generally found in north (Sept-Îles archipelago, Lannion bay, Granite rose coast) and western Brittany (Ushant Island, Brest bay). The rarest assemblages were found in both Ushant and Sept Îles islands. Besides, all sites had assemblages with no rare species (IRR equal to 0).

Discussion

This study provides a method for a quantitative assessment of the importance for conservation of subtidal rocky habitats, based on the rarity of sessile and low-mobility species. Although the reference database has limitations and needs to be completed with more observations, with this method we evaluated for the first time, the rarity patterns of invertebrate phyla of the subtidal rocky benthos around Brittany. Our study highlighted three main findings. Firstly, rarity appeared to be unevenly distributed among the studied phyla, which emphasised the need to fit rarity metrics according to each phylum. Secondly, the distribution of rarity in subtidal assemblages of species was not clearly congruent among phyla, which prevented any attempt to identify indicator phyla. Thirdly, according to multi-phyla indices of rarity, assemblages of high rarity were not concentrated within a few sites: they were distributed across all of our studied sites, and were spread throughout our study region. These three findings have important implications both for future research investigating the rarity of different invertebrate phyla, as well as for the conservation and monitoring of subtidal rocky assemblages.

A flexible approach to weigh species rarity

In this study, we investigated the rarity of species and assemblages of species of eight different phyla. We chose to analyse rarity on the basis of the most widespread definition of rarity in conservation literature, i.e. rare species are the 25% species with the lowest occurrence [START_REF] Flather | Species rarity: definition, causes and classification[END_REF][START_REF] Gaston | Rarity[END_REF]. Such a cut-off criterion was chosen because it ensures that all phyla have equal contribution to the ranking of assemblages (Appendix C), thus preserving phylogenetic diversity in the targeted rare species. This cut-off criterion was proven to be appropriate for other marine invertebrate taxa, particularly because it standardises the definition and enables comparisons across assemblages and phyla [START_REF] Benkendorff | Multiple measures are necessary to assess rarity in macro-molluscs: a case study from southeastern Australia[END_REF]. With this criterion, we showed that cut-off values differed between phyla: several phyla had very low cut-offs (e.g., Mollusca), whereas other phyla had relatively high cutoffs (e.g., Bryozoa, Chordata, Cnidaria and Echinodermata). Because of the differences in rarity cut-offs between phyla, a flexible weight function was required to fit the weight assignation with respect to the rarity of each phylum. The function we applied here explicitly integrated the rarity cut-off point, thereby ensuring that assemblages of different phyla with different rarity cut-offs could be compared [START_REF] Leroy | Integrating multiple scales in rarity assessments of invertebrate taxa[END_REF][START_REF] Leroy | Improving occurrence-based rarity metrics in conservation studies by including multiple rarity cut-off points[END_REF]. While we advocate here the use of a phylumspecific definition of rarity, one can argue that the conservation of the rarest species, regardless of their phylogenetic origin, is primordial. In such a case, the use of a single criterion for all phyla may be preferred, which can be done with the application of a single rarity cut-off for all phyla.

Discrepancy in cross-taxon congruence prevents the identification of surrogate phyla

We did not observe any striking congruence in the rarity of assemblages among the studied phyla of the subtidal rocky areas of Brittany. In other words, the spatial distribution of rarity in Brittany differed between the 8 phyla of marine subtidal organisms we studied. The major consequence of this finding is that no single phylum or group of phyla can be chosen as a surrogate to monitor subtidal rocky areas. Furthermore, given the diversity of rarity patterns among phyla, the choice of any restricted set of phyla to monitor subtidal rocky areas for conservation would neglect the rare species of omitted phyla. Consequently, our approach did not permit us to reveal any easy taxonomic shortcut to monitor subtidal rocky areas for conservation around Brittany, which is worrying given the costs of sampling multiple phyla in such biota.

Interestingly, the identification of surrogate phyla by cross-taxon congruence (not limited to rarity) in marine species assemblages yielded contrasting results in the literature. On the one hand, the same incongruence was reported for marine biota similar to the ones we studied [START_REF] Hirst | Surrogate measures for assessing cryptic faunal biodiversity on macroalgaldominated subtidal reefs[END_REF], as well as related marine biota such as intertidal zones in the United Kingdom [START_REF] Reddin | Between-taxon matching of common and rare species richness patterns[END_REF] or structurally less similar biota such as coral reefs [START_REF] Jimenez | Multi-taxa coral reef community structure in relation to habitats in the Baa Atoll Man and Biosphere UNESCO Reserve (Maldives), and implications for its conservation[END_REF] and tropical seabeds [START_REF] Sutcliffe | Biological surrogacy in tropical seabed assemblages fails[END_REF]. On the other hand, studies on a similar set of phyla found that molluscs constituted an appropriate surrogate phylum for rapid assessments of biodiversity [START_REF] Smith | Rapid assessment of invertebrate biodiversity on rocky shores: where there's a whelk there's a way[END_REF]; and neither annelids, arthropods nor molluscs were appropriate surrogate phyla for estuarine conservation [START_REF] Shokri | Annelids, arthropods or molluscs are suitable as surrogate taxa for selecting conservation reserves in estuaries[END_REF]. Thus, the ability to identify surrogate phyla and the identity of these surrogate phyla appear region-dependant. In accordance with [START_REF] Sutcliffe | Biological surrogacy in tropical seabed assemblages fails[END_REF], we therefore recommend that the effectiveness of surrogate phyla requires testing when defining surrogates in a new region. However, even if we recognise that studying rare species requires time-consuming sampling protocols and is thus problematic because of limited funds available for monitoring, we disagree with [START_REF] Sutcliffe | Biological surrogacy in tropical seabed assemblages fails[END_REF] regarding their proposal to exclude rare species for monitoring biodiversity, because rarity is one of the only measures directly characterising species extinction risk for lesserknown taxa [START_REF] Flather | Species rarity: definition, causes and classification[END_REF][START_REF] Leroy | Using biodiversity databases for the conservation of invertebrate taxa: Assemblage-level indices of rarity and species distribution models[END_REF].

Spatial distribution of rarity

On the basis of the multi-phyla index of relative rarity, we showed that assemblages of high rarity were distributed across all the sites around Brittany. Indeed, each site exhibited at least one assemblage of high rarity. Nevertheless, assemblages with rare species occurred more frequently in northern Brittany (Granite rose coast, Sept-Îles archipelago, Lannion bay) and also western Brittany (Brest bay, Ushant island). These sites concentrating potentially vulnerable species can therefore be targeted for further investigation for their conservation. Interestingly, this spatial pattern is similar to the spatial pattern of genetic diversity in two brown macroalgae: [START_REF] Robuchon | Contrasting genetic diversity patterns in two sister kelp species co-distributed along the coast of Brittany, France[END_REF] have shown, using a similar spatial distribution of samples, that the highest genetic diversity was found in both northern and western Brittany, compared to northeastern and southern Brittany. This peculiar pattern may be explained by the seemingly more stable temperatures in north-western Brittany than elsewhere [START_REF] Gallon | Twenty years of observed and predicted changes in subtidal red seaweed assemblages along a biogeographical transition zone: inferring potential causes from environmental data[END_REF]. Northwestern Brittany is also characterised by a widespread and continuous rocky bottom, separated from the more fragmented rocky bottom of the other regions of Brittany by large sandy beaches [START_REF] Cabioch | Contribution à la connaissance des peuplements benthiques de la Manche occidentale[END_REF][START_REF] Méléder | Predictive modelling of seabed habitats: Case study of subtidal kelp forests on the coast of Brittany, France[END_REF][START_REF] Raffin | Bases biologiques et écologiques de la conservation du milieu marin en mer d'Iroise[END_REF][START_REF] Retière | Contribution à la connaissance des peuplements benthiques du golfe normanno-breton[END_REF]. Indeed, habitat fragmentation increases extinction probability [START_REF] Roberts | Extinction risk in the sea[END_REF]. Stable environmental conditions and widespread habitat may therefore offer suitable conditions for species with narrow ecological requirements, an intrinsic cause of species rarity [START_REF] Flather | Species rarity: definition, causes and classification[END_REF], and this may explain the observation of higher concentrations of rarity in north-western Brittany. In addition, this area is the location of the transition zone between two biogeographical provinces (the Lusitanian province in the south and the Boreal province in the north; [START_REF] Cox | Biogeography: an ecological approach[END_REF]Moore 2000, Spalding et al. 2007). Consequently, north-western Brittany is an area where the range limits of both Lusitanian and Boreal species overlap. This overlap may explain the concentration of rarity in this transition zone, since species which are rare in Brittany may be at the edge of their geographical ranges (e.g., [START_REF] Leroy | Integrating multiple scales in rarity assessments of invertebrate taxa[END_REF]. Therefore, north-western Brittany is likely to concentrate rare species, not only because this area exhibits environmental characteristics favourable for specialist species, but also because it constitutes a transition zone between two biogeographical regions, thus increasing the likelihood of having species at the edge of their range.

However, the more general pattern showing that rare species occur everywhere around Brittany is more difficult to explain. Among the few studies investigating the causes of spatial patterns of rarity in the marine benthos, [START_REF] Ellingsen | Rare species, habitat diversity and functional redundancy in marine benthos[END_REF] showed that on soft substrates, habitat characteristics might play an important role. Indeed, they highlighted the fact that the number of rare species was strongly correlated to the number of habitats and environmental variability, suggesting that the number of rare species increased with both within-and between-site heterogeneity, and that these relationships may arise from habitat-specific species with restricted ranges. The coastline of Brittany exhibits a huge mosaic of benthic habitats [START_REF] Bajjouk | Typologie d'habitats marins benthiques: Analyse de l'existant et propositions pour la cartographie. Habitats côtiers de la région Bretagne-Note de synthèse n° 2, Habitats du circalittoral[END_REF][START_REF] Guillaumont | Typologie d'habitats marins benthiques : analyse de l'existant et propositions pour la cartographie (Habitats côtiers de la région Bretagne)-Note de synthèse[END_REF] and displays an important medium-scale (> 100km) environmental heterogeneity [START_REF] Ayata | How does the connectivity between populations mediate range limits of marine invertebrates? A case study of larval dispersal between the Bay of Biscay and the English Channel (North-East Atlantic)[END_REF][START_REF] Gallon | Twenty years of observed and predicted changes in subtidal red seaweed assemblages along a biogeographical transition zone: inferring potential causes from environmental data[END_REF].

Another interesting hypothesis suggests that local water flow direction and velocity drives diversity by mediating the delivery (larval recruitment) of rare species [START_REF] Palardy | Water flow drives biodiversity by mediating rarity in marine benthic communities[END_REF]. The high complexity of flows around Brittany, generating multiple distinct hydrographic areas [START_REF] Ménesguen | Observation and modelling of natural retention structures in the English Channel[END_REF][START_REF] Pingree | Celtic Sea and Armorican current structure and the vertical distributions of temperature and chlorophyll[END_REF][START_REF] Salomon | Courants résiduels de marée dans la Manche[END_REF], may also explain the observed spatial distribution of rarity. To summarise, our study area is characterised by a complex diversity of habitats and environmental conditions, which might explain why rare assemblages seem to occur everywhere in Brittany. However, this hypothesis requires further testing by correlating rarity patterns with environmental patterns, and also the investigation of a higher number of locations following a nested sampling design.

Caveats and future directions

Estimating rarity is highly dependent on uneven sampling intensities among sites, and crosstaxon analyses are highly dependent on uneven sampling between phyla. As a consequence, we decided to make several corrections to improve the quality of the database before conducting cross-taxon analyses of rarity. For example, we defined a quality threshold for our sampling.

While the choice of the threshold value was arbitrary, it has the advantage of ensuring that sampling intensity was similar among phyla and among sites, by the removal of phyla and sites with the poorest sampling intensities. Another correction we applied was the use of a corrected index of occurrence to account for sampling discrepancies among sites. This index has the advantage of addressing the issue of occurrence overestimation in sites with many stations vs. sites with few stations. However, because of its mathematical formulation, rare species from sites with few stations will have higher occurrences, and thus lower weights than rare species from sites with many stations (e.g. a species found in one of four inventories has an occurrence of 0.25, while a species found in one of 20 inventories has an occurrence of 0.04). As a consequence, less sampled sites may be penalised by this correction, which may partly explain the lower overall indices of Morlaix Bay (N = 4), Morbihan Gulf (N = 6) and Etel Ria (N = 6). We believe that this correction is conservative, and the results may be improved in the future by targeting new sampling stations in these sites in order to achieve comparable sampling intensities with other sites. Overall, the pros derived from these corrections are that we estimate that our main findings regarding the discrepancy in cross-taxon rarity and the spatial distribution of rare assemblages are reliable. These results are the first step to refine the identification of important assemblages or habitats for conservation. Once this step is performed, further investigations can be undertaken on a subset of selected assemblages or habitats, which greatly reduce the costs of sampling in these habitats, inflated by the absence of surrogate phyla. To this aim, we recommend investigating species abundances that can be included in the Index of Relative Rarity [START_REF] Leroy | Rarity: Calculation of Rarity Indices for Species and Assemblages of Species[END_REF][START_REF] Leroy | Complementarity of rarity, specialisation and functional diversity metrics to assess community responses to environmental changes, using an example of spider communities in salt marshes[END_REF], in order to target assemblages containing species regionally rare, but locally abundant.

We focused only on animal taxa, and thus omitted the primary producers which are macroalgae (Heterokontophyta, Chlorophyta and Rhodophyta) despite the fact that they represent a significant proportion of the total biodiversity and biomass of subtidal rocky communities. However, our omission of macroalgae is at least partly offset by the indicator that has been specifically developed for the assessment of macroalgae under the European Water Framework Directive (Le [START_REF] Gal | Quality Index of Subtidal Macroalgae (QISubMac): A suitable tool for ecological quality status assessment under the scope of the European Water Framework Directive[END_REF]. Nevertheless, we recommend further sampling of macroalgae to include them in future studies investigating biodiversity surrogates and/or rarity patterns in subtidal assemblages of rocky shores in Brittany.

Biodiversity patterns of marine benthos are known to be scale-dependant for species richness and assemblage structure (e.g. [START_REF] Smale | Subtidal macroalgal richness, diversity and turnover, at multiple spatial scales, along the southwestern Australian coastline[END_REF], as well as for rarity [START_REF] Ellingsen | Rare species, habitat diversity and functional redundancy in marine benthos[END_REF]. This scale-dependency also holds for the effectiveness of biodiversity surrogates (e.g. [START_REF] Smale | Monitoring marine macroalgae: The influence of spatial scale on the usefulness of biodiversity surrogates[END_REF][START_REF] Sutcliffe | Biological surrogacy in tropical seabed assemblages fails[END_REF]. In our analysis, we detected this scale-dependency: rarity of Cnidaria was correlated to rarity of Porifera between sites but not within sites. Given that evidence regarding scale-dependency in biodiversity patterns of marine benthos and the probable role of environmental variability in explaining spatial distribution of rarity, we recommend that scaledependency should be tested specifically in future studies investigating biodiversity surrogates and/or rarity patterns in subtidal assemblages of rocky shores in Brittany. In addition, the database used to estimate species rarity in this paper is limited to Brittany, and thus does not permit us to distinguish between species which are intrinsically rare or are rare because they are at the edge of their geographical range. Therefore, future work to improve the explanation of the observed patterns should also investigate multi-scale patterns from a larger spatial perspective [START_REF] Leroy | Integrating multiple scales in rarity assessments of invertebrate taxa[END_REF].

Furthermore, we were not able to investigate temporal patterns because of the limited sampling costs. Yet, changes have been shown to occur in subtidal algal communities over the temporal span of our database samplings [START_REF] Gallon | Twenty years of observed and predicted changes in subtidal red seaweed assemblages along a biogeographical transition zone: inferring potential causes from environmental data[END_REF]. Among changes, it can be expected that some rare species may have become more common (such as warm-water species at their northern range limits), and common species may have become rare (such as cold-water species close to their southern range limits). Nevertheless, since (i) these changes have been rather slowly occurring, (ii) western and north-western Brittany have been identified as the most stable zones for other phyla (red seaweeds) [START_REF] Gallon | Twenty years of observed and predicted changes in subtidal red seaweed assemblages along a biogeographical transition zone: inferring potential causes from environmental data[END_REF] and (iii) our metric is averaged at the community level, we assume that these changes would not have major impacts on our findings. We nonetheless recognize that, given the accelerating environmental changes, the temporal aspect is a major challenge for such understudied and costly-to-sample habitats, and this challenge should be rapidly tackled in future studies.

Concluding remarks

Our initial aim was to identify surrogate phyla to indicate rarity of other phyla, to be used for the monitoring of subtidal assemblages of rocky shores in Brittany. Unfortunately, such surrogates could not be identified because of the observed discrepancy in rarity patterns among phyla. This finding has an important implication for the conservation of subtidal rocky habitats as it means that there is no shortcut to monitor their rarity: working on a subset of phyla implies a biased evaluation of biodiversity. Nevertheless, the method proposed here allowed us to successfully evaluate the rarity of sessile and low-mobility species as one facet of the conservation needs of subtidal rocky habitats. This method is thus a working indicator for the implementation of European directives as the MSFD, WFD and Habitat Directive.

Figure 1 .

 1 Figure 1. Box-and-whisker plots of A. species occurrences and B. rarity weights for each phylum. Each point stands for a species. Points are transparent; hence, the overlap of several points results in darker points. Box-and-whiskers represent the non-outlier range of values with horizontal bars representing from bottom to top: first quartile (bold horizontal segments in blue), median (bold horizontal segments), third quartile. For each phylum, rarity weights were calculated with a cut-off point corresponding to the first quartile of species occurrence.

Figure 2 .

 2 Figure 2. Box-and-whisker plots of phylum-specific Indices of Relative Rarity (IRR) of the 133 inventories of our database. Box-and-whiskers represent the non-outlier range of values with horizontal bars representing from bottom to top: first quartile, median, third quartile.

Figure 3 .

 3 Figure 3. Cross-taxon congruence of Indices of Relative Rarity (IRR) among the eight studied phyla. The lower half shows the scatter plots of IRR ranks of inventories among phyla. The upper half shows values of Spearman's rank correlation coefficients between

Figure 4 .

 4 Figure 4. Box-and-whisker plots of multi-phyla Indices of Relative Rarity (IRR) based on phylum-specific rarity cut-off points for the 133 inventories of our database. Box-andwhiskers represent the non-outlier range of values with horizontal bars representing from bottom to top: first quartile, median, third quartile. The dashed line represents the

Figure 5 .

 5 Figure 5. Spatial distribution of Indices of Relative Rarity (IRR) of sampled assemblages around Brittany. IRR were coloured according to four equal classes based on quartiles of all IRR values: (i) Lowest 25%: sites with IRR below the first quartile; (ii) Low-mid 25%: sites with IRR above the first quartile and below the median; (iii) Upper-mid: sites with IRR above the median and below the third quartile; (iv) Top 25%: sites with IRR above the third quartile.

Table 1 .

 1 Sampling completeness of the different sites of the database. Averaged estimated richness is the average richness based on three estimators: Chao2, ICE and Jack1 (all values in Appendix A). Completeness index: species richness divided by estimated species richness. Sampling intensity: number of inventories.

		Species richness	Average estimated richness	Completeness index	Sampling intensity
	Database	540	574	0.94	137
	Rance estuary	204	343	0.60	4
	St Malo bay	309	354	0.87	10
	Granite rose coast 310	341	0.91	16
	Sept Îles archipelago 355	409	0.87	12
	Lannion bay	323	363	0.89	13
	Morlaix bay	246	304	0.81	4
	Brest sea	256	282	0.91	11
	Brest bay	275	320	0.86	12
	Crozon peninsula	168	194	0.87	13
	Cape Sizun	218	263	0.83	10
	Ushant island	217	253	0.86	20
	Morbihan gulf	299	337	0.89	6
	Etel ria	265	311	0.85	6
	Mean values	265	314	0.84	11

Table 2 .

 2 Completeness and species richness of each phylum of the database. Completeness is the ratio between observed and estimated species richness[START_REF] Soberón | Assessing completeness of biodiversity databases at different spatial scales[END_REF]. For each phylum, the completeness is based on three estimators (Chao2, ICE and Jack1; all values in Appendix A) and is averaged across all sites. Numbers in brackets are standard deviations.

							Echino-
		Annelida	Arthropoda Bryozoa	Chordata	Cnidaria	dermata	Mollusca	Porifera
	Total richness	32	37	88	61	102	27	88	101
	Average completeness 0.83 (0.13) 0.92 (0.06) 0.85 (0.08) 0.85 (0.12) 0.85 (0.05) 0.86 (0.14) 0.84 (0.11) 0.85 (0.06)
	Average observed
	species richness per
	inventory	4.7 (3.3)	14.3 (2.9)	20.6 (15.7) 14.2 (9.5)	19.5 (10.6) 7.5 (2.9)	21.1 (5.5)	15.7 (10.8)
	Average observed
	species richness per
	site	11.4 (6.3)	17.2 (8.5)	42.9 (22.7) 29.6 (15.5) 46.3 (21.6) 13.3 (6.1)	31.4 (17.1) 39.6 (21.1)

Acknowledgements

We sincerely thank Thibault Schvartz from CREOCEAN for his initial work on the compilation of the database. We also thank two anonymous reviewers for useful comments and suggestions on an earlier draft of the manuscript. We thank Bernard Hugueny for useful discussions on the results of our study. BL was funded by the Service du Patrimoine Naturel of the Museum National d'Histoire Naturelle.

Supplementary Material


Pedro Sanzana
email: ppsanzana@uc.cl

Jorge Gironás
email: jgironas@ing.puc.cl

Isabelle Braud

Flora Branger
email: flora.branger@irstea.fr

Fabrice Rodriguez
email: fabrice.rodriguez@ifsttar.fr

Ximena Vargas
email: xvargas@ing.uchile.cl

Nancy Hitschfeld

José-Francisco Muñoz
email: jfmunoz@ing.puc.cl

Sebastián Vicuña
email: svicuna@ing.puc.cl

Urban and peri-urban drainage network representation using Geo-PUMMA in the Mercier (France) and El Guindo (Chile) catchments Représentation du réseau d'assainissement urbain et péri-urbain par l'utilisation de Geo-Pumma dans les bassins versants du Mercier (France) et de l'El Guindo (Chili)

Keywords: Peri-urban catchments, Drainage extraction, Computer-assisted mesh generation, GRASS-GIS

Les processus hydrologiques dans les zones urbaines et péri-urbaines sont significativement affectés par des éléments artificiels comme les canaux, les tuyaux, les rues et les réseaux d'eau pluviale. Nous présentons Geo-PUMMA, un ensemble d'outils SIG en Python pour GRASS qui représentent explicitement les différentes composantes des zones périurbaines par un maillage vectoriel irrégulier. Geo-PUMMA rassemble informations spatialisées telles que le cadastre, les types de sols, la géologie et des modèles numériques de terrain. Geo-PUMMA fournit des Unités de Réponse Hydrologique (HRUs) et des Eléments Hydrologiques Urbains (UHEs) et permet d'améliorer le maillage initial pour une meilleure représentation des chemins de l'eau. Geo-PUMMA permet également l'extraction de propriétés morphologiques du bassin telles que la fonction largeur, aire drainée ou d'imperméabilisation. Geo-PUMMA est appliqué à deux bassins péri-urbains: le Mercier (10% d'urbanisation) situé près de Lyon, France, et El Guindo (40% d'urbanisation) situé dans la région du piémont andin, Chili. La comparaison des fonctions largeur et aire drainée entre les maillages initiaux et finaux montre une amélioration de la représentation de réseau de drainage à petite et moyenne échelle (120-150 m pour El Guindo et 80-150 m pour le Mercier).

INTRODUCTION

Urban development and associated increase of soil imperviousness has significantly changed the urban hydrology cycle and modified the drainage network in non-natural areas. This change also affect peri-urban basins, characterized by a high level of heterogeneity, a mixture of natural, rural and urban areas (Braud et al., 2013). In this type of basins, typical elements or infrastructures can modify the direction of the superficial and sub-superficial flows, generating complex environments and increasing the types of connectivity between pervious and impervious areas that must be represented.

Several GIS tools have been developed to represent and visualize landscapes and extract information from them.. Most of these tools were developed to represent regional and medium size areas, and do not consider urban and peri-urban elements. Classical methodologies of drainage extraction and basin delineation use flow direction algorithms such as the D8 and Multiple Flow Directions (MFD), which consider Digital Elevation Models (DEM). These algorithms give good results in natural and non-flat areas at regional or median scales, but are not suitable for urban and peri-urban typical scales (<0.1-10 km 2) where surfece and subsurface infrastructure can modify the natural flow directions dramatically, and even change substantially the limits of urban and suburban basins. Thus to the best of our knowledge, a specific tool that allows creating a hydrological model mesh in urban and periurban catchments is still not available. In this work we develop Geo-PUMMA, a GIS tool that allows generating a good quality mesh for realistic urban and peri-urban hydrological connectivity and water pathways representation.

Geo-PUMMA uses a vectorial approach that preserves the principal physiographic units of urban and peri-urban basins, defined as polygons with irregular shape. Geo-PUMMA allows perparing urban and peri-urban terrain for semi-distributed or distributed hydrological models.Thus, it is possible to consider explicitly the infrastructures encountered in this complex environment (e.g. works of arts, retention infrastructures, aqueducts, piping, and streets or avenues). The representation of urban areas is based on Urban Hydrological Elements (UHE) (Rodríguez et al., 2008), while the rural and peri-urban areas are depicted using Hydrological Response Units (HRU) (Flüguel, 1995).

We apply Geo-PUMMA to the Mercier catchment (close to Lyon in France), and the Estero El Guindo catchment, a piedmont catchment located in Santiago, Chili. These two catchments differ in their landscape and climate conditions.

MATERIALS AND METHODS

Geo-PUMMA presentation

Geo-PUMMA was initially developed to prepare the geospatial information for the PUMMA hydrological model (Jankowfsky et al., 2014), but can also be used for any distributed hydrological model based on irregular modelling units. Thus previous work was put together in a set of tools programmed in Python language (python.org) using the GRASS platform (grass.osgeo.org) and QGIS (www.qgis.org). A methodology is proposed for helping in the spatial pre-processing of the HRUs and UHEs

The proposed methodology has four steps. Step 1 corresponds to data collection, digitalization and quality improvement of all the geospatial maps of relevance for hydrological processes in urban and peri-urban areas.

Step 2 corresponds to the representation of the urban area, in which all the UHEs are delineated and characterized using attributes such as average height, area, imperviousness, green area, and distance from the centroid of each unit to the closest sewer or street. In step 3 we improve the initial HRUs segmentation to address geometric constraints and insure accurate numerical modelling. To do this we use a triangulation and dissolution process based on geometrical indicators. In this step, the HRUs can be process to obtain more homogeneous properties such as slope, topographic index, aspect or other topographic parameter of interest that can be obtained from high resolution digital terrain models. Finally, in Step 4 the drainage network is extracted after accounting for the hydrological connectivity between the different units, using a recursive algorithm for identifying surface flow directions. This drainage network is composed of channelized and nonchannelized elements. Finally, we describe and assess the obtained drainage networks by means of geomorphologic properties.

Case studies, meshes generation and comparison

The Mercier catchment is part of Yzeron peri-urban watershed (150 km 2 ) located southwest of Lyon, France. The catchment has an area of 7 km 2 and is covered with forests, crops and urbananization (about 10%), although a significant portion is rural. The peri-urban Estero El Guindo basin (6.5 km 2) is located in the foothills of Santiago, Chili, in a location where urban sprawl is taking place quite rapidly. The natural area is covered by native vegetation (3.3 km 2 ) and the urban area (3.2 km 2 ) it is distributed along and close to the Estero El Guindo.

We applied Geo-PUMMA to the Mercier and El Guindo catchments. We generated three meshes to identify geometrical indices that allow getting the best representation drainage network, keeping the number of elements as low as possible to limit computing cost. The different maps were created using different values of 2 geometrical indices: the Convexity Index (CI) and Form Factor (FF). The first index measures the level of convexity of each element, and the second one measures how thin and needle-shaped elements are. The three meshes developed include the Initial Mesh (i.e. the mesh coming directly from the intersection of the land use, soil type, subbasins and geology layers), the Reference Mesh (a extremely dense mesh to represent the terrain as much as possible, which is created using high values of the geometric indices, i.e. CI = 0.975; FF = 0.50, and a maximum polygon area of 0.2 ha), and the Recommended Mesh (i.e. a mesh that compromise the initial and reference meshes, which is obtained using geometric indices that do not greatly increase the number of HRUs, but improve the water pathways representation to a large extent). Indeed, in this study we attempt to identify threshold values for the geometric indices to obtain this realistic yet simple representation of the terrain. After a process of trial and errors we obtained values of CI = 0.75 and FF = 0.2, and a maximum polygon area of 2 ha.

To analyze and compare the drainage network extracted from each mesh, the width-function, W(x), and area-function, A(x) were used (Shreve, 1969; Rodriguez-Iturbe y Rinaldo, 1997). These properties describe the flowpath organization and drainage contribution at different locations from the outlet x, which allows reducing a 2D drainage structure to a 1D mathematical function. These functions and the corresponding Cross Spectral Density (CPSD) functions were used to compare the results obtained using the various meshes.

RESULTS AND DISCUSSION

Figure 1 presents the results of the application of Geo-PUMMA to the two basins: the Mercier (left) and El Guido (rigth) in terms of model mesh and drainage extraction. The figure presents the three drainage networks corresponding to the three kinds of meshes described above. For both catchments, major changes are observed in the upper part of the catchments, where natural (resp. forest) land use were present and were initially digitized with larger polygons. These are the areas which were further segmented to get more representative water pathways. The mesh improvement in those areas at the largest distance from the outlet is visible in Figure 2 (top) with a large increase of NS from the initial to the recommended mesh (except for the area function of El Guindo). The increase in NS is larger for the width function than for the area function. The CPSD (Figure 2, bottom) applied to the widthfunction allows comparing the spatial similarity between the reference width-function and area-function versus the recommended and initial drainage network. Figure 2 shows improvement from 80 m to 60 m for CPSD of W(x) of the Mercier catchment and from 120 m to 60 m for CPSD of W(x) of El Guindo catchment. In the case of the area-functions, Figure 2 shows improvement from 150 m to 55 m for CPSD of A(x) of the Mercier and from 150 m to 100 m for CPSD A(x) of El Guindo. For both functions, NS is greatly improved when moving from the initial to the recommended mesh, expect for the area function of El Guindo, which requires further analysis. Generally, we find that the performance of the recommended mesh is satisfactory in terms of water pathways representation.

To conclude, the results show that Geo-PUMMA allowed getting the drainage network considering the main physiographic units of the natural and urban. The output layers and their corresponding databases provide useful information for any distributed hydrological model that requires a detailed representation of the conditions of land use and the hydrological connectivity. Our results show that it is not necessary to use the reference mesh as input for modeling, and that the recommended mesh can be used safely without increasing heavily the final number of HRUs to represent the topography and to obtain a realistic drainage network. The current version of Geo-PUMMA was carried out in GRASS 6.4 in a virtual machine with Ubuntu 14 (64b). The use of Geo-PUMMA and its methodology is recommended to users with spatial and hydrological modelling knowledge, which allow taking decisions when making simplifications on the HRU.

LIST OF REFERENCES

Braud, I., Fletcher, T.D., Andrieu, H., 2013. Preface: Hydrology of peri-urban catchments: processes and modelling, Journal of Hydrology, 485, 1-4. Flügel, W.A., 1995. Delineating hydrological response units by geographical information system analyses for regional hydrological modeling using PRMS/MMS in the drainage basin of the River Bröl, Germany. Hydrological Processes 9, 423-436. Jankowfsky, S., Branger, F., Braud, I., Rodriguez, F., Debionne, S., Viallet, P., 2014. Assessing anthropogenic influence on the hydrology of small peri-urban catchments: development of the object-oriented PUMMA model by integrating urban and rural hydrological models. Journal of Hydrology, 517, 1056-1071. Rodriguez, F., Andrieu, H., and Morena, F. (2008). A distributed hydrological model for urbanized areas -model development and application to case studies. Journal of Hydrology, 351:268-287. Rodriguez-Iturbe I., Rinaldo A., (1997). Fractal River Basins: Chance and Self-organization. Cambridge University Press, New York, pp. 1-6. Shreve R.L., (1966). Statistical Law of Stream Numbers, J. Geology, Vol (74), pp. 17-37.

Figure 1 ,

 1 Figure 1, Initial, Reference and Recommended Mesh for Mercier (a) and El Guindo (c). Drainage Network from Initial (b), Recommended (e) and Reference (f) Mesh for Mercier and Drainage Network from Initial (d), Recommended (g) and Reference (h) Mesh for El Guindo.

Figure 2 :

 2 Figure 2: (a) Width-Function, W(x), and (b) Area-Function, A(x), of three meshes of Mercier. (ca) Width-Function, W(x), and (d) Area-Function, A(x), of three meshes of El Guindo. In both drainage networks came from initial mesh (segmented gray line), recommended mesh (red) and reference (black). (e) CPSD of Width-Function and (f) CPSD of Area-Function of Mercier. (g) CPSD of Width-Function and (h) CPSD of Area-Function of El Guindo. In both cases, comparing reference vs reference (Py0y0, black line), reference vs initial (Py0y1, gray segmented line) and reference vs recommended (Py0y2, red line). At six figures, NS coefficients are indicated.

Acknowledgment to Projects MAPA (IDRC 107081-001), ECOS-CONICYT C14U02, and FONDECYT 1131131.