Alienor Jeliazkov
email: alienor.jeliazkov@gmail.com

Romain Lorrillière

Aurélien Besnard

Josette Garnier

Marie Silvestre

François Chiron

| Aurélien Besnard

Cross-scale effects of structural and functional connectivity in pond networks on amphibian distribution in agricultural landscapes

| INTRODUC TI ON

One of the main challenges in biodiversity research is better predicting species distribution in the face of unprecedented environmental changes and habitat loss. Improving our models requires going beyond the simple habitat models [START_REF] Guisan | Predicting species distribution: Offering more than simple habitat models[END_REF] and explicitly integrating ecological connectivity, as this is a crucial aspect of population persistence in the landscape [START_REF] Bennett | Linkages in the landscape: The role of corridors and connectivity in wildlife conservation[END_REF][START_REF] Hilty | Corridor ecology: The science and practice of linking landscapes for biodiversity conservation[END_REF].

Understanding the role of connectivity in species distribution is essential, but it remains a challenging task because it requires encompassing the diversity of scales at which processes shape species distribution [START_REF] Levin | The problem of pattern and scale in ecology[END_REF][START_REF] Noss | Landscape connectivity: Different functions at different scales[END_REF]. Environmental control, metapopulation dynamics, and seasonal migration are all key processes with different and, sometimes, interacting operating scales.

Moreover, connectivity can be characterised and modelled according to various paradigms [START_REF] Calabrese | A comparison-shopper's guide to connectivity metrics[END_REF][START_REF] Fischer | Landscape modification and habitat fragmentation: A synthesis[END_REF]. Landscape connectivity sensu lato was defined as "the degree to which the landscape facilitates or impedes individual movement between resource patches" [START_REF] Taylor | Connectivity is a vital element of landscape structure[END_REF], with certain differences between structural and functional connectivity [START_REF] Tischendorf | On the usage and measurement of landscape connectivity[END_REF][START_REF] Uezu | Effects of structural and functional connectivity and patch size on the abundance of seven Atlantic Forest bird species[END_REF]. Structural connectivity refers to the physical structure of the landscape (i.e. shape, number of resource patches, and distance between them) and can be independent of any trait of the organism of interest. On the contrary, functional connectivity accounts for the ecological traits of the organism [START_REF] Tischendorf | On the usage and measurement of landscape connectivity[END_REF]: maximum distance of dispersal, energetic cost of displacement across different types of land covers, mortality risks associated with crossing hostile habitats, etc. Based on speciesspecific characteristics, functional connectivity often appears as a better proxy of actual ecological connectivity between habitats that results from biological processes [START_REF] Fischer | Landscape modification and habitat fragmentation: A synthesis[END_REF].

In the present study, we focus on amphibians in pond ecosystems. Several studies have demonstrated the importance of pond connectivity for species and population persistence (e.g. [START_REF] Fortuna | Spatial network structure and amphibian persistence in stochastic environments[END_REF][START_REF] Ribeiro | The pond network: Can structural connectivity reflect on (amphibian) biodiversity patterns?[END_REF]. As ground-dwelling animals with inter-seasonal mobility, due to the biphasic life cycle of some species (i.e. aquatic and terrestrial), amphibians are dependent on both pond conditions and habitat connectivity.

In addition, some amphibians function according to metapopulation dynamics and have strong breeding site fidelity (e.g. [START_REF] Smith | Dispersal and the metapopulation paradigm in amphibian ecology and conservation: Are all amphibian populations metapopulations?[END_REF] and high species turnover (e.g. [START_REF] Trenham | Regional dynamics of wetland-breeding frogs and toads: Turnover and synchrony[END_REF]. These traits indicate that amphibian populations can be sensitive to both local and regional environmental conditions [START_REF] Skelly | Distributions of pond-breeding anurans: An overview of mechanisms[END_REF][START_REF] Werner | Turnover in an amphibian metacommunity: The role of local and regional factors[END_REF] and to the spatial configuration of the ponds [START_REF] Cushman | Effects of habitat loss and fragmentation on amphibians: A review and prospectus[END_REF].

There have been different ways to account for connectivity when studying amphibian distribution. Pond connectivity is often estimated through buffer approaches by calculating metrics such as the number of ponds, or the presence of barriers in circular zones around a focal pond [START_REF] Pellet | Multiscale determinants of tree frog (Hyla arborea L.) calling ponds in western Switzerland[END_REF][START_REF] Simon | Spatial correlates of amphibian use of constructed wetlands in an urban landscape[END_REF][START_REF] Zanini | Landscape effects on anuran pond occupancy in an agricultural countryside: Barrierbased buffers predict distributions better than circular buffers[END_REF]. Based on this approach, [START_REF] Prugh | An evaluation of patch connectivity measures[END_REF] showed that the number of adjacent occupied habitats, as potential sources of propagules, better predicts amphibian occupancy and colonisation than the simple number of suitable habitats. Some studies also assessed species distribution according to network-centred approaches. A network is considered as a group of inter-dependent ponds that can be related to each other via spatial and temporal processes, and exchanges of information. The explicit description of pond network connectivity can be either structural [START_REF] Fortuna | Spatial network structure and amphibian persistence in stochastic environments[END_REF][START_REF] Ribeiro | The pond network: Can structural connectivity reflect on (amphibian) biodiversity patterns?[END_REF] or functional [START_REF] Joly | Habitat fragmentation and amphibian conservation: Building a tool for assessing landscape matrix connectivity[END_REF][START_REF] Ray | Modeling spatial distribution of amphibian populations: A GIS approach based on habitat matrix permeability[END_REF]. Comparing network-centred and buffer approaches, [START_REF] Ray | Modeling spatial distribution of amphibian populations: A GIS approach based on habitat matrix permeability[END_REF] showed that the functional network approach performed better than the structural buffer approach in predicting amphibian occurrence. The former takes into account the energy costs associated with movement through particular types of habitats and, thus, the area and shape of potential migration zones, whereas the latter assumes a circular and homogeneously permeable migration zone. However, these comparative analyses do not take into account the local condition effects, although some studies have demonstrated the additive effects of pond characteristics and network variables [START_REF] Decout | Integrative approach for landscape-based graph connectivity analysis: A case study with the common frog (Rana temporaria) in human-dominated landscapes[END_REF][START_REF] Ribeiro | The pond network: Can structural connectivity reflect on (amphibian) biodiversity patterns?[END_REF], especially using valuable advances based on graph theory [START_REF] Clauzel | Impact assessment of a high-speed railway line on species distribution: Application to the European tree frog (Hyla arborea) in Franche-Comté[END_REF].

Another important aspect to consider when investigating the potential processes underlying a species' distribution, which is subject to different scales of influence, is the cross-scale interactions. The idea of cross-scale interactions [START_REF] Peters | Crossscale interactions and changing pattern-process relationships: Consequences for system dynamics[END_REF] dynamics. Both structural and functional connectivity should be considered as complementary features of pond connectivity when dispersal traits are likely to be less limiting. From a management perspective, our results encourage conservation plans to combine functional network-centred strategies with the preservation of local conditions. We give guidelines as to the features of such networks to help maintain amphibian populations in agricultural landscapes.

K E Y W O R D S
agricultural areas, amphibian conservation, landscape connectivity, pond networks, scales [START_REF] Schneider | The rise of the concept of scale in ecology: The concept of scale is evolving from verbal expression to quantitative expression[END_REF]) is a consequence of the hierarchy theory [START_REF] Allen | Hierarchy: Perspectives for ecological complexity[END_REF]. A relationship at one given level (whether spatial, temporal, or organisational) can be constrained by the characteristics of the upper level and can gain new properties from the combination of lower-level characteristics (O 'Neill, Johnson, & King, 1989). In particular, processes operating at different spatial scales can interact with each other and generate particular biological patterns on a given scale [START_REF] Lawler | A variance-decomposition approach to investigating multiscale habitat associations[END_REF]. Analysing cross-scale interactions thus allows one to link patterns emerging at different scales and to relate these patterns to potential processes [START_REF] Peters | Crossscale interactions and changing pattern-process relationships: Consequences for system dynamics[END_REF] (see an illustration of cross-scale interactions in Supporting Information Appendix S1).

Some studies investigated the interaction between pond-and landscape-level effects on amphibians, as the alteration of local conditions could, for instance, limit the positive effect of connectivity by impeding the settlement of colonisers [START_REF] Joly | Habitat fragmentation and amphibian conservation: Building a tool for assessing landscape matrix connectivity[END_REF]. Complex interactions between local habitat selection (e.g. for reproduction, growth), seasonal migration through the landscape, and metapopulation processes (dispersal) were shown to have important effects on both species distribution and community patterns [START_REF] Ficetola | Landscape-stream interactions and habitat conservation for amphibians[END_REF][START_REF] Werner | Cross-scale interactions and the distributionabundance relationship[END_REF]. Several studies have investigated the additive (e.g. [START_REF] Denoël | Conservation of newt guilds in an agricultural landscape of Belgium: The importance of aquatic and terrestrial habitats[END_REF] and interactive effects [START_REF] Ficetola | Landscape-stream interactions and habitat conservation for amphibians[END_REF][START_REF] Guerry | Amphibian distributions in a landscape of forests and agriculture: An examination of landscape composition and configuration[END_REF] of pond and landscape variables on amphibian populations. However, no framework has been available yet to investigate the interaction effects (cross-scale effects) between pond local conditions and functional and structural connectivity at pond and network scales together in a single study.

We propose here a hierarchical analytical framework based on a nested structure with three spatial scale levels [START_REF] Schaldach | An integrated approach to modelling land-use change on continental and global scales[END_REF], later on referred as to scales: pond, pond in network, and network. This nested structure is used to build the predictors of species distribution and approaches different scales of influence, while the response scale is fixed; the pond scale. This original approach allows us to describe the habitat and connectivity features, and to test their relation with species distribution, as well as all their within-scale and cross-scale interactions [START_REF] Peters | Crossscale interactions and changing pattern-process relationships: Consequences for system dynamics[END_REF].

In an agricultural zone, we estimated the occurrence probability of three amphibian species (Bufo bufo, Rana dalmatina, and Lissotriton vulgaris), which have different habitat requirements, distributional patterns, and dispersal abilities. (1) We tested the link between amphibian distribution, and structural and functional connectivity metrics to assess their relative influence. (2) We assessed the scaledependence and cross-scale effects of the environment on amphibian distribution. (3) For every species, we identified networks that could be considered as relevant management units and ranked their connectivity features in order of importance, which would consequently constitute action targets with regard to amphibian conservation in an agricultural landscape.

Our working hypotheses are as follows:

(1a) Generally, functional measures of network connectivity are expected to explain amphibian species distribution better than structural ones because the former take into account the species' perception of the landscape, the availability of potential sources of resources, shelters (e.g. for hibernation; [START_REF] Neill | Hibernation of amphibians and reptiles in Richmond County, Georgia[END_REF], and the biological risks involved in migration/dispersal processes (e.g. [START_REF] Janin | Use of stress-hormone levels and habitat selection to assess functional connectivity of a landscape for an amphibian[END_REF]. (1b) Assuming that amphibian population can benefit from the nearby presence of potential sources of propagules, we generally expect a positive effect of pond connectivity to any network (e.g. pond centrality) on the species' probability of occurrence. (1c) As a corollary, we expect network expansion, and network resistance to have negative effects, by making the potential sources further apart, and by altering the connectivity between these sources, respectively.

(2) We already showed the importance of local conditions (i.e. fish presence, water quality, and neighbouring proportion of woodland cover) for these specific species in this study area [START_REF] Jeliazkov | Level-dependence of the relationships between amphibian biodiversity and environment in pond systems within an intensive agricultural landscape[END_REF]. We here focus on the relative importance of network and pond-in-network characteristics in addition to, and in interaction with local conditions. We expected that the influences from the pond, pond-in-network, and network scales could interact by amplifying or suppressing each other. (2a) We expect local conditions to interact with connectivity: e.g. a pond with poor local conditions, even if well connected, will not allow the species to settle [START_REF] Joly | Habitat fragmentation and amphibian conservation: Building a tool for assessing landscape matrix connectivity[END_REF]. (2b) We also expect the network-scale variables to interact with each other; e.g. the spatial expansion of the network could have a positive effect in the case of low network resistance, because this would provide a higher number of available and reachable ponds but this effect could become negative if the resistance is too high.

(3) Because of species-specific dispersal abilities, we expect contrasting responses depending on the species, with the best dispersers being less limited by connectivity compared to poor dispersers.

We expect different thresholds to exist for the species' responses to network variables. In landscapes where the suitable habitats are initially scarce and fragmented, the positive relationships with connectivity metrics may only be detected for the highest values of habitat availability or accessibility [START_REF] Metzger | The structural connectivity threshold: An hypothesis in conservation biology at the landscape scale[END_REF].

| ME THODS

| Study area and pond survey

The study area (c. 430 km 2 ) is located in the French Seine-et-Marne department (eastern Ile-de-France region; 48.6°N/3.2°E). This site is representative of an intensive agricultural landscape, with 80% of the area covered by cereal crops.

In this study area, we located ponds based on the French national BD TOPO ® v2 (IGN, 2008) and two complementary regional pond censuses (2009-2012) (Département de Seine-et-Marne & SNPN, 2012) (spatial precision 5 m).

We selected 160 out of 300 ponds in the study area on the basis of two criteria. First, we selected permanent ponds to avoid major differences in ecological functioning, which is linked to the hydrological regime. Second, we chose ponds within clusters of different pond densities in which inter-pond distances did not exceed 2,000 m. This distance allowed us to: (1) to remain in the range of the known maximum migration distances of most of the amphibian species under study [START_REF] Jeliazkov | Level-dependence of the relationships between amphibian biodiversity and environment in pond systems within an intensive agricultural landscape[END_REF]; and (2) have a sufficient number of distinct clusters for analyses (n = 15), and a reasonable gradient of cluster densities, i.e. number of ponds per cluster (n ∊ {1, 2,3,4,5,6, 7,12,13,20,27,50}). These clusters are not the networks themselves but were used as initial patterns to guide the sampling. Very few interstitial ponds were missed or unsampled in our area. We omitted temporary ponds as well as those destroyed during the survey.

We acknowledge that temporary ponds potentially act as important stepping stone habitats for some species and their presence could influence the analysis of amphibian distribution. In our case, however, most of the ponds were historically artificial, and there were no temporary ponds within the sampling clusters (BD TOPO ® v2 IGN, 2008) and potentially only very few in the vicinity of clusters (Département de Seine-et-Marne & SNPN, 2012;[START_REF] Jeliazkov | Level-dependence of the relationships between amphibian biodiversity and environment in pond systems within an intensive agricultural landscape[END_REF]. Thus, their effect is probably negligible compared to the effects we investigated.

We assessed the occurrence of three amphibian species-B. bufo, R. dalmatina, L. vulgaris-at the 160 sites in 2011 and 2012 using a standardised protocol (52 ponds in 2011, 65 in 2012, and 35 in both years). Since there was no substantial year effect in the speciesoccurrence data [START_REF] Jeliazkov | Level-dependence of the relationships between amphibian biodiversity and environment in pond systems within an intensive agricultural landscape[END_REF], we merged the data of the two annual campaigns into a single visit dataset for the analyses [START_REF] Chase | Connectivity, scale-dependence, and the productivity-diversity relationship[END_REF].

We conducted two fortnight nocturnal sessions per year during the main reproductive seasons (March and June). We combined three complementary methods to detect, identify, and count individuals of the target taxa at different stages, while minimising habitat destruction: records of calling males, observations using lamps along the pond banks, and dip-netting. See [START_REF] Jeliazkov | Level-dependence of the relationships between amphibian biodiversity and environment in pond systems within an intensive agricultural landscape[END_REF] for more details.

We also recorded the presence/absence of potentially predatory fish based on visual observations, dip-netting, and/or information provided by pond owners.

Furthermore, the physical-chemical characteristics of the water were measured four times (twice per year): dissolved oxygen, water temperature, pH, concentrations of ammonium, ammonia, nitrates, nitrites, silica, phosphates, and chlorophyll a (see [START_REF] Jeliazkov | Level-dependence of the relationships between amphibian biodiversity and environment in pond systems within an intensive agricultural landscape[END_REF]. Based on these parameters, we calculated an average water quality index (WQI) ranging from 0% (very poor quality) to 100% (very good quality). The WQI was adapted from [START_REF] Pesce | Use of water quality indices to verify the impact of Cordoba City (Argentina) on Suquia River[END_REF][START_REF] Sánchez | Use of the water quality index and dissolved oxygen deficit as simple indicators of watersheds pollution[END_REF] as it has proven to be a good indicator of water quality with regard to amphibian response.

The availability of terrestrial habitats around the ponds was represented by the proportion of wood cover within a 200-m-radius buffer zone surrounding the ponds (for more information on the methods, see [START_REF] Jeliazkov | Level-dependence of the relationships between amphibian biodiversity and environment in pond systems within an intensive agricultural landscape[END_REF]. These variables were used as pond-scale predictors of amphibian occurrence (see next sections).

| Study species

We studied two anuran and one urodele species that have different habitat requirements, distributional patterns, migration behaviours, and dispersal abilities (for a literature-based synthesis, Supporting Information Appendix S2). The three species have the following general sequence of seasonal events in common; prenuptial migration of the adults from summering sites to the breeding pond; reproduction in the breeding pond; egg-laying followed by the development of tadpoles/larvae; postnuptial migration of adults to the wintering sites; dispersal of juveniles; and wintering period. However, they have different phenology and use different types of habitat to complete these stages.

Bufo bufo-the common toad-generally prefers cool, wet, nutrientrich and wooded environments (ACEMAV coll., 2003) but it can be found in a variety of habitats [START_REF] Arnold | Le guide herpéto -Amphibiens et reptiles d'Europe[END_REF]. The wintering sites are generally close to the breeding site (<500 m, up to 1,000 m; ACEMAV coll., 2003). After hibernation, the adults massively migrate to the reproduction zone that is often-but not necessarily-the same as the previous year [START_REF] Reading | Breeding pond fidelity in the common toad, Bufo bufo[END_REF]. The reproduction results in c. 3,000-8,000 eggs per female per season. In summer, juveniles leave the water and shelter in ground cracks, or under the vegetation around the pond shore before being able to disperse across terrestrial habitats. After reproduction, the adults migrate to summering sites up to 1,000 m away from the breeding site [START_REF] Sztatecsny | Into thin air: Vertical migration, body condition, and quality of terrestrial habitats of alpine common toads, Bufo bufo[END_REF] but can travel up to 3 km [START_REF] Smith | Dispersal and the metapopulation paradigm in amphibian ecology and conservation: Are all amphibian populations metapopulations?[END_REF] Supporting Information Appendix S2). They have few predators (snakes, herons, otters) and specific defence mechanisms against predation from fish [START_REF] Gunzburger | Critical literature review of the evidence for unpalatability of amphibian eggs and larvae[END_REF][START_REF] Kats | Antipredator defenses and the persistence of amphibian larvae with fishes[END_REF].

Rana dalmatina-the agile frog-prefers woodland habitats, especially dense deciduous forest, although it can also be found in woodlots and some grasslands. It can stay in dry areas outside the reproduction season [START_REF] Arnold | Le guide herpéto -Amphibiens et reptiles d'Europe[END_REF]. It coexists with other amphibian species but usually avoids fishponds. In summer, the adults have distinct vital ranges of about tens of meters square in wooded or grassland areas.

The distance between the vital range and the reproduction site reaches 1 km (ACEMAV coll., 2003). The reproduction period is usually short and results in 450-1,800 eggs per female per season. Rana dalmatina migrates to terrestrial habitat in autumn and hibernates in the ground close to the reproduction site (Supporting Information Appendix S2).

Lissotriton vulgaris-the smooth newt-thrives in open habitats [START_REF] Denoël | Landscape-level thresholds and newt conservation[END_REF] such as agricultural areas, gardens, and field edges, but there is no real consensus about its actual relation to woodlands, whether positive or negative [START_REF] Denoël | Similar local and landscape processes affect both a common and a rare newt species[END_REF]; Supporting Information Appendix S2). It may have a high fidelity to the reproduction site [START_REF] Smith | Dispersal and the metapopulation paradigm in amphibian ecology and conservation: Are all amphibian populations metapopulations?[END_REF]. It prefers small and shallow ponds with vegetation and ditches. The reproduction season can end with a postnuptial migration or with a longer stay in the water, sometimes until December (ACEMAV coll., 2003).

The egg-laying is much less productive than the two other species with c. 200-300 eggs per female per season. The hibernation of the adults usually takes place in the ground.

| Definition of structural and functional networks

To accurately characterise the habitat connectivity between the sampled ponds and delimit the networks, we combined three geographic information system (GIS) regional databases: MOS (IAU IDF, 2008) providing detailed information about urban infrastructures and arable land, Ecomos (IAU IDF, 2005) focusing on semi-natural areas such as forest and grassland, and Ecoline (IAU IDF, 2012) surveying all semi-natural features of the landscape such as single trees, hedgerows, groves, and field margins at a very fine spatial resolution (1:2,500) and precision (20 cm) (Supporting Information Appendix S3).

To create network maps, we first converted all the vector layers to raster format with a 3 × 3 m pixel resolution and then merged them to a single raster layer referred to as landscape (Supporting Information Appendix S4). The grain size used for connectivity analyses needs to be appropriately chosen [START_REF] Baguette | Landscape connectivity and animal behavior: Functional grain as a key determinant for dispersal[END_REF][START_REF] Galpern | Finding the functional grain: Comparing methods for scaling resistance surfaces[END_REF]. In our case, the grain of 3 × 3 m was a good compromise between the species' perception abilities [START_REF] Rothermel | Migratory success of juveniles: A potential constraint on connectivity for pond-breeding amphibians[END_REF] and the high precision level (both in spatial and thematic resolutions) provided by the landscape database Ecoline, thus allowing a maximum amount of valuable information to be conserved.

All GIS treatments, including generating the networks, were performed with ArcGIS for Desktop9.3/9.4 ® (under Basic license, formerly ArcView), Spatial Analyst, and ETGeoWizard10 ® .

To assess the resistance of the landscape features, we attributed resistance coefficients to the main cover types based on previous studies dealing with the same three amphibian species. Following an extensive literature search focussing precisely on the three species under study, we identified 10 studies that provided direct, indirect, or expert estimates. Values for mainland covers such as woodlands, crops, grasslands, urban areas, and roads were taken from the following references: [START_REF] Cushman | Habitat fragmentation effects depend on complex interactions between population size and dispersal ability: Modeling influences of roads, agriculture and residential development across a range of life-history characteristics[END_REF][START_REF] Decout | Integrative approach for landscape-based graph connectivity analysis: A case study with the common frog (Rana temporaria) in human-dominated landscapes[END_REF][START_REF] Janin | Assessing landscape connectivity with calibrated cost-distance modelling: Predicting common toad distribution in a context of spreading agriculture[END_REF][START_REF] Joly | Habitat fragmentation and amphibian conservation: Building a tool for assessing landscape matrix connectivity[END_REF][START_REF] Lenhardt | An expert-based landscape permeability model for assessing the impact of agricultural management on amphibian migration[END_REF][START_REF] Patrick | Multi-scale habitat-resistance models for predicting road mortality "hotspots" for turtles and amphibians[END_REF][START_REF] Ray | Modeling spatial distribution of amphibian populations: A GIS approach based on habitat matrix permeability[END_REF][START_REF] Van Buskirk | Permeability of the landscape matrix between amphibian breeding sites[END_REF][START_REF] Vos | Matrix permeability of agricultural landscapes: An analysis of movements of the common frog (Rana temporaria)[END_REF][START_REF] Zetterberg | Making graph theory operational for landscape ecological assessments, planning, and design[END_REF]. Note that Van Buskirk (2012) was only used to inform the relative resistance of the different habitat types as it provided measures that were not directly comparable to ours. Because urodeles have usually proved to be worse migrators than anurans (e.g. [START_REF] Semlitsch | Are small, isolated wetlands expendable?[END_REF][START_REF] Smith | Dispersal and the metapopulation paradigm in amphibian ecology and conservation: Are all amphibian populations metapopulations?[END_REF], we amplified the resistance values of unfavourable habitats for the newt species [START_REF] Ray | Modeling spatial distribution of amphibian populations: A GIS approach based on habitat matrix permeability[END_REF]. Although we found resistance values for most areal elements, there was no information for the fine elements censed in Ecoline except for hedgerows [START_REF] Lenhardt | An expert-based landscape permeability model for assessing the impact of agricultural management on amphibian migration[END_REF]. Thus, based on a thorough reasoning about the potential biological risks of the species crossing various habitats, we classified these elements according to the cover structure (strata height and density), hierarchised the associated biological risks (i.e. predation, desiccation, poisoning, physical barrier, and direct mortality), and attributed a value of resistance accordingly (Table 1).

The landscape raster was then reclassified into a cost raster using resistance values. The higher the probability for a species to migrate through a cell, the lower the cost.

Using the cost distance algorithm from ArcGIS Spatial Analyst, we built a map of the least accumulated costs (weighted distances) from any pixel of our study area to the nearest pond. This allowed us to define potential migration areas (PMA) around the ponds that formed functional networks when merged together according to a set of species-specific maximal migration distances (MMD; see the next section). Then, to compute the least-cost paths between prospected ponds, we used the cost path algorithm and recorded the length and mean accumulated cost of the paths (Supporting Information Appendix S4).

| Analytical framework for comparing effects of structural and functional connectivity and assessing cross-scale effects

The structural and functional approaches were all constrained by the same set of species-specific MMD (ACEMAV coll., 2003;[START_REF] Semlitsch | Biological criteria for buffer zones around wetlands and riparian habitats for amphibians and reptiles[END_REF][START_REF] Smith | Dispersal and the metapopulation paradigm in amphibian ecology and conservation: Are all amphibian populations metapopulations?[END_REF]:

• 1,000, 1,500, and 2,000 m for B. bufo;

• 300, 500, and 1,000 m for R. dalmatina;

• 100, 200, and 800 m for L. vulgaris.

We used three different approaches to build networks and measure pond connectivity:

1. The structural approach (S) refers to potential pond connectivity [START_REF] Calabrese | A comparison-shopper's guide to connectivity metrics[END_REF] and was based on inter-pond Euclidean distances in the limit of the MMD (Figure 1a).

2.

The area-functional approach (AF) was based on the merging of the PMAs around the ponds resulting from the mapping of the landscape resistance (Figure 1b). PMAs are proxies of habitat connectivity for migration between ponds by adult individuals [START_REF] Ray | Modeling spatial distribution of amphibian populations: A GIS approach based on habitat matrix permeability[END_REF] and/or for dispersal. They also reflect habitat availability and quality around ponds, including potential sources of resources (e.g. [START_REF] Janin | Beyond occurrence: Body condition and stress hormone as integrative indicators of habitat availability and fragmentation in the common toad[END_REF] and shelters.

3.

The path-functional approach (PF) involved considering the physiological constraint of cost accumulation along the least-cost paths that linked the ponds to each other from among all the possible distances in the limit of the MMD (Figure 1c). In this case, a dispersing animal accumulates costs during its displacement in relation to potential biological risks encountered over the distance travelled in the limit of the MMD (e.g. [START_REF] Bartelt | Modeling amphibian energetics, habitat suitability, and movements of western toads, Anaxyrus (=Bufo) boreas, across present and future landscapes[END_REF]. In some cases, the path-functional networks can appear less limited than the area-functional ones in terms of the number of connected ponds, but this implies a unique path and a directional inter-pond movement. As the accumulated cost measures have been shown to be more realistic than the least-cost path measures for assessing habitat connectivity [START_REF] Etherington | Least-cost path length versus accumulated-cost as connectivity measures[END_REF], the path-functional measure of network connectivity was expected to better explain the species distribution compared to the structural measure and to be complementary with the area-functional measure. The resistance value of pond has been fixed at 1 to ensure that they are the main attractors in the friction maps, that is required to study the connectivity between ponds precisely.

We thus represented influences deriving from three scales: pond, pond-in-network, and network scales. We considered the following eight variables as potential predictors of amphibian occurrence probability: fish presence, water quality, proportion of woodland in the direct surrounding of the pond, distance to the nearest occupied pond, pond degree centrality (according to the three different connectivity approaches S, AF, and PF), network expansion (based on S), network resistance (based on AF), and network cumulated cost (based on PF) (see details in Table 2).

| Statistical analyses

The objectives of the analyses were to assess: ( 1 For all models, we first checked for multicollinearity within our set of explanatory variables. The three variables of degree centrality (PondCentral_S/AF/PF) were collinear. To ensure the comparability between models while solving the problem of collinearity, we kept the degree centrality calculated using the area-functional approach and excluded the two other measures in all models (Supporting Information Appendix S5a).

To analyse the relative and interactive effects of pond-, pondin-network-, and network-level variables on species distribution, we We assessed the explanatory power of the final models using the pseudo-R 2 calculated as (TotalDeviance -ResidualDeviance)/ TotalDeviance.

To ensure model parsimony, we performed model simplification and estimated the percentage of deviance explained by each variable from the total deviance of the model. This allowed us (a) to rank variables according to the relative strength of their association with the species occurrence probability. We also (b) tested and selected all pairwise interactions for their relative importance (according to the residual variance generated by the linear model). As the boosted regression tree supports non-linear relationships, a graphical visualisation of the estimated species' response across scales allowed us (c) to identify potential thresholds.

All statistical analyses were performed with the statistical program R 3.3.1 (R Core Team, 2016).

| RE SULTS

| Species distribution and network description

Species occurrence and PMAs are illustrated in Figure 2. In general, the network number varied from eight to 39 depending on the species, migration distance, and network type (i.e. structural, area-functional, and path-functional). A substantial number of ponds turned out to be isolated, especially in the case of the areafunctional network, which was generally amplified as the MMD decreased (Supporting Information Appendix S6).

| Effects of functional versus structural network connectivity measures on species occurrence

When focussing on the network scale variables, all three connectivity metrics were predictive of species occurrence. However, the functional measures of network connectivity better explained species occurrence than the structural measure. The mean resistance within the area-functional network better explained the occurrence of R. dalmatina and L. vulgaris (Figure 3b,c), whereas the mean accumulated cost within the path-functional network better explained the occurrence of B. bufo (Figure 3a).

The probability of occurrence of both B. bufo and R. dalmatina showed a negative relationship with the network-scale variables (Figure 4a,b). In particular, the thresholds in the fitted function of occurrence probability in response to the network expansion indicated that the response to the mean network expansion became critical from the values of 600 m, and 100 m for B. bufo, and R. dalmatina, respectively. Both species' occurrence probability also responded negatively to the mean network resistance, with a critical threshold value of average resistance at 30% (Figure 4a,b). By contrast, L. vulgaris' occurrence showed a positive relationship with the mean network resistance from a threshold value of 50% (Figure 4c).

However, these thresholds substantially varied depending on the MMD (Supporting Information Appendix S7).

| Scale-dependence and cross-scale effects of pond and network variables

On average, the pond variables explained a slightly higher deviance than the network variables in the cases of R. dalmatina and L. vulgaris, whereas pond and network variables were equally ranked in the case of B. bufo (Figure 3). However, this pattern largely depended on the TA B L E 2 Description of the pond-, pond-in-network-, and network-scale variables of which the main and interaction effects were tested for species occurrence. S, AF, and PF respectively refer to the structural, area-functional, and path-functional approaches of connectivity 3).

The WQI had a positive relationship with species occurrence, generally for a minimum value of 70%, corresponding to the average level of good water quality (Figure 4). The distance to the nearest pond occupied by the focal species (DistOccPond) showed a strong negative relationship for L. vulgaris, also notable for R. dalmatina (Figure 4).

Boosted regression tree analysis highlights the interaction effects among pond-, pond-in-network-, and network-level variables, although most of the interaction effects had a relatively weak contribution to the deviance explained by the different predictors (Table 4a). The pond variables especially interacted with the higher-scale variables, i.e. with network and pondin-network variables in 40% and 25% of cases, respectively.

Pond-scale variables also interacted with each other (20%) (Table 4b). In particular, we observed a sharper decrease in the occurrence of L. vulgaris in response to fish when the water quality was better (Figure 5a). Pond variables can also interact with pond-in-networkscale variables such as DistOccPond (Figure 5b,c); e.g. the more connected the focal pond is, the stronger the response of L. vulgaris to water quality (Figure 5b). The same observation can be made for R. dalmatina in the case of fish presence (Figure 5c). Interactions between pond-and network-scale variables showed how the positive response to connectivity can be moderately reduced by the negative effects of local conditions (Figure 5d,e). Finally, we found only one small interaction between network variables for R. dalmatina, which involved network expansion and network resistance (Figure 5f).

| D ISCUSS I ON

| Structural versus functional characteristics of connectivity in relation to amphibian distribution

At the network influence scale, we showed the overall power of functional connectivity variables in explaining amphibian species distribution in line with previous studies [START_REF] Joly | Habitat fragmentation and amphibian conservation: Building a tool for assessing landscape matrix connectivity[END_REF][START_REF] Ray | Modeling spatial distribution of amphibian populations: A GIS approach based on habitat matrix permeability[END_REF]. In particular, area-functional connectivity better explained the distribution of R. dalmatina and L. vulgaris compared to path-functional and structural connectivity. This was expected (hypothesis 1a) since this connectivity metric probably includes the effects of both aquatic and terrestrial habitat availability/quality (breeding sites, source of resources, shelters) and habitat connectivity (potential migration and dispersal routes). Another possible explanation is that the area-functional approach makes no particular assumption about the movement direction; it thus encompasses all possible paths an animal could take within the potential migration area. Although the movement directions for migration and dispersal were shown to be non-randomly distributed for several species [START_REF] Marty | Directionality of preand post-breeding migrations of a marbled newt population (Triturus marmoratus): Implications for buffer zone management[END_REF][START_REF] Sinsch | Movement ecology of amphibians: From individual migratory behaviour to spatially structured populations in heterogeneous landscapes[END_REF], the directional movement might be highly disrupted in areas disturbed by human activities such as intensive agriculture [START_REF] Mazerolle | Landscape resistance to frog movements[END_REF].

However, in the case of B. bufo at shorter migration distance, structural connectivity explained almost as much of the species occurrence variation as functional connectivity. This suggests that in this range of migration distance, both network resistance and network expansion might be relevant predictors of toad distribution.

Given the relatively better dispersal ability of B. bufo compared to the two other species (ACEMAV coll., 2003;[START_REF] Smith | Dispersal and the metapopulation paradigm in amphibian ecology and conservation: Are all amphibian populations metapopulations?[END_REF],

we may be able to detect both the effects of network quality and network size. This result suggests that at the network scale, the functional and structural features of the connectivity should not be treated exclusively, but rather in complementarity to get a more complete view of the ecological connectivity in a system.

At the pond-in-network influence scale, the distance to the nearest occupied pond remained a good predictor of species occurrence, in particular for L. vulgaris, and outperformed the functional pond centrality, contrary to our prediction in (1b). [START_REF] Moilanen | Simple connectivity measures in spatial ecology[END_REF] reviewed that the distance to potential source populations was an effective predictor of habitat use, colonisation, extinction, and population genetic differentiation in the particular case of amphibians in highly fragmented habitats. The sensitivity of amphibians regarding the distance to the nearest occupied pond is often attributed to metapopulation processes and individuals' movements between close ponds, but it can also be due to habitat availability for adults [START_REF] Van Buskirk | Local and landscape influence on amphibian occurrence and abundance[END_REF]. In our case, the metapopulation hypothesis cannot be rejected because the distance to the nearest occupied pond (considering ponds as potential sources of propagules)

showed greater effect than the degree centrality (considering ponds as potential available habitats in the surrounding, i.e. ponds that are not necessarily occupied by the species). However, to test this hypothesis in our framework, we would need to improve the parameterisation of our networks so that they better reflect the actual ecological requirements of this species.

| Cross-scale effects of pond and network variables on species distribution

| Main effects of pond and network variables across spatial scales

The pond-scale effects on both anurans' occurrence usually dominated the other scales. Habitat preference, including for reproduction, may thus be one major limiting process underlying their population persistence in this landscape (e.g. [START_REF] Van Buskirk | Local and landscape influence on amphibian occurrence and abundance[END_REF].

By contrast, L. vulgaris was especially responsive to pond isolation with a sharp negative relationship, suggesting that its distribution has a strong spatial auto-correlation and may rely on source-sink processes [START_REF] Harrison | Local extinction in a metapopulation context: An empirical evaluation[END_REF] or on rescue effects.

At the network scale, reduced connectivity had the expected negative effects on pond occupancy probability for both anurans (1c), but a positive effect on L. vulgaris in terms of both network resistance and accumulated cost for joining the focal pond. In a study about the links between landscape permeability and amphibian migration, [START_REF] Lenhardt | An expert-based landscape permeability model for assessing the impact of agricultural management on amphibian migration[END_REF] showed that woodland areas had a good habitat suitability for L. vulgaris and thus used low cost values for woodland areas. We used this study to set the resistance values of woodlots for this species. However, in other studies, L. vulgaris occurrence has proven to be rather favoured by open landscapes, especially grasslands, and even positively responded to the distance to the nearest woodland [START_REF] Denoël | Landscape-level thresholds and newt conservation[END_REF]. This could explain the positive link we find for network resistance, although there are still contradicting results about the effect of the distance to the nearest forest on this species (Table 2 in [START_REF] Denoël | Similar local and landscape processes affect both a common and a rare newt species[END_REF]. If woodlands constitute actual barriers for L. vulgaris, our evaluation according to which the tree-filled areas were favourable habitats might be reconsidered [START_REF] Garshelis | Delusions in habitat evaluation: Measuring use, selection, and importance: Measuring use, selection, and importance[END_REF]. The results for L. vulgaris are also difficult to interpret as they include many isolated ponds for which the measure The relative ranking of influences from pond, pond-in-network, and network scales varies depending on the species, which gives some insights about the specific differences in the dominating scales of influence. For instance, B. bufo is more responsive to network-scale influences whereas L. vulgaris is rather responsive to pond-in-networkscale influence. This could be explained by strong differences in their dispersal limitations or in their spatial orientation (e.g. [START_REF] Sinsch | Homeward orientation of displaced newts (Triturus cristatus, Lissotriton vulgaris) is restricted to the range of routine movements[END_REF]. We need further investigations to conclude about what are the actual processes that underlie these species' distribution.

TA B L E 3 Ranking of the three scales of influence (pond, P; pond-in-network, PN; network, N) on species occurrence according to the total amount of deviance explained by the variables from each tested scale, species, and maximal migration distance (MMD); in grey, the best models in terms of deviance explained

| Interaction effects of pond and network variables across spatial scales

The interaction effects were relatively small compared to the main effects. Nevertheless, our results support the existence of crossscale effects in amphibian-environment relationships [START_REF] Werner | Cross-scale interactions and the distributionabundance relationship[END_REF], especially between pond and network scales. As expected in (2a), the positive relationship with network connectivity was usually amplified in the presence of good water quality and the absence of fish. (Reciprocally, a poor water quality and the presence of fish would annihilate the positive relationship with connectivity by impeding the settlement of population-ponds would then become sinks.) Although the detection of interactions at the network scale was more arduous-probably due to an increase of complexity-we showed that the network expansion could amplify the negative effect of network resistance as predicted in (2b). Overall, these results highlight that the different processes shaping species distribution are likely to operate and interact with each other across a gradient of scales.

| Methodological insights and challenges for a cross-scale expertise of ecological connectivity

The overall performance of our models showed substantial sensitivity to the species maximum migration distances modelled. This sensitivity gives some indications about the relevance of the range of maximum distances used for each species. It generally suggested that for the anuran species, the distance of c. 1,000 m was relevant to capture the potential effects of both dispersal limitation (connectivity metrics), habitat suitability (local conditions), and possibly habitat complementarity [START_REF] Janin | Assessing landscape connectivity with calibrated cost-distance modelling: Predicting common toad distribution in a context of spreading agriculture[END_REF]. We may have tested ranges at the limit and overestimated the potential migration distances for B. bufo (e.g. [START_REF] Sztatecsny | Into thin air: Vertical migration, body condition, and quality of terrestrial habitats of alpine common toads, Bufo bufo[END_REF].

However, the ranges we tested remain overall in the same order of magnitude than the ones found in most of the radio-tracking studies [START_REF] Smith | Dispersal and the metapopulation paradigm in amphibian ecology and conservation: Are all amphibian populations metapopulations?[END_REF]. The fact that the models based on 300 m and 1,000 m performed equally well for R. dalmatina, also in accordance with [START_REF] Lemckert | Variations in anuran movements and habitat use: Implications for conservation[END_REF] greater distances than adult habitat use in some anuran species (e.g. [START_REF] Angelone | Where movement happens: Scale-dependent landscape effects on genetic differentiation in the European tree frog[END_REF].

Note that even though the migration distances were chosen within a reasonable range based on previous works (Supporting Information Appendix S2), the ranking of predictors also changed depending on these distances. This suggests that using only one distance would have led to biased or incomplete conclusions about the dominating scales of influence. Given our still scattered knowledge about movement patterns and distances for amphibian species, it seems reasonable to test as many distances as possible.

We used biological risk-based criteria related to vegetation stratum to classify landscape features and assign them resistance coefficients. This approach, despite deriving from expert knowledge, is more transparent, standardised, and evaluable than a classical, fuzzy expert-based approach, as it relies on a framework of explicit hypotheses along with the relevant literature. Other approaches directly estimated ecological connectivity based on biological data from either species distribution [START_REF] Janin | Assessing landscape connectivity with calibrated cost-distance modelling: Predicting common toad distribution in a context of spreading agriculture[END_REF], population genetic structure [START_REF] Stevens | Gene flow and functional connectivity in the natterjack toad[END_REF], or radio-tracking [START_REF] Smith | Dispersal and the metapopulation paradigm in amphibian ecology and conservation: Are all amphibian populations metapopulations?[END_REF]. When such data are available in good quality and with reasonable spatial coverage, it seems highly recommendable to use them for connectivity analyses, as they avoid the need to subjectively choose resistance values. Otherwise, the risk-based approach is a good alternative for guiding cost-effective, applied conservation (Supporting Information Appendix S8).

| Thresholds in local and connectivity effects: Implications for amphibian conservation in agricultural landscapes

The detection of cross-scale effects of the environment on species occurrence proves that our analysis design gives a reasonable picture of the intrinsic scales [START_REF] Wu | Concepts of scale and scaling[END_REF] at which processes may operate in shaping species distribution. Our approach integrates multi-level processes [START_REF] Peters | Crossscale interactions and changing pattern-process relationships: Consequences for system dynamics[END_REF]. It should be considered for more effective pond management [START_REF] Semlitsch | Critical elements for biologically based recovery plans of aquatic-breeding amphibians[END_REF] to design management units that have ecologically meaningful contours and spatial cohesion at the landscape scale.

As predicted in (3), we showed the existence of thresholds in the species occurrence probability in response to predictors from different scales. The existence of such threshold values for habitat connectivity is consistent with the principle of percolation theory [START_REF] Stauffer | Introduction to percolation theory[END_REF]. Critical connectivity thresholds emerge when the abundance and arrangement of favourable patches, as well as the displacement abilities of the species, become limiting.

The responses to local conditions such as water quality differed quantitatively among the three species, thus reflecting the differential habitat selection of these species. Connectivity results also highlighted potential optimal pond networks for these species in agricultural contexts. For B. bufo, the appropriate network would

| CON CLUS ION

Our work provides an original analytical framework for investigating the multifactorial relationships between amphibian distribution and environment, and highlighting cross-scale effects and nonlinearities in amphibian response to habitat quality and connectivity. We integrated in one single study all the aspects that have usually been addressed separately when studying amphibian distribution. We thus addressed the different scales of influence from pond to network scales, their relative importance and the way they interact, the buffer-based versus network-based and the structural versus functional approaches of connectivity.

We emphasised the relevance of network approaches as a complement to pond-centred approaches for amphibian conservation [START_REF] En Boothby | Pond conservation: Towards a delineation of pondscape[END_REF][START_REF] Roe | Heterogeneous wetland complexes, buffer zones, and travel corridors: Landscape management for freshwater reptiles[END_REF]. We encourage conservation plans to adopt functional network-centred strategies based on potential migration areas as relevant ecological units [START_REF] Ray | Modeling spatial distribution of amphibian populations: A GIS approach based on habitat matrix permeability[END_REF]. We give clear guidelines as to the features of such networks to help maintain amphibian populations in agricultural landscapes. Our results and approach can help stakeholders to find possible compromises between the constraints inherent to agricultural activities and the conservation of biodiversity in these areas.

ACK N OWLED G M ENTS

The Fédération d'Ile-de-France pour la Recherche en Environnement (FIRE FR-3020) is greatly acknowledged for the PhD grant awarded to Alienor Jeliazkov. We also thank Annabelle We are grateful to the landowners, farmers, and authorities from Seine-et-Marne.

TA B L E 1

 1 Risk-based classification of landscape features according to their resistance values (see the main text for the list of references used to assign resistance values to the main land-cover types)

  ) the relative importance of the variables from the different spatial scales; pond, pond-in-network, and network, and in particular, of the different structural, and functional connectivity metrics;(2) the interactions between the variables from the different scales; and (3) the potential thresholds with regard to species distribution. To do so, we modelled the occurrence of the three species in response to all variables of the different scales of influence under the three MMDs (thus, one model per species per MMD = nine models).

  used the boosted regression tree method[START_REF] De'ath | Boosted trees for ecological modeling and prediction[END_REF] (package {dismo} in R Hijmans, Phillips, Leathwick, & Elith, 2016) (Supporting F I G U R E 1 Example of a pond network under the three different connectivity approaches: (a) structural, (b) area-functional, and (c) path-functional, for Bufo bufo and a maximal migration distance (MMD) of 1,000 m. Mean accumulated cost values are weighted distances based on the resistance values of the different types of habitats crossed by the least-cost path and on the pixel resolution (3 × 3 m). The cost unit is metre-equivalent, e.g. to reach one pixel of woodland (of resistance value = 5) through a straight line (length = 3 m), the cost equals 5 × 3 m = 15 metre-equivalentInformation Appendix S5b). We modelled the residual distribution using the logit link function for binary responses. To obtain a reasonable number of trees (>1,000) and good convergence (see[START_REF] Elith | Inventaire des mares d'Ile-de-France. Etat des lieux 2011 et perspectives[END_REF], we used a bag fraction of 0.5 and a learning rate of 0.001 for all models.

  focal pond to the nearest pond occupied by the focal species PondCentral_S Degree centrality in structural network, i.e. the number of edges linking the focal pond to the other ponds of the structural network PondCentral_AF Degree centrality in area-functional network, i.e. the number of edges linking the focal pond to the other ponds of the area-functional network PondCentral_PF Degree centrality in path-functional network, i.e. the number of edges linking the focal pond to the other ponds of the path-functional network Network NetExpan_S The mean inter-pond Euclidean distances within the structural network to which the focal pond belongs NetResist_AF The mean resistance calculated within the area-functional network to which the focal pond belongs NetCost_PF The mean accumulated cost along the least paths connecting the ponds within the path-functional network to which the focal pond belongs MMD. The pond-level relative influence decreased with the increase in MMD in the cases of R. dalmatina and L. vulgaris, while the opposite relationship was observed in the case of B. bufo (Table

F I G U R E 2

 2 Potential migration areas obtained from the resistance map analyses for the three species according to their respective maximal migration distances (MMD): (a) Bufo bufo, (b) Rana dalmatina, and (c) Lissotriton vulgaris. Species occurrence is represented by red (presence) and black (absence) circles F I G U R E 3 Ranking of local and network variables obtained from boosted regression tree analyses for the three species -(a) Bufo bufo, (b) Rana dalmatina, and (c) Lissotriton vulgaris -and for each migration distance considered after model simplification. The variables are coded according to Table 2. DevPerc = the percentage of deviance explained by the final model

F I G U R E 4

 4 Plots of the boosted regression tree-fitted relationships of species occurrence in response to pond-(white), pond-in-network-(light grey), and network-scale (dark grey) variables (Table2) for (a) Bufo bufo in the 1,000 m model, (b) Rana dalmatina in the 1,000 m model, and (c) Lissotriton vulgaris in the 800 m model. Only the variables selected by the final models are shown. The relative percentages of explained deviance are shown in parenthesis of network resistance actually corresponds to a measure of habitat quality rather than connectivity.

  , and might indicate different movement patterns. For instance, natal dispersal has been shown to operate at F I G U R E 5 Plots illustrating the major interaction effects selected and estimated by boosted regression tree models of species occurrence involving (a) only pond, (b, c) pond and pond-in-network, (d, e), pond and network, and (f) network variables for Lissotriton vulgaris under maximal migration distances (MMD) of 100 m (b), 200 m (a), Rana dalmatina under MMD of 300 m (c) and 1,000 m (d, f), and Bufo bufo under MMD of 1,000 m (e)

  have a mean resistance below 40 corresponding to half-woodland and half-crop areas, comprising at least seven ponds (one focal pond connected to six ponds) with a WQI above 70% and separated by <500 m, or by least-cost paths measuring <1,500 m and crossing low resistance features (e.g. woodlands, continuous shrub hedges).Similarly, a potentially viable pond network for R. dalmatina could, as a minimum, have a mean network resistance below 30, include fishless ponds separated by <200 m, and be immediately surrounded by more than 20% of woody areas. Finally, in the case of L. vulgaris, one should favour patchy networks in more open areas with inter-pond distances not exceeding 800 m. These simple statements allow us to highlight the threshold values that could be used as general guidance to orientate pondscape (Boothby, 1997) management in agricultural areas according to a network strategy.

  Boucher (NATUREPARIF), Vincent Vignon (OGE), Pierre Joly (UMR 5023), Dirk Schmeller (UFZ), SNPN, and Département de Seineet-Marne for their help, advice, and support. Thanks are due to Pierre Rivallin (Seine-et-Marne Environnement) for providingearly training with species identification and to Benjamin Mercier (UPMC) for the water analyses. We are also grateful to the student interns, Fabien Michel, Flora Guillier, Victor Court, Pietro Viacava, Nicolas El Battari, Elodie Rey, and others for their help on the field.

Type of cross-scale interaction Number of interactions detected across all models Proportion of interactions detected across all models

  Full list (a) and summary (b) of the most important interactions detected using the boosted regression tree procedure modelling species occurrence in response to pond-(P), pond-in-network-(PN), and network-level (N) variables for Bufo bufo (Buf) under MMD of 1,000, 1,500, and 2,000 m; Rana dalmatina (Ran) under MMD of 300, 500, and 1,000 m; and Lissotriton vulgaris (Tri) under MMD of 100, 200, and 800 m. The interaction size refers to the relative contribution of the interaction to the explained deviance compared to a model where no interaction is considered. In a), Highlighted in grey are the interactions plotted in Figure5

	Species	MMDmin	MMDmed	MMDmax
	Bufo bufo	N ≫ P ≫ PN	P > N > PN	P ≫ N > PN
	Rana dalmatina	P ≫ PN ≫ N	P ≫ N > PN	N > P > PN
	Lissotriton	PN > P > N	PN > P > N	PN > P > N
	vulgaris
	TA B L E 4 (a)
						Type of cross-scale
	Model	Predictor 1	Predictor 2	Interaction size	interaction
	Buf1000	NetExpan_S	WQI	1.9	P × N
		NetResist_AF	WQI	0.61	P × N
	Buf1500	NetCost_PF	DistOccPond	0.38	PN × N
		NetCost_PF	WQI	0.36	P × N
		NetCost_PF	WoodProp	0.13	P × N
	Buf2000	WoodProp	WQI	0.27	P × P
		NetCost_PF	WQI	0.24	P × N
	Ran300	DistOccPond	Fish	1.13	P × PN
		WoodProp	WQI	0.91	P × P
	Ran500	NetExpan_S	Fish	0.27	P × N
		DistOccPond	Fish	0.17	P × PN
		NetResist_AF	DistOccPond	0.12	PN × N
	Ran1000	NetCost_PF	Fish	0.86	P × N
		NetResist_AF	NetExpan_S	0.32	N × N
	Tri100	DistOccPond	WQI	5.41	P × PN
		WQI		Fish	4.92	P × P
	Tri200	DistOccPond	WQI	6.8	P × PN
		WQI		Fish	5.35	P × P
	Tri800	DistOccPond	WQI	5.31	P × PN
		WQI		Fish	3.72	P × N
	(b)
	P × N	8		40
	P × PN	5		25
	P × P	4		20
	PN × N	2		10
	N × N	1		5
	Total	20		100

S U PP O RTI N G I N FO R M ATI O N

Additional supporting information may be found online in the Supporting Information section at the end of the article.


Hanan Hindy

David Brosset
email: david.brosset@ecole-navale.fr

Ethan Bayne
email: e.bayne@abertay.ac.uk

Amar Seeam
email: a.seeam@mdx.ac.mu

Xavier Bellekens
email: x.bellekens@abertay.ac.uk

Improving SIEM for Critical SCADA Water Infrastructures Using Machine Learning

Keywords: Cyber-Physical Systems, Machine Learning, SCADA, SIEM

Network Control Systems (NAC) have been used in many industrial processes. They aim to reduce the human factor burden and efficiently handle the complex process and communication of those systems. Supervisory control and data acquisition (SCADA) systems are used in industrial, infrastructure and facility processes (e.g. manufacturing, fabrication, oil and water pipelines, building ventilation, etc.) Like other Internet of Things (IoT) implementations, SCADA systems are vulnerable to cyber-attacks, therefore, a robust anomaly detection is a major requirement. However, having an accurate anomaly detection system is not an easy task, due to the difficulty to differentiate between cyber-attacks and system internal failures (e.g. hardware failures). In this paper, we present a model that detects anomaly events in a water system controlled by SCADA. Six Machine Learning techniques have been used in building and evaluating the model. The model classifies different anomaly events including hardware failures (e.g. sensor failures), sabotage and cyber-attacks (e.g. DoS and Spoofing). Unlike other detection systems, our proposed work helps in accelerating the mitigation process by notifying the operator with additional information when an anomaly occurs. This additional information includes the probability and confidence level of event(s) occurring. The model is trained and tested using a real-world dataset.

Introduction

An increasing number of networked control systems are being deployed for the monitoring, control and response of physical infrastructure. Deprecated infrastructure is being replaced by SCADA systems that are compatible with IoT Fig. 1: Network Architecture devices, enabling remote monitoring and control of previously isolated systems. Through this upgrade, industries are able to achieve higher reliability and flexibility of deployed systems. However, by enabling external internet access, industries are introducing the increased risk of cyber-attacks to potentially critical infrastructure [START_REF] Ten | Cybersecurity for critical infrastructures: Attack and defense modeling[END_REF].

Vulnerabilities in critical infrastructure protection (CIP) is due to five factors. The first is the lack of open protocols. Current programmable logic controllers (PLCs) are often controlled through proprietary network protocols that have not been scrutinized by security professionals and act as a black box, hence increasing the attack surface [START_REF] Calderón Godoy | Integration of sensor and actuator networks and the scada system to promote the migration of the legacy flexible manufacturing system towards the industry 4.0 concept[END_REF]. The second is due to the lack of segmentation in network infrastructure [START_REF] Jiang | Research of paired industrial firewalls in defensein-depth architecture of integrated manufacturing or production system[END_REF]. Current PLCs are connected to control networks for the data to be analysed in real time by higher management. Whilst this may increase profit margins and decrease the reaction time to market changes, the lack of segmentation allows for transversal network attacks. The third one is due to the number of off-the-shelf equipment being set up in the network [START_REF] Bujari | Standards, security and business models: key challenges for the iot scenario[END_REF]. By including off-the-shelf hardware, the attack surface increases, and introduces other potentially weak links (i.e. Entry points). The fourth one is due to the lack of training received by the operators and the inability of the operator to distinguish between incidents and cyber-attacks [START_REF] Bellekens | Cyber-physical-security model for safety-critical iot infrastructures[END_REF][START_REF] Ahmed | A scada system testbed for cybersecurity and forensic research and pedagogy[END_REF]. Finally, the fifth one is due to the exponential growth of organized cyber-crime and nation state-sponsored cyber-warfare to destabilize countries [START_REF] Jensen | Computer attacks on critical national infrastructure: A use of force invoking the right of self-defense[END_REF][START_REF] Tan | Cyber deterrence in singapore: Framework & recommendations[END_REF].

Critical infrastructure network Security is intrinsically different from computer networks as the interaction between the nodes is done in real time at a physical level. Threats affecting computer networks and their associated intrusion detection systems (IDS) are discussed in [START_REF] Hindy | A taxonomy of malicious traffic for intrusion detection systems[END_REF][START_REF] Hindy | A taxonomy and survey of intrusion detection system design techniques, network threats and datasets[END_REF]. Numerous efforts have been made to apply current cyber-security solutions to CIP networks, however, the solutions proposed are often not suited for the infrastructure and do not consider the underlying risks posed by compromised sensor data [START_REF] Gupta | Editorial security, privacy, and forensics in the critical infrastructure: advances and future directions[END_REF]. As a results, despite the advances made in the design and implementation of cyber-security solutions, there is little research being made with the goal of improving SIEM and increase the resilience of CIP against cyber-incidents. In recent years, the number of attacks against critical infrastructure has significantly increased. The Stuxnet malware is an example of a high-profile SCADA attack, which was first discovered in an Iranian power plant in 2010. Stuxnet was used to infect and reprogramme the PLCs while hiding the changes made using a custom-made PLC rootkit [START_REF] Langner | Stuxnet: Dissecting a cyberwarfare weapon[END_REF]. Another example is the Maroochy attack against a sewage system, causing 800.000L of sewer waste to be released in waterways and parks. The Maroochy attack was perpetrated by an insider, taking advantage of insecure communication between the pumping station and the central SCADA control system [START_REF] Brenner | Eyes wide shut: The growing threat of cyber attacks on industrial control systems[END_REF]. The Ukrainian power grid attack is another example of a large SCADA attack, where 225,000 people were affected by the loss of power. The attacker managed to switch off 30 substations for over 6 hours. In total, over 74 MWh of electricity was not supplied to the energy grid, representing a total of 0.015% of the daily electricity consumption of the country. The attack was later attributed to "Sandworm", a Russian group of hackers known for advanced persistent threats [START_REF] Lee | Analysis of the cyber attack on the ukrainian power grid[END_REF]. These attacks demonstrate that when attacking a CIP, attackers often study the system extensively to perpetrate the attack as stealthily as possible, tailoring their strategies towards the particular system.

In this manuscript, we aim at improving Security Information and Event Management (SIEM) for critical infrastructure using machine learning to identify patterns in the data reported by PLCs in a water system controlled by SCADA. The contributions of this paper are as follows:

-We provide a new SIEM methodology that leverage sensor data and the event-driven nature of cyber-physical systems. -We categorize the attacks in 14 categories, based on the data emerging from the dataset. Out prototype features 3 experiments using machine learning algorithms for the detection of the said anomalies. -We propose the usage of probabilistic model to decrease the mitigation time.

-We conduct a thorough evaluation of SIEM performances through a operational scenarios in the third experiment. The two best predictions are proposed to the operator allowing him to alleviate cyber-attacks when an event is detected. The anomaly detection provides an accuracy of 95.64% with an 85% confidence.

The remainder of this paper is organized as follows; Section 2 discusses equipment used to generate the dataset. Section 3 provides an overview of the dataset, the scenarios and the features, while Section 4 presents the results obtained by applying machine learning technique to differentiate between different scenarios. Section 5 discusses the results and limitations discovered. Section 6 presents related work in this field. The paper ends by concluding our findings in Section 7.

Experimental Set-up

In this section we describe the architecture of SCADA controlled critical infrastructure (CI) used to gather the dataset described in Section 3. Current CI systems are interconnected and are increasingly prone to faulty operations and cyber-attacks. Risks present are caused by the fragmentation of technology incorporated into the CI, vulnerabilities within hardware and software components and through physical tampering of equipment by malicious actors. The design of the CI presented in this research represents a real world, real-time system that is capable of working under normal and abnormal conditions. Additionally, the presented system can be exposed to the aforementioned vulnerabilities.

Figure 1 provides a high level overview of the network of the system. The system is composed of two tanks-the main tank and a secondary tank. Each tank can contain either fuel or water and can be set to two distinct modesacting either as a distributor or as storage. The main and secondary tank has a capacity of 9 and 7 litres respectively. The main tank is composed of four sensors connected to a PLC. The sensors are then connected to both Pump1 and Pump2. Both pumps control the flow of water between the main tank and the secondary tank. The main tank utilizes the four sensors to measure the level of liquid, whilst the secondary tank monitors fluid volume with an ultrasound sensor. Data gathered from the sensors are transferred to the control and monitoring network using the Modbus protocol.

Modbus Protocol

The Modbus protocol operates at the application layer of the open systems interconnection (OSI) model. It enables communication between interconnected network nodes based on a request and reply methodology. The protocol requires little processing overhead, which makes it a sensible candidate for communications between PLCs, Sensors, or Remote Terminal Units (RTUs).

The Modbus protocol works independently of protocols implemented in other layers of the OSI model, hence, it can be used both on routable network or used for serial communications [START_REF] Huitsing | Attack taxonomies for the modbus protocols[END_REF].

The data being stored by the Modbus protocol in slave devices can be categorized in four different ways as listed below:

-Discrete Input Read-Only Access, it provides Physical I/O -Input Register Read-Only Access, it provides Physical I/O -Holding Register Read-Only Access, it provides Read-Write Data -Coil Read-Write Access, it provides application data Each table can contain up to 9,999 values, however, some devices can allow up to 65,536 addresses across all tables. Each vendor has its own specification of the Modbus data tables, and the set-up often requires the operator to read the vendor-specific documentation.

Sensors and PLC

As shown in Figure 1 the system is composed of a single PLC. The PLC is a compact Base Twido PLC, with 40 discrete I/O, 24 discrete inputs, and 14 relays. The system is also composed of a Modicon M238 logic controller. The Fig. 2: System Architecture system can be sent instructions through the touch screen or via a remote system connected to the network.

The system is also composed of five sensors as shown in Figure 1 . Four discrete level sensors are located in the main tank. The first sensor (S0) indicates a low level in the tank (1.25L), the second sensor (S1) indicates a measure of less than 3.35L, the third sensor (S2), indicates a level of 8L while the last sensor (S3) indicates a full tank measure (9L).

The secondary tank is fitted with a single ultrasound sensor. The ultrasound sensor is a Schneider-Electric cylindrical M18 ultrasonic sensor. The sensor is fitted at the top of the tank and is used to measure the distance of the liquid surface to the top of the tank. The sensor can also be used to detect the presence or absence of liquid in the tank.

Operation Mode

The system is controlled automatically through the PLC to avoid the spillage of liquid stored in tanks. Figure 2 illustrates the physical system. For the purpose of this research, the primary tank is filled from a recovery tank, however, the recovery tank in our system can be substituted by other liquid source, such as a river, a fuel line, etc.

Once the primary tank is filled from the recovery tank using pump2, the PLC activates Pump1 to transfer liquid to the secondary tank to avoid spillage.

To simulate a constant liquid consumption from the tanks, valves at the bottom of both the primary and the secondary tanks can be opened. The PLC monitors tank volumes by reading all sensor data registers at 0.1 second intervals. The PLC will automatically refuel the primary tank using Pump 2 when the volume of liquid goes below 1.25L. Similarly, Pump2 will be instructed to stop when a total volume of 9L is reached. The secondary tank will be refuelled when ultrasound sensors detect the liquid level to be below 2.1L. Furthermore, the PLC will only deactivate pump 1 once the ultrasound sensor detects that the secondary tank has reached a level of 6.3L.

Dataset

This section describes the dataset gathered by the IC and the scenarios recorded. The dataset used in this manuscript has been published separately in a dataset publication [START_REF] Laso | Dataset of anomalies and malicious acts in a cyber-physical subsystem[END_REF]. The CI presented in section 2 was used to gather data for the dataset and outputted into CSV format. The dataset comprises actuator and sensor readings that the PLC recorded periodically at 0.1 second intervals. Within the data collected, PLC registers 2 through 4 provided output data describing the state of the system that was used for analysis.

Table 1 provides an overview of the different registers used by the PLC. As shown, Register 2, provides the bits indicating the binary state of the discrete sensors. A population count can be done on the register to retrieve the state of each sensor separately. Register 3 contains the state of the pump, either as active or inactive, while Register 4 contains the step value from 0 to 10,000 of the ultrasound sensors (e.g. Step 3,000 represents 2.1L of liquid in the tank).

Scenarios

The dataset consists of 14 different scenarios as shown in Table 2. Each scenario covers one of 5 operational scenarios representing the potential threat (i.e. sabotage, accident, breakdown, or cyber-attack) as well as 6 affected components. The affected components are system components that are directly affected by the anomaly.

The recorded data is organized in 15 different CSV files of variable duration based on the type of incident recorded (i.e. a sabotage incident may take less time than a distributed denial of service).

Pre-processing

This subsection highlights the pre-processing of the data obtained. The preprocessing is composed of six steps.

Extract Instances

Firstly, each scenario instances are extracted from the log file. An instance is represented by the recording of the register values at a specific time. Each log file has 10 rows per instance. Each row contains The Date, Time, the Register Number, and the Register Value of the PLC.

Calculate rate of change of Register 4

The value of Register 4 is the most significant. However, the significance is not in its single value but in how the value changes over time. This significance can be formulated in the rate of change over time. Figure 3

Rate of change

i = reg4 i -reg4 i-10 time i -time i-10 (1)

Apply Threshold

Table 2 shows the number of recorded instances per scenario. As shown in the table, the events are not statistically distributed over the scenarios. Therefore, the most instances scenario will bias the classification model output. A threshold is applied to take only the first N instances of each file. N should satisfy two conditions a) reduce the gap between instances count per scenario, b) keep the variation of instances count per scenario.

Serialization

A single file is needed for the training process. All Instances are, therefore, serialized.

Normalization 6. Split

Finally, the data is normalized and split into 80% training and 20% testing sets.

Experiments and Results

This section outlines three experiments conducted and how accurately threats could be identified. Note that the code (pre-processing and experiments) is available in a Github repository 1 .

The aim of the experiments is to alert, and provide the operator with the most likely affected components, in order to decrease the time to apply the corrective action.

Six machine learning techniques are used for the classification models, namely, Logistic Regression (LR) Gaussian Naïve Bayes (NB) , k-Nearest Neighbours (k-NN) , Support Vector Machine (SVM) , Decision Trees (DT) and Random Forests (RF) [START_REF] Vanderplas | Python data science handbook: Essential tools for working with data[END_REF]

Parameters

All experiments are conducted using the following computation parameters:

-Training Set : Testing Set = 80% : 20% -Evaluation: The overall accuracy is calculated as follows:

First Experiment

The aim of the first experiment is to alert the user when an anomaly occurs without specifying the associated scenario. The aim is to provide the operator of the CI with a binary output. Figure 4 shows the classification results of the six machine learning algorithms applied to data. The algorithms provide a binary output. As shown, the highest accuracies reached are 93.42%, 93.48% and 94% using DT, RF and k-NN respectively.

Providing an alert that an anomaly exists in the data recorded by the sensors is considered important, however, as the alert is provided in a binary fashion, the operator is -in this case -unable to identify the anomaly at first sight and take corrective action. To this end, a second experiment was established with the aim of proving more information to the operator. Figure 5 shows the classification results of the different machine learning algorithms. The highest accuracies reached are 82.56% and 82.8% using k-NN and RF respectively. The result shows a trade-off between the binary classification offered in the first experiment and by providing the operator with more detailed information. As shown trade-off limits, the accuracy of the system accounts for a large number of false positive, potentially misleading the operator during normal operations.

Third Experiment

The third experiment aims at notifying the operator with an occurring scenario. To this end, the third experiment is divided in three operational trials. Trial One:

In the first trial, models are trained to classify the different scenarios. The operator is notified with the scenario occurring. Figure 6 shows the results of these trial. As shown, the accuracy is reduced and demonstrate a high false negative rate, with the highest accuracy only reaching 81.89% To overcome this high false positive rate, the scenarios were reviewed and the following conclusions were made: (a) The scenarios are co-related. (b) When models are trained to output the probability of an event belonging to another of the 14 scenarios, a maximum of 4 scenarios have non-zero probabilities. For example, Table 3 shows the count of instances having 1, 2, 3 or 4 probable scenarios. In the second row of Table 3, 61.64% of the instances are classified to a single scenario, 31.51% are classified to be one of 2 possible scenarios and only 6.85% are classified to be one of 3 possible scenarios.

As a result of this analysis, trial two to four are based on reporting two probable scenarios to the operator, reducing the uncertainty of our approach.

Trial Two:

In the second trial, the operator is notified with two potential attack scenarios instead of a single one. The two scenarios are the highest probabilities provided by the algorithms. Figure 7 demonstrates that this technique increased the accuracy to reach 95.55% and 95.64% using RF and k-NN receptively. By providing the operator with a probability in the attack scenario he is able to act accordingly and alleviate the attack, hence reducing the overall response time needed. Table 4 shows an example of correctly classified instances when two probable scenarios are considered. The numbers are calculated in reference to the k-NN classifier. For example, in the first row, '2 Floating Objects' scenario misclassified instances are shown . 53 are misclassified as 'Plastic Bag' sabotage. However, 48 of them can be correctly classified by considering the second probable scenario. In the same manner, in Table 4 row 4, 85 instances of the 'Plastic Bag' scenario are misclassified to be '2 Floating Objects' (38 instances), 'Sensor Failure' (24 instances), 'Spoofing' (18 instances), 'Normal' (3 instances) and 'Wrong Con- While the operator is notified with two probable scenarios, it can also be misleading. For this reason our third attempt provides a confidence measure to the operator increasing his situational awareness.

Trial Three and Four:

In the third and fourth trials, a single scenario is reported to the operator unless the probability is less than the confidence percentage. Figure 8 shows the results for 75% confidence. In the first case, a single scenario is reported if its probability is greater than or equal to 0.75, otherwise, if the threshold is reached a second scenario is reported as well. The accuracy provided is 91.84%. This accuracy rises when using an 85% confidence. Figure 8 shows a maximum accuracy of 95.64% using k-NN.

It can be seen from Figure 7 and 8 that DT accuracy is the same. This is due to the DT output, as a single scenario for each instance is available, therefore it is not possible to output two probable scenarios to the operator. Finally, the result shows that the fourth operational trial is the most convenient for the operator as a single scenario is provided, unless the threshold is reached. The accuracy provided by the fourth operational scenario also reaches 95.64%, hence reducing the uncertainty of the attack currently at hand, and the ability of the operator to alleviate the attack. In this section, we discuss the main takeaways based on the experiments evaluated in Section 4. Additionally, the limitations of our models are listed.

Key takeaways:

-Using high confidence: Probabilities are considered to solve the trade-off between the co-related scenarios and reducing the amount of information provided to the operator (i.e. multiple scenarios). This implies reporting the probability of having a scenario occurring. For this piece of work, a confidence interval of 85% is evidenced to provide the best results. -Sensors location and architecture: The position of the sensors and the architecture of the SCADA system affects the collected data and, consequently, affects the model detection accuracy. Moreover, the collection of data greatly affect the pre-processing stages, it is therefore important to gather data consistently. -Longer recording periods needed: In order to increase the accuracy of the machine learning model, the operators should have the recordings of long periods of time for all sensor data. This helps have multiple instances over time and therefore, more instances for training. Moreover, it is recommended to have similar recording duration for the different scenarios. In the case of massively different recording lengths, files need to be serialised in order to avoid bias towards the scenario with the longest recording time. -Importance of scenarios: The more recorded scenarios used for training, the more robust the detection model. With more scenarios, the model can learn variations in scenarios. Moreover, the operator will be provided with higher confidence of reported anomalies.

Limitations:

-Limited Number of Scenarios: The experiments were conducted using 15 log files which covered 14 scenarios. Any new scenario other than the ones covered in the log files will not likely be reported to the operator. -Model Evaluation: Models are evaluated using six machine learning techniques. Hence, the results are limited to these algorithms. -Architecture: The data is limited to the scale of the CI used in this study.

More complex architecture -such as introducing more tanks and sensorsmay not yield the same anomaly detection accuracy. -Real-Time testing: The model is built and tested based on the dataset, however, testing the model real-time is essential to measure the performance.

Related Work

The area of critical infrastructure protection and Security Information and Event Management has a rich history of books, deployments and lessons learned from both universities and organizations deploying tools and techniques to protect their environments. Robert Mitchell and Ing-Ray Chen [START_REF] Mitchell | A survey of intrusion detection techniques for cyber-physical systems[END_REF] presented a survey on recent Cyber-Physical Systems (CPS) intrusion detection systems. They classified CPSs into two categories based on potential detection technique: (1) Knowledge-Based or Misuse Detection and (2) Behaviour-Based. While the survey provides a general overview of potential threats, it focuses essentially on the network traffic rather than on data provenance and data accuracy. This survey, however, provides the stepping stone of the field. In [START_REF] Amin | Cyber security of water scada systemspart i: Analysis and experimentation of stealthy deception attacks[END_REF] and [START_REF] Amin | Cyber security of water scada systemspart ii: Attack detection using enhanced hydrodynamic models[END_REF] Amin et al. provide an overview a security threat assessment of a networked CI including different layers (i.e. supervisory and control networks). The manuscripts present a grey-box approach where the hacker possesses a certain level of knowledge of the system and is able to perpetrate a deception attack against the system in order to enable liquid spilling from the canal the IDS was tested on. Cheng et al. [START_REF] Cheng | Orpheus: Enforcing cyber-physical execution semantics to defend against data-oriented attacks[END_REF] provide a technique to alleviate control-oriented attacks, code-injection attacks or code-reuse attacks on embedded devices. The highlight the lack of existing mechanisms and present 'Orpheus', a program behaviour model, taking advantage of the eventdriven nature of embedded devices controlling critical infrastructure. Mathur [START_REF] Mathur | On the limits of detecting process anomalies in critical infrastructure[END_REF] discusses the limitations of the detecting a incidents in critical infrastructures by analysing processes. The manuscript refers to two methods, the first described is the CUSUM [START_REF] Cárdenas | Attacks against process control systems: risk assessment, detection, and response[END_REF]. CUMSUM is a statistical method allowing to detect anomalies in time series, corresponding to a specified process. The technique requires two parameters to operate, the 'bias' and thethreshold'. The CUMSUM provides the cumulative sum of the deviation for the process measured. By plotting the predicted and observed state, the operator is able to identify the state of the facility and identify changes in behavior. The second methodology is based on State Entanglement (SE) [START_REF] Adepu | Distributed detection of single-stage multipoint cyber attacks in a water treatment plant[END_REF]. SE combines the states of multiple components of a system to construct a state space. The state space act as a blacklist, highlighting prohibited states during normal operations.

Conclusion and Future Work

This work focused on building an anomaly detection and a SIEM tool for a SCADA water system. The model was evaluated using a real-world dataset covering 14 anomaly scenarios including normal system behavior. The presented scenarios covered a wide range of events, ranging from hardware failure to sabotage. Three experiments were conducted using 6 Machine Learning techniques.

The experiments varied based on the level of information reported to the operator. The First experiment allowed anomaly event to be reported to the operator as a binary output. While events where being detected as either anomalies or normal operation, the operator was unaware of the type of anomaly occurring. The second experiment reported the affected component, providing the operator with information relating to a single or multiple sensor data. Finally, the third experiment -which provided the best results -reported the anomaly with an accuracy level, helping the operator to take subsequent correcting steps.

The overall evaluation showed that k-NN, Decision Tree and Random Forest outperformed 'Gaussian Naïve Bayes, SVM and kernel SVM. Moreover, k-NN results demonstrated the highest accuracy amongst all algorithms in the three experiments. The experiments achieved the accuracy of 94% for the binary output and 95.64% for the scenarios classification. As aforementioned, the scenarios are co-related, therefore, the operator -in the third experiment -is notified with the most probable scenarios/anomaly occurring. Moreover, a confidence level was used to provide the operator with the best information available.

To further enhance the detection accuracy, the following should be considered. Increasing the dataset instances to enhance the training process and building hybrid model to classify subgroups of events.

Fig. 3 :

 3 Fig. 3: Rate of change of register 4 value over different scenarios

  [START_REF] Jr | Applied logistic regression[END_REF] [23][START_REF] Steinwart | Support vector machines[END_REF] [20][START_REF] Barber | Bayesian reasoning and machine learning[END_REF].

  OverallAccuracy =T P + T N T P + T N + F P + F N Where:• True Positive (TP): Number of anomaly instances which scenario is correctly detected • True Negative (TN): Number of normal instances which are correctly detected • False Positive (FP): Number of normal instances which are detected as one of the anomaly scenarios • False Negative (FN): Number of anomaly instances which are detected as normal scenarios

Fig. 4 :Fig. 6 :Fig. 7 :

 467 Fig. 4: Anomaly Classification Accuracy Fig. 5: Affected Component Classification Accuracy

Fig. 8 :

 8 Fig. 8: Scenarios Classification Accuracy (Alerting with one/two scenario(s) based on 75% and 85% confidence)

Table 1 :

 1 Registers Extracted Bits Representation

	Reg. No. 2 4	Bit No. 4 5 6 7 16-bits	Value Discrete Sensor 3 Discrete Sensor 2 Discrete Sensor 1 Discrete Sensor 0 Depth Sensor	Reg. No. 3	Bit No. 0 1 5 4	Value Pump 2 Pump 1 Pump 2 Valve Pump 1 Valve

Table 2 :

 2 Dataset Scenarios, Operational Scenarios and Affected Components

		Scenario	Affected Component	Operational Scenario	No. of instances
	1	Normal	None	Normal	5519
	2	Plastic bag		Accident/Sabotage	10549
	3 4	Blocked measure 1 Blocked measure 2 Floating objects in	Ultrasound Sensor	Breakdown/Sabotage	226 144
	5	main tank (2 objects)		Accident/Sabotage	854
		Floating objects in
	6	main tank (7			733
		objects)
	7	Humidity			157
	8	Discrete sensor failure	Discrete sensor 1	Breakdown	1920
	9	Discrete sensor failure	Discrete sensor 2		5701
	10 11	Denial of service attack Spoofing	Network	Cyber-attack	307 10130
	12	Wrong connection		Breakdown/Sabotage	6228
	13	Person hitting the tanks (low intensity) Whole subsystem Person hitting the	Sabotage	347
	14	tanks (medium			281
		intensity)
	15	Person hitting the tanks (high intensity)			292

  shows the rate change in Register 4 value over 15 different scenarios.

		10000
	Register 4 Value	2000 4000 6000 8000	Person htting low Person htting med Person htting high Blocked measure 1 Blocked measure 2 2 Floating objects 7 Floating objects Discrete Sensor 1 Discrete Sensor 2 Wrong connection DoS Spoofing Plastic bag Humidity
		0
		0	50 100 150 200 250 300 350 400 Time (in sec)

Table 3 :

 3 Distribution of probabilistic classification of scenarios

	Algorithm	Maximum number of probable scenarios per instance	Number of instance with 1 probable scenario	Number of instance with 2 probable scenarios	Number of instance with 3 probable scenarios	Number of instance with 4 probable scenarios
	DT	1	3053 100%	--	--	--
	k-NN	3	1882 61.64%	962 31.51%	209 6.85%	--
	RF	4	1737 56.9%	1004 32.89%	297 9.71%	15 0.5%

Table 4 :

 4 Co-relation of scenarios that are misclassified based on highest probable scenario and correct with the second probable one (Calculated based on k-NN experiment)

		Instances count		Scenario (Y) The count of instances classified
		where X			as Y while the correct is X
	Scenario (X)	Is Not 1 st Scenario	Is 2 nd Scenario	2 Floating Objects	7 Floating Objects	Normal	Plastic Bag	Sensor Failure	Spoofing	Wrong Connec-tion
	2
	Floating	53	48	-	-	-	53	-	-	-
	Objects
	7
	Floating	5	5	-	-	-	-	-	-	5
	Objects
	Normal	80	58	-	-	-	1	68	-	11
	Plastic Bag	85	75	38	-	3	-	24	18	2
	Sensor Failure	183	137	-	-	75	53	-	5	50
	Spoofing	45	25	-	-	-	37	8	-	-
	Wrong
	Connec-	102	72	-	5	23	5	69	-	-
	tion
	nection'(2 instances). 75 can be correctly classified with the second probable
	scenario in consideration.

https://github.com/AbertayMachineLearningGroup/machine-learning-SIEMwater-infrastructure