Mapping Ecological Infrastructure in a Cross-Border Regional Context

Abstract

Facing the decline of biodiversity worldwide, the conservation of the remaining natural and semi-natural areas is fundamental. To do so, the concept of green infrastructure has gained attention recently. This case study presents the method developed to identify the green infrastructure in a cross-border, urbanized territory between Switzerland and France in the area of influence of the city of Geneva. The first part of the methodology consists of calculating and mapping the inputs aggregated in four pillars: (i) the distribution of habitats as well as the predicted distribution of hundreds of plant and animal species, (ii) the supply of five ecosystem services, (iii) the functional connectivity for three animal species and the light pollution and (iv) five indices of landscape’s structure. These inputs are then used to run a prioritization model to identify the areas with the highest ecological interest according to these weighted inputs. The cross-border situation of this case study had impacts on the way the input data were gathered and weighted and on the way the output was created to consider the expectations of the three main local authorities involved, without creating any legal obligations on the implementation of the green infrastructure. As a positive sign of the usefulness of these results, the resulting maps were immediately transferred to the land use planners in charge of developing ambitious visions of the “Grand Genève” territory for 2050 in alignment with 10 objectives of ecological transition as recently agreed and signed by local authorities. The method presented in this article is flexible and includes a broad description of biodiversity, supporting a reliable network of areas with high ecological values for conservation purposes and human well-being.

1. Introduction

Global biodiversity is under a major crisis at every level and its genetic, species and ecosystem diversity is declining rapidly [1,2]. The destruction of natural habitats for agriculture and/or urbanization is the main cause of its decline and is directly linked to our way of occupying terrestrial and marine surfaces as well as to our consumption patterns [3,4,5,6,7]. Biological diversity decline might ultimately alter ecosystem functions such as productivity, stability and resilience, jeopardizing our food and water security as well as our socio-economic well-being [8,9,10,11,12,13,14,15,16,17]. Thus, the conservation of the remaining natural and semi-natural areas is fundamental, especially in urbanized environments where urbanization represents an additional pressure.

Green infrastructure (GI) is defined as a network of (semi-)natural areas allowing structural and functional connectivity of the landscape where biodiversity and ecosystem services are concentrated. The concept of GI fits perfectly the modern view of nature’s conservation that emphasizes the cohabitation of people and nature with sustainable and resilient interactions [18]. This new paradigm encompasses the common health of human societies and natural systems, highlighting our direct dependence upon ecosystems as described in the “One Health” concept [19,20]. GI is usually described as an interconnected network of (semi-)natural areas designed to deliver wide range of ecological, social and economic benefits [21,22,23], although several definitions have been used [24,25]. It is usually made up of large areas concentrating biological diversity and ecosystem service supply, linked with corridors allowing structural and functional connectivity [22,26]. GI is highly relevant in urban areas because it gives an ecological value to each element of the territory, focusing on the multifunctionality of the landscape. GI also integrates the nature-based solutions to mitigate the effects of global changes [26,27,28]. It is promoted at the European scale but also at the federal and cantonal scale in Switzerland [29,30,31]. In France, GI contributes directly to the policy on green and blue networks (http://www.trameverteetbleue.fr (accessed on 1 July 2023). Internationally, GI fits Target 3 of the urgent actions that need to be taken over the decade to 2030 from the 15th Conference of the Parties to the Convention on Biological Diversity, which proclaims to “Ensure and enable that by 2030 at least 30 per cent of terrestrial areas of particular importance for biodiversity and ecosystem functions and services, are effectively conserved and managed through ecologically representative, well-connected and equitably governed systems of protected areas and other effective area-based conservation measures, while ensuring that any sustainable use is fully consistent with conservation outcomes.” [32].

There is no consensus on the methodology nor the inputs that should be used in order to identify and map a GI [24]. This has led to confusions where the term “green infrastructure” was used in very divergent ways while several concepts and terms were referring to the same idea (e.g., ecological network, green corridors, green prints, etc.) [25]. For example, in highly urbanized environments, GI is often used as a greening method or as architectural elements such as green walls or green roofs [33]. In other situations, GI is restricted to areas supporting ecosystem supply only or protected lands [34,35,36,37,38,39]. More details of how GI is used in the scientific literature as well as the methods employed to identify it and the associated limits are available in Honeck et al. (2020) [26]. This literature review identified a methodological gap in the identification of GI where most of the articles do not consider all aspects of biodiversity conservation and GI’s definition.

The methodology employed here is based on the “three pillars” approach that has already been applied in Geneva, Switzerland [40]. This approach allows the consideration of all aspects of biodiversity conservation and respects the initial definition of GI [26]. The method is adapted here to aggregate the inputs in four pillars, the third and the fourth being initially grouped: (1) the diversity and distribution of species and habitats using species distribution models and a land use–land cover (LULC) map of the territory, (2) the supply of ecosystem services and (3) the functional and (4) the structural connectivity of the landscape. A spatial prioritization tool is used to select the network of areas with the highest ecological interest. The novelty of this article is the presentation of an application of the theoretical approach developed in Honeck et al. (2020) [26] on a cross-border territory, emphasizing the various methods used to calculate and integrate 2437 inputs for the identification of the GI. This exhaustive work can be used as a baseline for any territory aiming at mapping its own GI.

This article is focused on the prioritization of the GI network covering 30% of the territory of a regional cross-border agglomeration between France and Switzerland named “Greater Geneva”. This region is located in the European Alps and in the area of economic influence of the city of Geneva. In January 2023, the elected representatives of Greater Geneva signed the Charter for Greater Geneva in Transition with the desire to make the ecological transition the backbone of cross-border cooperation, recognizing that the erosion of life, the depletion of natural resources and climate degradation are our greatest threats [41]. The Charter sets out 10 strategic commitments to respect both the social floor and the ecological ceiling. The GI is fully in line with Objective 3 on Biodiversity of the Charter that aims at stopping the loss of natural habitats by 2050. It is also expected to have a positive impact on all other objectives.

This particular transboundary setting generated several issues that were addressed by the research questions presented below.

(i) What are the difficulties of gathering input data across borders?; (ii) What is the distribution of prioritization value across the studied areas?; (iii) What is the best 30 percent of the territory?; (iv) How can we accommodate the desire of each administrative entity to identify its own best 30%?; (v) What share of the identified GI is already protected? (vi) What are the difficulties in establishing GIs across borders?

2. Materials and Methods

2.1. Study Area

Greater Geneva is a cross-border territory between Switzerland and France of approximately 2000 km² located around the city of Geneva. In its strict limits, it integrates three administrative entities grouped in two Swiss cantons (“Genève” and “Vaud” with the District of Nyon), and two French Departments (“Ain” and “Haute-Savoie”) with the “Pôle Métropolitain du Genevois Français”. This peculiar territory induces difficulties in compiling data because the taxonomy, the methods and the data availability vary from one administrative entity to another. However, the territory has a biogeographic consistency and is delimitated by mountain ranges, the Alps in the south-west and the Jura in the north-east, which justifies the GI assessment at this scale beyond borders. The region is particularly dynamic, and the population is growing rapidly due to the attractivity of the city of Geneva. The territory is dominated by urbanized areas and crops in the lowlands and forests and pastures in the mountainous areas.

2.2. Method for Mapping the Ecological Infrastructure

According to the definition of GI, the inputs used to identify and map it have to consider several aspects of biodiversity and ecosystem service conservation [22,26,42]. For clarity, they were grouped into four main pillars: (1) the diversity pillar that includes the assessment of species and habitat distributions based on models and the aggregation of available LULC data; (2) the ecosystem service pillar that aims at mapping the supply of five regulating ecosystem services; (3) the connectivity pillar that includes maps of the functional connection for three animal species and light pollution; and (4) structural indices of the landscape based on the LULC categories. Once the inputs have been prepared at a spatial resolution of 25 m, they are included in a spatial prioritization tool set to classify every pixel of the territory according to its relative importance for biodiversity and ecosystems service conservation [42,43] (Figure 1). The method and the theoretical background used here were developed and explained exhaustively in two papers [26,40]. The inputs were selected according to the available data for describing the four pillars, their collinearity and their ecological meaning. The weight attributed to each of them was discussed in the working team and with the main stakeholders. The following sections explain the methods used to calculate and map the inputs.

2.3. Pillar 1: Species and Habitat Distribution

The assessment of the distributions of many species of plants and animals allows the identification of important areas for the conservation of specific richness. Furthermore, the inclusion of all distributions into the spatial prioritization process allows the selection of areas that are of the highest importance for rare species. The distribution of natural habitats also plays an important role in nature’s conservation by providing food and shelter to animals and plants but also by the maintenance of their ecological functions.

2.3.1. Natural Habitats

The distribution of habitats and more globally the LULC information are highly important for spatial planning at the regional scale but also for species distribution models (SDMs) [44].

The LULC map was created based on the compilation of the French and Swiss geomatic information sources, respectively named “Institut national de l’information géographique et forestière” (IGN, https://www.ign.fr/institut) and “Topographic Land Model” from SwissTopo (TLM3D, https://www.swisstopo.admin.ch/fr/geodata/landscape/tlm3d.html) (accessed on 1 July 2023). The data available across the study area were heterogenous in typology and geometry so both sources of information were used to homogenize the various maps into one. The geometry of the IGN map was extracted to divide the territory by administrative parcels that were transformed into polygons. Then, the habitat maps were transformed into five-meters points and added to the polygons where the most represented habitat was selected for each parcel. The dense urban environments, as well as the roads, railways, highways, rivers and running water, were then added. Finally, the diffuse urban environment class was created based on the presence of vegetation in the urban classes using NDVI information. More details of the method are presented in Figure 2.

The categories that were used in the prioritization process represent 19 (semi-)natural habitats, with the urban ones being excluded from this analysis. Each selected category was extracted as a unique input, ensuring the selection of at least a part of each (semi-) natural habitat in the final GI network by the prioritization process.

2.3.2. Species Distribution Modeling

SDM allows the creation of a covering map of habitat suitability based on the georeferenced observations of species’ individuals and predictive variables [45,46,47]. Several methods exist and have been used extensively in conservation [48,49,50,51]. Species’ occurrences were compiled from French and Swiss botanical conservatories and monitoring programs. Only observations between 2000 and 2020 and with a precision below 25 m were kept. To better conserve endangered species, we compiled the red list statuses from the different entities and selected the most threatened status. This ensures that the threats species are facing are not under-evaluated. At the end of the process, 585 species of animals and 1816 plants were selected. Predictive variables were selected based on their collinearity, ecological meaning and modeling performances in the study area and are presented in Table 1 [52]. The resolution of these variables is 25 m.

The chosen modeling algorithm was MaxEnt [53,54] (version 3.4.1) because it is widely used in SDM [55] and known to perform well especially with presence-only data [56,57]. The models were run using “Dismo” [58], “ENMeval” [59] and “sdm” [60] packages in R [61]. The default settings were kept except for the beta multiplier that was set to 2.00 to avoid over-fitting [62,63]. For each model, the occurrences were randomly split with 75% used for calibration and 25% for evaluating the model’s performances 10 times in a raw. Ten thousand background data were randomly created for each model. Then, for each species, a final model was calibrated with all occurrences available to map habitat suitability with all the information. More details about the modeling method and data selection process can be found in Sanguet et al. (2022) [52].

Table 1. Predictive variables used in the SDM.

Variables	Description	Origin
Temperatures	Mean annual temperature	Worldclim, R
Precipitations	Annual precipitations	Worldclim, R
Exposition	Northness index	ArcMap 10.2.1
Slope	Continuous slope	ArcMap 10.2.1
Solar radiations	Mean seasonal solar radiation	ArcMap 10.2.1
Landscape dominance	Index of landscape domination	ArcMap 10.2.1, modified from Weiss, 2001 [64]
Cambisol	Cambisol proportion in the surrounding soil	Hengl et al., 2017 [65]
Podzol	Podzol proportion in the surrounding soil	Hengl et al., 2017 [65]
Closed forests	Distribution of deciduous and coniferous forests	LULC map
Open forests	Distribution of opens forests and barrens	LULC map
Urban areas	Distribution of highly and moderately dense urban areas	LULC map
Transportation	Distribution of railways, paths, highways and roads	LULC map
Disturbed vegetation	Distribution of urban and wooded disturbed vegetation	LULC map
Natural meadows	Distribution of dry, alpine and extensive meadows	LULC map
Agriculture	Distribution of crops, vineyards and orchards	LULC map
Wetlands	Distribution of wet meadows, riverbeds and wet forests	LULC map

Table 1. Predictive variables used in the SDM.

Variables	Description	Origin
Temperatures	Mean annual temperature	Worldclim, R
Precipitations	Annual precipitations	Worldclim, R
Exposition	Northness index	ArcMap 10.2.1
Slope	Continuous slope	ArcMap 10.2.1
Solar radiations	Mean seasonal solar radiation	ArcMap 10.2.1
Landscape dominance	Index of landscape domination	ArcMap 10.2.1, modified from Weiss, 2001 [64]
Cambisol	Cambisol proportion in the surrounding soil	Hengl et al., 2017 [65]
Podzol	Podzol proportion in the surrounding soil	Hengl et al., 2017 [65]
Closed forests	Distribution of deciduous and coniferous forests	LULC map
Open forests	Distribution of opens forests and barrens	LULC map
Urban areas	Distribution of highly and moderately dense urban areas	LULC map
Transportation	Distribution of railways, paths, highways and roads	LULC map
Disturbed vegetation	Distribution of urban and wooded disturbed vegetation	LULC map
Natural meadows	Distribution of dry, alpine and extensive meadows	LULC map
Agriculture	Distribution of crops, vineyards and orchards	LULC map
Wetlands	Distribution of wet meadows, riverbeds and wet forests	LULC map

2.4. Pillar 2: Ecosystem Service Supply

The preservation of ecosystem services of regulation and support to biodiversity contributes to conserving the good functioning of ecosystems. Furthermore, the preservation of ecosystem services’ supply, as part of nature’s contributions to people and nature-based solutions, helps mitigate the detrimental effects of climate change or the consequences of extreme meteorological events. Other types of ecosystem services such as resource production or cultural services were not included in this work because they might induce the selection of areas with low quality or detrimental effects on biodiversity. Five ecosystem services were modeled and mapped using InVEST (version 3.12.0): the suitable areas for pollinators, the atmospheric carbon storage, the nutrient delivery ratio, the sediment delivery ratio and the leaf area index.

2.4.1. Suitable Areas for Pollinators

Pollinators are highly important for crop pollination and, as a consequence, for our food provision. Identifying their most suitable habitats to be integrated into the final GI network participates in maintaining the populations’ health and pollinators’ availability. The index calculated here represents a potential abundance of pollinators for each pixel of the resulting map based on their ecology, considering the quality and attractivity of the habitats for feeding and nesting habits. The “Crop Pollination” program was used in the software InVEST. Two tables are needed to run the model. The guild table contains the characteristics of 20 wild bee species while the biophysical table associates habitats with wild species’ traits and habits. Both tables were created based on a literature review and local expert knowledge to adjust the value to the local context (

Table A1. Biophysical table used for the InVEST pollination model.

Code	Category	Nesting Cavity Availability Index	Nesting Ground Availability Index	Floral Resources Spring Index	Floral Resources Summer Index	Floral Resources Autumn Index
1	Highways	0	0	0	0	0
2	Other wetlands	0	0	0.3	0.5	0.2
3	Path	0	0	0	0	0
4	Railways	0	0	0	0.7	0.3
5	Crops	0	0	1	0	0
6	Still water	0	0	0	0	0
7	Running water	0	0	0	0	0
8	Coniferous forests	0	0	0	0	0
9	Deciduous forests	0	0	0	0	0
10	Wet forests	1	0	0.7	0.3	0
11	Open forests	0.5	0.5	0.4	0.5	0.1
12	Riverbeds	0	1	0.3	0.6	0.1
13	Gravel pits	0	1	0.2	0.6	0.2
14	Barrens	0	1	0.1	0.8	0.1
15	Rocks and cliffs	0	1	0	0	0
16	Alpine grasslands	0	1	0.1	0.8	0.1
17	Extensive grasslands	0	1	0.6	0.2	0.2
18	Wet meadows	1	0	0.3	0.5	0.2
19	Dry meadows	0.3	0.7	0.6	0.2	0.2
20	Rivers	0	0	0	0	0
21	Roads	0	0	0	0	0
22	Stream	0	0	0	0	0
23	Dense urban	0	0	0	0	0
24	Diffuse urban	0.5	0.5	0.6	0.3	0.1
25	Wooded disturbed vegetation	0.5	0.5	0.8	0.2	0
26	Urban vegetation	0	0	0.7	0.3	0
27	Orchards	0	1	0.9	0.1	0
28	Vineyards	0	1	0.9	0.1	0

and

Table A2. Guild table used for the InVEST pollination model.

Species	Nesting Suitability Cavity Index	Nesting Suitability Ground Index	Foraging Activity Spring Index	Foraging Activity Summer Index	Foraging Activity Autumn Index	Alpha	Relative Abundance
Andrena carantonica	0	1	0.7	0.3	0	512	1
Andrena chrysosceles	0	1	0.6	0.4	0	260	1
Andrena cineraria	0	1	0.6	0.4	0	300	1
Andrena dorsata	0	1	0.4	0.6	0	650	1
Andrena flavipes	0	1	0.3	0.7	0	1150	1
Andrena fulva	0	1	1	0	0	315	1
Andrena haemorrhoa	0	1	0.8	0.2	0	373	1
Andrena minutula	0	1	0.6	0.4	0	112	1
Andrena nitida	0	1	0.8	0.2	0	288	1
Bombus hortorum	0.5	0.5	0.3	0.6	0.1	604	1
Bombus hypnorum	1	0	0.3	0.6	0.1	288	1
Bombus lapidarius	0.5	0.5	0.2	0.6	0.2	1500	1
Bombus pascuorum	0.5	0.5	0.3	0.4	0.3	2300	1
Bombus pratorum	0.5	0.5	0.4	0.5	0.1	674	1
Bombus terrestris	0.5	0.5	0.3	0.5	0.2	1500	1
Lasioglossum calceatum	0	1	0.4	0.4	0.2	1000	1
Lasioglossum malachurum	0	1	0.3	0.5	0.3	600	1
Lasioglossum morio	0	1	0.2	0.4	0.4	69	1
Lasioglossum politum	0	1	0.2	0.6	0.2	14	1
Osmia rufa	1	0	0.6	0.4	0	600	1

in Appendix A). The optional farm map was not used. The model produces one map for each season (winter excluded) which were added to map the total suitability of the landscape.

2.4.2. Atmospheric Carbon Storage

Carbon dioxide (CO₂) is a greenhouse gas massively rejected into the atmosphere by human activities and is the main cause of the observed global warming. The preservation of natural habitats known to store carbon avoids their destruction and thus the release of the stored carbon dioxide into the atmosphere. Preserving forest growth compensates for a part of our emissions and participates sequestrating carbon into organic matter and the soil. The mapping of this ecosystem service uses a biophysical table linking each habitat category to its carbon storage capacity. The table, named “Carbon pools” in the “Carbon Storage and Sequestration” program of InVEST, was adapted from the available data in InVEST’s documentation and is available in Appendix A (

Table A3. Biophysical table used for the InVEST carbon storage model.

Code	Category	C_above	C_dead	C_below	C_soil
1	Highways	0	0	0	0
2	Other wetlands	6.5	0	0	68.23
3	Path	0	0	0	0
4	Railways	0	0	0	0
5	Crops	3.51	0	0	49.91
6	Still water	0	0	0	0
7	Running water	0	0	0	0
8	Coniferous forests	134.51	8.10	8.70	55.4
9	Deciduous forests	134.51	8.10	8.70	55.4
10	Wet forests	134.51	8.10	8.70	55.4
11	Open forests	134.51	8.10	8.70	55.4
12	Riverbeds	7.16	0	0	26.31
13	Gravel pits	0	0	0	0
14	Barrens	6.5	0	0	68.23
15	Rocks and cliffs	0	0	0	0
16	Alpine grasslands	4.98	0	0	59.4
17	Extensive grasslands	4.98	0	0	59.4
18	Wet meadows	6.4	0	0	68.23
19	Dry meadows	4.98	0	0	59.4
20	Rivers	0	0	0	0
21	Roads	0	0	0	0
22	Stream	0	0	0	0
23	Dense urban	0	0	0	0
24	Diffuse urban	0	0	0	0
25	Wooded disturbed vegetation	20.45	0	8.7	55.4
26	Urban vegetation	15.43	0	0	53.4
27	Orchards	22.73	0	0	64.76
28	Vineyards	5.36	0	0	53.4

). Only the storage was measured and not the sequestration.

2.4.3. Nutrient Delivery Ratio

The “Nutrient Delivery Ratio” (NDR) program of InVEST calculates the flow of nutrients into the rivers and other water bodies or their retention in the soil’s upper layers. Excessive nutrient accumulation in water could impacts aquatic ecosystems composition and functioning. The NDR program models the landscape’s load and retention of nitrogen and phosphorous based on a biophysical table linking the LULC categories to their nutrient retention and load abilities, but also on the digital elevation model, a nutrient runoff proxy such as the annual precipitations, and the distribution of watersheds. The values of the biophysical table were adapted to the local context from the literature [66,67,68] as well as from the available data on the InVEST documentation [69] and is available in the Appendix A (

Table A4. Biophysical tables used for the InVEST sediment and nutrient delivery models (SDR and NDR).

Code	Category	SDR		NDR
		usle_c	usle_p	load_n	eff_n	crit_len_n	proportion_subsurface_n	load_p	eff_p	crit_len_p	usle_c
1	Highways	0	0.5	16.6	0.1	25	0	1.18	0.01	25	0
2	Other wetlands	0	1	16.7	0.5	150	0	0.01	0.8	150	0
3	Path	0	0.5	16.6	0.1	25	0	1.18	0.01	25	0
4	Railways	0	0.5	16.8	0.1	25	0	1.18	0.01	25	0
5	Crops	0.1384	0.9942	155	0.65	150	0.5	44.6	0.4	150	0.1384
6	Still water	0	1	11.4	0.05	25	0	0.01	0.01	25	0
7	Running water	0	1	17.1	0.05	25	0	0.01	0.01	25	0
8	Coniferous forests	0.0012	1	16.8	0.9	150	0	2.36	0.6	150	0.0012
9	Deciduous forests	0.0012	1	16.8	0.9	150	0	2.36	0.6	150	0.0012
10	Wet forests	0.0012	1	16.7	0.5	150	0	0.01	0.8	150	0.0012
11	Open forests	0.0012	1	14.3	0.9	150	0	2.36	0.6	150	0.0012
12	Riverbeds	0	1	7.3	0.1	25	0	1.18	0.01	25	0
13	Gravel pits	0	1	7.3	0.1	25	0	1.18	0.01	25	0
14	Barrens	0.0219	1	10.3	0.9	150	0	2.36	0.6	150	0.0219
15	Rocks and cliffs	0	1	7.3	0.1	25	0	1.18	0.01	25	0
16	Alpine grasslands	0.0903	0.9942	12.6	0.65	150	0.5	10.5	0.5	150	0.0903
17	Extensive grasslands	0.0903	0.9942	74	0.65	150	0.5	27.55	0.5	150	0.0903
18	Wet meadows	0.0435	1	16.7	0.5	150	0	0.01	0.8	150	0.0435
19	Dry meadows	0.0435	1	9.2	0.1	150	0.5	2.36	0.6	150	0.0435
20	Rivers	0	1	17.1	0.05	25	0	0.01	0.01	25	0
21	Roads	0	0.5	16.6	0.1	25	0	1.18	0.01	25	0
22	Stream	0	1	17.1	0.05	25	0	0.01	0.01	25	0
23	Dense urban	0	0.5	17	0.1	25	0	19.1	0.01	25	0
24	Diffuse urban	0	0.5	16.4	0.4	25	0	19.1	0.01	25	0
25	Wooded disturbed vegetation	0.0219	1	14.7	0.4	150	0	2.36	0.6	150	0.0219
26	Urban vegetation	0.02652	0.5	15.7	0.4	150	0.4	2.36	0.6	150	0.02652
27	Orchards	0.1232	0.9942	74	0.65	150	0.5	44.6	0.4	150	0.1232
28	Vineyards	0.3527	0.9942	74	0.65	150	0.5	44.6	0.4	150	0.3527

). The Borselli K parameter was set to 2, the subsurface critical length to 200 m, the subsurface maximum retention efficiency to 0.8 and after several tests the threshold flow accumulation was set to 140. This last parameter adjusts the modelling of the location of temporary rivers based on the digital elevation model. The value was selected after several tests to better fit known permanent and temporary rivers on the territory. Several results were produced, and the effective retention map was kept. It represents the relative capacity of each pixel to retain nutrients.

2.4.4. Sediment Delivery Ratio

The “Sediment Delivery Ratio” program in InvEST models the flow of sediments and thus the erosion of the landscape. Erosion might induce a higher risk of landslides and a loss of organic matter in the soil. The preservation of areas reducing the risk of erosion is a nature-based solution and allows the mitigation and avoidance of natural hazards. The data and settings required for this model are relatively similar to the NDR model. The biophysical table was adapted from the existing data in the InVEST documentation [69] and the literature [70,71] and is available in Appendix A (Table A4). The values link each LULC category with its ability to reduce the loss in sediment and its management by humans to reduce the erosion. The erosivity and erodibility maps were downloaded from the European Soil Data Centre [72] and projected in the territory at 25 m resolution. The settings used were the following: threshold flow accumulation = 140, Borselli K parameter = 2, maximum SDTR value = 0.8, Borselli ICO parameter = 0.5, maximum L value = 122. The results are composed of several maps and the avoided sediment export was selected. It gives a value to each pixel according to their ability to avoid sediment export.

2.4.5. Leaf Area Index

Vegetation cover reduces the temperature and filters the air. Due to mitigating the effects of climate change and regulating the micro-climate, it is especially interesting to preserve green spaces in urban environments. This ability could be mapped by the leaf area index based on remote sensing images of the territory. The normalized differentiation vegetation index (NDVI) allows mapping of the greenness of a landscape and has been largely used to classify vegetation types [73,74]. The average maximum value of the NDVI in the territory was calculated based on the remote sensing images from Landsat-5, Landsat-7 and Landsat-8 and compiled in the Swiss Datacube [52,73,74,75]. Then, a formula was applicated to the raster to calculate the leaf area index (1) [76].

$$0.57\ast e^{(2.33\ast NDVI)}$$

(1)

2.5. Pillar 3: Functional Connectivity

Functional connectivity ensures spatial connections between habitats and maintains species movements which are especially important for their resilience against climate change [77] and for gene-flow. According to their characteristics, shape, surface or fragmentation level, natural habitats’ quality and functions vary [78]. Functional connectivity corresponds to the relative ease of mobility in the landscape for a species and depends on the intrinsic characteristics of the species and the landscape [79,80]. Indeed, the same landscape could be used very differently from one animal to another. Thus, a well-connected territory should allow various kind of species movements such as daily movement, large-scale migration and dispersion, ultimately permitting gene flow across populations. Hence, the functional connectivity was studied by the identification of both the corridors and the areas constraining species’ movements for three animal species as well as by the mapping of light pollution, an essential factor for nocturnal species’ movements.

2.5.1. Combined Connectivity and Corridors

Three species with various spatial behaviors and using different habitat types were selected for these analyses: Cervus elaphus L. (red deer), Capreolus capreolus L. (roe deer) and Lepus europaeus P. (Brown hare). Two maps representing the global connectivity of the landscape as well as the constraining areas were modeled for each species. They are both based on two main inputs: the species reservoirs (or core areas) and a resistance matrix. Potential reservoirs of wild populations were calculated using species’ resistance and connectivity maps of the Greater Geneva region [81,82]. These maps were created using species’ habitat requirements and expert knowledge based on the LULC information of the territory. Ecological barriers to species’ movements such as buildings and fenced highways are taken into account, which means the resulting model excludes portions of the territory that are not used or accessible to the species.

The first map represents the energetic cost of crossing the habitats located on the animal’s path and thus the probability for it to move across the landscape. Cumulative costs were calculated using a workflow in GIS software between the species’ reservoirs, using a matrix allocating a resistance value to each LULC category of the territory. This method corresponds to the surface generalization of least-cost path models [82]. Hence, suitable habitats close to the reservoirs of the considered species’ population are easier to cross because of its low energetic cost, while unsuitable habitats are more energy intensive.

The second connectivity map emphasizes the areas constrained by the urban occupation. In other words, it represents how few or many alternative ways are available for wildlife to move from one point to another. CircuitScape (v0.1.0) software was used to compute this map based on the circuit theory [81,82,83]. The preservation of these constrained corridors ensures the connectivity of the landscape, even for the urbanized areas.

2.5.2. Light Pollution

The spatial organization of a landscape could be used differently depending on the animal. The spatial behaviors of nocturnal species do not only depend on the landscape’s structure but also on the artificial light of urban areas. Indeed, they need dark spaces to carry out their movements as well as other activities. Thus, the identification of light pollution allows the preservation of areas that are especially shaded and dark during the night and the identification of areas highly polluted by light. To do so, urban areas are transformed into light emitting spaces and the model is adjusted according to the altitude and the presence of forests or water bodies. The map was modeled by F. Tapissier in 2016 [84], before the current restrictions in the use of electricity and urban lighting. Many villages and urban areas now drastically reduce their nocturnal light, and the current map might over-represent current light pollution in the study area.

2.6. Pillar 4: Landscape Structure

Landscape structure, or structural connectivity, corresponds to the spatial arrangement of its LULC categories. The distribution and physical organization of the (semi-)natural habitats were assessed in the territory in order to identify areas with a high interest in conservation, based on five indices: the fragmentation (or the continuity) of natural areas, the soil permeability, the naturality of habitats, the diversity of (semi-)natural habitats and the identification of core natural areas.

2.6.1. Fragmentation

Natural habitats are considered fragmented when their distribution is discontinuous and patches are separated by ecological barriers, mostly anthropic land cover types and transportation networks [85]. The fragmented habitats have a reduced availability for species especially in an urban context. Indeed, sound, odors or light might prevent certain species from living at the margins of their natural habitat if it is surrounded by human-made infrastructures [86]. Furthermore, connected habitats favor species’ movements and migrations. One method to model and map habitat fragmentation is to calculate the MESH size that corresponds to the probability of two randomly picked pixels belonging to the same habitat patch [87]. To do so, the LULC map of the territory was transformed into two binary values, 1 for ecological anthropic barriers and 0 for (semi-)natural habitats. This raster was then used as input in the software Fragstat (v4.2) [88] and the “Effective mesh size (MESH)” program in the window “Aggregation” of the “Class metrics” category was selected. The moving window sampling strategy was used with a round radius of 200 m and a maximum of 50% border with no data. The resulting map was then modified to show the continuity of the natural habitats, which is the exact opposite of the fragmentation.

2.6.2. Soil Permeability

In an urban context, the environment is mostly impermeable, preventing water from being absorbed in the soil which increases the risk of flooding. This impermeability is mostly due to the use of concrete. Thus, saving permeable habitats is a nature-based solution to mitigate the effects of extreme weather events, maintain ecological functions linked to the water cycle and favor soil biodiversity. To map the permeability of the study area, the categories “highways”, “road” and “dense urban areas” were considered as impermeable while the other categories were permeable. This permeability layer favors the conservation of natural habitats in opposition to highly anthropic LULC categories.

2.6.3. Naturality

Naturality corresponds to the ecological quality of a habitat. A very anthropic LULC category would have a low naturality while a highly diverse, well-managed habitat with low disturbance would have a high naturality. This index allows ranking of the habitats according to their intrinsic quality in the study area. Using experts’ knowledge, all LULC categories were assigned a value between 1, corresponding to a very low naturality for urban areas, and 5, for the most interesting habitats. Then, a spatial focal statistic was applied to the raster using a 200 m radius to smooth the values and avoid class boundary effects. The map was exported at 25 m resolution.

2.6.4. Diversity of Natural Habitats

A high diversity of natural habitats favors a high species richness in a territory, especially when these habitats are equally distributed. Thus, the Shannon index [89] was used to calculate the diversity of natural habitats in the study area. To do so, the LULC categories were aggregated in seven classes based on their similarity, without considering dense urban areas and transportation. The classes are the following: meadows, lightly urban areas, natural cliffs and rocks, disturbed vegetation, forests, agriculture and wetlands. The model was run using Fragstat (v4.2), selecting “Shannon’s diversity index (SHDI)” in the “Diversity” window of the “Landscape metrics” class. The moving window sampling strategy was used with a round radius of 200 m and a maximum of 50% border/no data.

2.6.5. Core Areas

This input is complementary to the diversity of natural habitats because it identifies patches of habitat that are large enough to have a central core area free from the influence of neighboring habitats (edge effect). Indeed, some species need large areas of the same habitat to thrive. However, the influence of the neighboring habitats varies depending on their naturality and intrinsic characteristics. For example, an anthropic LULC category has a strong influence that could penetrate deeper into the natural habitat in the form of olfactive, chemical or light pollution. On the other hand, natural habitats have a lower influence between themselves, and their characteristics have a lower level of penetration. A core area would then correspond to an area free from any edge effect.

The modeling of the distribution of core areas necessitates two tables. The first one aggregates the habitats in categories based on their similarity. The second is a penetration matrix linking each category with each other with a distance of influence. The distance must be written using the raster’s metric system and be a multiple of its resolution. Here, the resolution of the raster is 25 m squared, so only multiples of 25 are accepted in the penetration matrix. The values used were based on experts’ knowledge and calibration tests. Then, core areas were mapped with Fragstat (v4.2) using “Core Area Median” in the “Core Area” window of the “Landscape Metrics” class. The moving window sampling strategy was used with a round radius of 200 m and a maximum of 50% border/no data.

2.7. Spatial Conservation Prioritization

Spatial conservation prioritization allows identification and mapping of the optimal compromise between the layers used as inputs, considering their allocated weight and other settings. We used the additive benefit function of Zonation 5 [43] to prioritize each pixel of the study area. The process starts by ranking all the pixels of the study area according to their ecological interest, then iteratively removes cells with the smallest marginal loss in terms of conservation value [43]. The additive benefit function prioritizes areas where many inputs have a high ecological value, thus selecting pixels with high richness over those with rare features. The resulting map is a raster in which pixels’ rank ranges between 0 for low conservation value and 1 for the relative highest ecological interest. Inputs should be used with the same logic in their pixels’ value, by which the most interesting areas that should be conserved should have a high value while low-quality areas have a low value.

The previously mentioned inputs were classified into four classes corresponding to the four pillars and attributed a weight depending on their quality and capacity to identify highly relevant areas for conservation. In total, 2437 inputs were used in the prioritization process, mostly species habitat suitability maps. The details of inputs and their associated weight are found in

Table 2. Summary of all the inputs used in the prioritization process in Zonation 5 and their associated weight.

Pillars	Inputs	Weight/Input	Weight/Group	Weight/Pillar
Species and habitat distributions	Plant species (n = 1480)	0.0116171	25	100
	Red-listed plant species (n = 336)	0.0232342
	Animal species (n = 450)	0.0347222	25
	Red-listed animal species (n = 135)	0.0694444
	Habitats (n = 19)	2.631578947	50
Ecosystem services	Areas for pollinators	4	20	20
	Carbon storage	4
	Nutrient delivery ratio	4
	Sediment delivery ratio	4
	Leaf area index	4
Connectivity	Global connectivity (n = 3)	5	40	80
	Constrained corridors (n = 3)	5
	Light pollution	10
	Fragmentation	8	40
	Permeability	8
	Naturality	8
	Diversity of natural habitats	8
	Core areas	8

. The weights were established empirically by several trials that were discussed among the authors and with the stakeholder group.

The study area is composed of three administrative entities: the French territory named “Pôle Metropolitain”, the Canton of Geneva and the District of Nyon in the Canton of Vaud. We compared the output of an analysis covering the entire region, with an analysis combining the results of three subregions. According to the stakeholders, the final GI network should cover 30% of each of these three regions, even though their size and ecological quality differ. To do so, the prioritization was made over the whole territory but the identification of the best 30% has been made on the separate entities. As the Canton of Geneva has already had a GI identified in a previous work with the same methodology, it was integrated in the final results to preserve the coherence with the previous analysis.

3. Results

3.1. Identification of the Habitats

The LULC map presented in Figure 3 consists of a raster of 25 m resolution and 27 categories. It represents broad natural and semi-natural habitats as well as various land uses and urban environments. Having access to a homogenous, precise and well-covering LULC map is fundamental for the analysis carried in this work but also for having a common baseline between stakeholders and the various administrative entities, especially in a cross-border territory.

3.2. Maps Used as Inputs to the Final Prioritization

The inputs used in the prioritization process represent interesting results by themselves because some areas might not be selected in the final GI while representing a high interest for one of the inputs. The 2437 maps calculated and modeled were aggregated into 41 maps, available for consultation by local authorities and planners in the following link: https://sitv-qual.geneve.ch/portal/apps/webappviewer/index.html?id=f5865a0162d64efb8434644372c022c6 (accessed on 1 July 2023).

All the input variables and the outputs of the analysis are available on the Yareta system for scientific archives (Supplementary Material S1).

3.3. Biodiversity Diagnosis and GI

The main output of Zonation is the biodiversity diagnosis of the study area and is presented in Figure 4. This map represents the ranking of all pixels according to their ecological interest based on the inputs. The mountainous regions located on the north-eastern and the south-western parts of the territory seem of particular interest for conservation, which was expected because the anthropic pressure is lower in these areas compared to the lowlands, located at the center of the map. Urban centers are highly visible, especially the agglomeration of Geneva located in the middle of the study area where the conservation value of the cells is low. Forested areas in the lowlands also seem of high interest.

From this covering map, the top 30% pixels from the French administrative entity as well as the Swiss District of Nyon were independently extracted and then merge together with the GI of the Canton of Geneva identified in a previous work. This methodology, presented in Figure 5, allows identification of the GI for each administrative entity as well as for the whole territory to fit all political agendas.

The GI identified is visible in Figure 6. It is mostly made up of large patches well connected with many smaller isolated patches orbiting around the main core areas. Most of the GI is covered by closed forests (35%) and natural meadows (25%), while cultivated and diffuse urban areas both cover 16% of the GI, wetlands only 5% and open forests 3%. Closed forests thus represent more than one-third of the GI which is not surprising because they cover a similar portion of the whole study area, which means that they are not over-represented. However, natural meadows only cover 14% of the study area which means that they are over-represented in the GI, as are the wetlands only covering 1.5% of the study area. These habitats are thus highly important for biodiversity conservation in the region. On the other hand, cultivated areas are under-represented because they cover more than one-fifth of the study area. Another way of verifying which habitat is especially important for conservation is to look at the proportions of these habitats integrated into the GI. As previously mentioned, slightly more than one-quarter of all closed forests is integrated in the GI while this proportion rises to 53%, 97% and even 100% of all natural meadows, open forests and wetlands, respectively. This is a clear signal that these habitats are crucial elements of a functional, effective and reliable network of conservation.

Figure 6 also shows the proportion of the GI that is already protected. This result demonstrates that there is still a lot of effort to be made to preserve the areas with a high biodiversity value in the region for the benefit of future generations in the context of steady economic and population growth. As the best identified 30% of the territory will probably not obtain a legal status of protection, other efficient area-based measures (OECM) should be considered for their preservation as proposed by the IUCN [90].

The GI was initially identified by selecting the best 30% of the main Zonation output pixels over the whole study area, but the GI cover over the three territorial entities was disequilibrated because the main natural areas are found in the mountains located on the French side. This initial version would represent the best possible areas to conserve over the whole study area as there is no modification brought to the result of the prioritization process. However, with the methodology presented here, each entity has its own 30% GI cover in its territory. This is more equitable for the administrative entities and valuable for spatial conservation planning in order to reach the objective of dedicating 30% of each territory to nature’s conservation. The aggregated version of GI in the study area redistributed cells mainly from the patches in the mountainous areas to the more urban and cultivated areas in the lowlands as presented in Figure 7. Interestingly, most of the pixels that were removed from the initial version are either isolated or located at the margins of the main patches. Thus, the GI distribution from the aggregated version does not fundamentally change compared to the first version. This result means that the main patches identified in the initial “optimal” version are still included in the aggregated version. Thus, the conservation value of the aggregated GI network is still similar to the initial version, which is verified by the high spatial overlap of 88.1% between the two versions.

4. Discussion

4.1. Selection of Inputs

The inputs used in this methodology were selected based on their complementarity to avoid collinearity as much as possible, their representativity of the natural processes and their ease to be reproduced and explained to the authorities in charge of implementing the results. The final selection is the result of the separate assessment of each pillar in order to test the most reliable inputs for the local context and with the data available [26,40,52].

The distribution of natural habitats is key information for the modeling of the majority of inputs. Having access to a high-resolution LULC map with many detailed categories is fundamental for this work and its quality determines the accuracy of the resulting GI. Identifying many natural habitats allows consideration of each one of them in the prioritization process and thus ensures their representativity and conservation in the GI. However, for most of the inputs and species distribution models, eight categories are sufficient [52]. All species of plant and animal for which enough precise occurrences were available were modeled, with no distinction for their native status. However, a higher weight was given to species considered vulnerable, endangered or critically endangered in one of the red lists of the three territories. The inclusion of as many species as possible in the prioritization process creates a more reliable GI network by integrating all available distributions.

The selection of the ecosystem services was mainly based on their modeling availability in the software InVEST coupled with the perspective of conserving biodiversity. Many other ecosystems services were available but were linked with the production of resources and energy, or with the cultural value of the landscape, which might have deleterious outcomes for nature’s conservation. Although conserving highly biologically diverse areas usually has a positive influence on the preservation of qualitative and quantitative ecosystem service supplies, the opposite is not true [42,91]. The modeling of ecosystem services highly depends on values given to the settings and to the biophysical tables, as well as the quality of input maps. Specific values for the LULC categories used in the biophysical tables are mostly impossible to find in the literature, especially when they have to be adjusted to the local context. Furthermore, many settings require relative values, which depend on the characteristics of the inputs and of the territory. Thus, experts’ knowledge is highly valuable in this type of work. They should examine the inputs with care and verify the credibility of the resulting maps to iteratively calibrate the settings and tables.

The functional connectivity of the landscape for the species Cervus elaphus L. was calculated with the help of GPS trackers placed on several individuals [81,92]. The data collected allowed to define a resistance matrix based on the species’ habitat preferences, and thus represent a fundamental input for modeling the global connectivity and constrained zones. It has not been carried out for the two other species where inputs were based on experts’ knowledge [82]. It would be highly interesting to use GPS trackers for more individuals and species, but this method is expensive and not always applicable, especially for small animals. However, to better understand the functional connectivity of the territory, more animals could be studied such as amphibians, reptiles or small mammals, although large mammals’ results might serve as an umbrella for other species. Even though small animals’ movements do not reach the full extent of the study area, assessing their connectivity might also result in identifying corridors at a finer scale, resulting in a better representativeness of animal connectivity in the prioritization process. This assessment, however, is more challenging because of the lack of fine-scale data and its computational requirements.

4.2. Prioritization

The prioritization process was carried out over the whole territory but the best 30% was selected for each administrative entity. One limit of this approach could be the loss of connection at the edges of the three territorial entities, because the methodology does not consider that one patch should be entirely selected if its area is located on both sides of a border. This would be an important loss for the global structural and functional connectivity of the network as borders have no impact on species movements. However, this problem was not observed here, which could be explained by the fact that most of the borders are in the urbanized lowlands. In these areas, natural habitats with high ecological values are not so common and, thus, are easily identified by the methodology. This implies that the cross-border patches of natural habitats are selected on both sides of the border. Another reason that could explain this pattern is that the habitats located in the highlands and the lowlands are different. For example, deciduous forests are mainly found in the lowlands and coniferous forests in the highlands. Thus, the prioritization process would preferentially select patches of the same habitat, especially if it is not found anywhere else in the territory. It means that having access to a detailed map of the LULC information including many natural habitats might prevent the final network from being too much impacted by the spatial limits of the analysis.

The GI identified here covers the most interesting areas for nature’s conservation over the whole study area, considering the local context and administrative entities. However, it does not imply that areas located outside of the GI should not be considered for protection as well. This result is an overview of an optimal protection network to conserve all the aspects of the biodiversity and should be seen as a common aim for the territory. At a smaller scale, the conservation of natural habitats is fundamental because they could host important functions and diversity at this scale that are not represented or that are repetitive at the regional scale of the study area. Nevertheless, identifying and protecting a network of natural habitats is part of the solution to halt biodiversity loss but should be complemented, for example by lowering our impacts outside protected areas. Indeed, it is preferable to conserve, protect and restore more than not enough.

4.3. Perspectives

The proposed GI is still theoretical, and its effective implementation raises questions, especially regarding the legislation and the inclusion of private lands. There are already many protected areas with various appellations and legal basis in this cross-border area. Integrating conservation areas into the GI is an interesting idea, but the selection of which types of conservation area are integrated or not should be discussed among the stakeholders as some of them do not have strict legal basis. However, the integration of such large patches of protected areas might change the final distribution of the GI due to the prioritization process because most protected lands are made up of only a few natural habitats. Another solution would be to use the GI as informative data for spatial and urban planning to avoid the destruction of areas of high ecological interest, but the effective positive impact of GI on nature’s protection might be lower.

The ponderation used in the prioritization process was selected according to the value and interest attributed to the inputs by a group of experts. However, a more inclusive approach might be interesting, using, for example, the “best–worst method” based on the multiple criteria decision-making process [93]. This method allows sorting of the inputs from the best to the worst one and vice versa based on their conservation value. At the end of the process, a weight is automatically calculated and can be attributed to each input and each pillar. This approach is interesting but should be used with caution. Indeed, being able to classify the interest given to maps representing ecological processes is a complicated task and necessitates knowledge about the study area, its ecology and ecosystems, as well as about the method used to calculate the inputs which participates in evaluating their quality. Thus, the inclusion of the stakeholders is highly important, and they should be part of the process and decisions to model and map the inputs in order to fully understand their intrinsic meaning.

5. Conclusions

Most of the complex analyses presented in this work rely on the cross-border map of natural habitats that was first created for the Greater Geneva region. Without this common input the rest of the analyses would have been made difficult. All the environmental data (topography, water and climate, species, etc.) are by nature crossing borders. Most of the difficulties faced during this work are linked to the cross-border situation of the study area which means that a lot of work is still needed to share common data, classifications and methods at a larger scale to allow such analysis. However, this article proves that innovative insights, positive for the regional biodiversity, could emerge when the scientific community and the stakeholders from the different administrative entities work together. Indeed, the resulting maps from this article were immediately transferred to the land use planners in charge of developing ambitious visions of the “Greater Geneva” territory for 2050 in alignment with 10 objectives of ecological transition as recently agreed and signed by local authorities [41].

This article has shown that the territory of “Greater Geneva” could be classified from its best to its worse pixel, in terms of ecological value, using a prioritization process based on the use of many inputs grouped into four main pillars. This approach is very useful to identify the best areas that should be considered as hotspots of biodiversity and ecosystem services but also to identify coldspots that can be good candidates for potential urban developments or for ecological restoration. The results have shown that selecting the best 30% of each administrative entity of the study area was better accepted by stakeholders and did not fundamentally change the quality of the GI. The share of the identified GI that is already under protection is relatively low, demonstrating that much effort would be needed to reach the international target of 30% by 2030. However, the use of OECM as advocated by the IUCN could be a good solution to provide a status to the newly identified areas.

The GI identified in this work is the result of several years of research on the pillars and inputs on the same territory. The study area is well prospected and working together with naturalists, experts and the authorities in charge of the ecological conservation of the territory represents an asset to tune the settings to fit experts’ knowledge, field observations and political agendas. The methodology developed here is adaptable to other territories, depending on the data availability, and we believe it allows selection of the most interesting areas from an ecological perspective representing all aspects of biodiversity, ultimately resulting in a highly relevant and robust GI.