Spatial Planning of Green Infrastructure for Mitigation and Adaptation to Climate Change at a Regional Scale

Green infrastructure has acquired greater importance in recent years in relation to climate change adaptation. Green infrastructure planning has been identified as a new and innovative means of land planning that can contribute to preventing the impacts of climate change. However, this has been explored more thoroughly in urban areas than at the regional scale. The present study proposes a methodology including multi-criteria evaluation techniques for assessing the ESS involved in the fight against climate change and for the spatial planning of multifunctional green infrastructure areas based on the results of this assessment. Application of the methodology for green infrastructure planning aimed at confronting climate change at landscape level in the region of Galicia (NW Spain) successfully delimited multifunctional green infrastructure zones. Results show that delimited zones have a higher provision potential for more ESS than protected natural areas and areas that are not part of the green infrastructure.


Introduction
The IPCC predictions for climate change (CC) in Europe foresee hotter and drier summers and alterations in precipitation regimes, with an increase in the frequency of heavy storms and rainfall events [1]. These changes will affect ecosystem dynamics and thus the composition of the species the ecosystems host. As a consequence, CC may exert negative impacts on biodiversity [1], as recognized in the EU biodiversity strategy for 2030 [2], which states that CC and biodiversity loss are intrinsically linked. Biodiversity and the provision of other ecosystem services (ESS) may, in turn, play an important role in preventing the negative effects of CC [3]. In this line, the EU biodiversity strategy [2] states that "nature is a vital ally in the fight against CC" and refers to nature based solutions and green infrastructure (GI) as key factors in mitigating the negative effects of CC. Consequently, the EU Commission communication "Green Infrastructure: Enhancing Europe's Natural Capital" [4], establishes the basis for the creation of GI in EU member states, in order to tackle biodiversity loss. This endeavor requires methodologies that can be applied at a regional scale to help policymakers and planners to delimit GI elements. In this paper, we will focus on the development of a methodology for delimiting GI aimed at helping to mitigate the negative effects of CC on biodiversity.
The concept of GI was first proposed by Rosenberg [5], who considered urban parks an extension of the infrastructure of cities, providing functions in addition to recreation, such as storm water infiltration. The notion of GI has since evolved and many definitions have been proposed [6]. The idea most frequently shared by these definitions is that GI is a network of green areas interconnected by are based on multi-criteria decision analysis [41,42] or consist of algorithms of heuristic optimization designed for different purposes and that are applied to delineate GI [23,43,44].
Methods of multi-criteria decision analysis include that proposed by Liquete et al. [42], which combines the provision potential of different ESS to identify GI areas for conservation (those with a high combined ESS provision potential) and GI areas for restoration (those with a moderate ESS provision potential). Kopperoinen et al. [40] also combined maps of the provision potential of several ESS to identify GI multifunctional areas.
In addition to maximizing the provision potential of several ESS, heuristic methods can incorporate criteria for delimiting GI that cannot be considered in methods based on multi-criteria decision analysis. This applies to criteria such as the size and compactness of the delimited areas or their proximity to other elements of the GI. These methods can therefore delimit larger and more regular GI areas [45] that reduce edge effects [28] and have a greater potential to provide a wider variety of ESS [37]. On the other hand, heuristic optimization methods require more computation time and specific knowledge for fitting the optimization criteria. In addition, they produce an average optimum for the criteria considered and therefore do not grant maximization of ESS provision potential in the areas delimited [45]. The available heuristic methods that can be used for GI planning include spatial conservation prioritization techniques for planning natural spaces [23,44] and land cover optimization algorithms for land planning [43]. In addition, the Marxan software was designed for reserve selection in conservation planning and has been used by Schröter and Remme [45] to delineate ESS hotspots and by Vallecillo et al. [44] to delineate GI areas under different scenarios considering several ESS.
None of the aforementioned methods delineate multifunctional GI areas that would provide a variety of ESS with a specific purpose, such as mitigating the impacts of CC on biodiversity and the population. Moreover, none of the methods delineate GI areas by considering how they will be integrated with other GI elements such as core areas and corridors. These two elements are key to developing a methodology for spatial planning of GI, as areas must be delimited according to their function and by taking into account how they complement other areas in zoning plans.
In the present study, we propose a methodology for delimiting multifunctional zones that will act as buffers to counteract the impacts of CC in core areas and corridors in GI. The method has been tested by applying it to GI zoning in the region of Galicia (NW Spain). With this aim, the ESS involved in CC mitigation and adaptation were identified and their provision potential was mapped for the whole region. A multi-criteria evaluation method was used to combine the resulting ESS provision potential maps to produce a suitability map for multifunctional areas. The cells with the highest suitability values were selected, and the largest highly suitable cell patches adjacent to core areas and corridors in the GI were then delineated as multifunctional buffer zones. The resulting zones were analyzed to check whether they meet the criteria of having a high potential for providing several ESS and serving as buffer zones to prevent the impacts of CC in core areas and corridors in the GI. The methodology successfully delimited multifunctional GI elements that can counteract the impacts of CC since results show that they have a high potential for providing more CC related ESS than GI corridors and core areas as well as zones outside the GI. The delimited areas were found to be rather irregular due to landscape heterogeneity in the region. Suggestions are made for ways of incorporating landscape structure for delineating GI, in order to produce more regular GI elements and thus increase ESS provision.
In the following sections, the study area is introduced together with the criteria that are followed to select the CC related ESS to map. Then, the methodologies that were used to map the provision potential of the selected ESS and to delimit the GI multifunctional buffer zones are presented. The obtained ESS provision potential maps and the delineated multifunctional buffer zones are analyzed in the results section. With this aim, zonal statistics were used to estimate which land covers have the highest ESS provision potential and whether the delineated GI multifunctional buffer zones have a ESS provision potential higher than other areas. Landscape statistics were also used to assess whether the delineated Sustainability 2020, 12, 10525 4 of 22 GI zones are compact and located close to the areas that demand the ESS they produce. Finally, the main conclusions from the results analysis are presented.

Study Area
The region of Galicia, in north-western Spain (Figure 1), occupies an area of approximately 2.9 million hectares and has a population of around 2.7 million. Some 60% of the inhabitants of Galicia own land. Each landowner has on average 1.7 ha of land divided in 7 plots, resulting in an average plot size of 0.25 ha [46]. Around 30% of the population settlements in Spain are in Galicia [47], although the surface area of the region only accounts for 5.7% of the total surface area of the country [48]. The Galician population is disperse, but most of the population is concentrated on the Atlantic shore of the region. Inland Galicia is mainly rural and is currently suffering a population decline, leading to abandonment of agriculture and afforestation of agricultural land with fast growing exotic species [49]. This situation is leading to deterioration of ecosystems and therefore the biodiversity the ecosystems host [50] and the services they provide.
although the surface area of the region only accounts for 5.7% of the total surface area of the country [48]. The Galician population is disperse, but most of the population is concentrated on the Atlantic shore of the region. Inland Galicia is mainly rural and is currently suffering a population decline, leading to abandonment of agriculture and afforestation of agricultural land with fast growing exotic species [49]. This situation is leading to deterioration of ecosystems and therefore the biodiversity the ecosystems host [50] and the services they provide.
Galician soils are mainly sandy and have a high organic matter content, due to the mild and rainy oceanic climate that predominates in most of the region [51]. Average summer temperatures are rising and summer droughts are more frequent due to CC [52]. This has led to an increase in the frequency of wildfires in the region [53]. In the future, soil organic matter may evolve faster due to higher temperatures and therefore its capacity to retain water and nutrients will be reduced, exacerbating the effect of drought and leading to changes in the species composition of ecosystems. This adds to the deterioration of ecosystems due to rural population decline and is further exacerbated by CC.
The aforementioned trends increase the need for GI planning that will help to mitigate the impact of CC on natural values. The zoning of GI is rather complicated in the region due to land cover fragmentation, fostered by smallholdings and a varied terrain composed of hills and small valleys.

Selection of Ecosystem Services
Delineation of multifunctional buffer zones in GI is based on the capacity of an area to provide ESS involved in the reduction of CC impacts on biodiversity. In order to identify these ESS, the potential impacts of CC in the study area were identified by reviewing the IPCC reports.
According to the 5th assessment report of the IPCC [1], in the coming decades, CC will have the following consequences in Europe:

-
More frequent extreme precipitation events in northern and continental Europe, which will increase erosive processes and floods. -Increased drought and wildfire risk in the south of Europe. -Leaching of soil nutrients and pollutants in northern and continental Europe due to an increase in winter precipitation, which will lead to water pollution, poorer water quality, and reduced soil fertility. -Increased incidence of forest and agriculture pests and diseases. Galician soils are mainly sandy and have a high organic matter content, due to the mild and rainy oceanic climate that predominates in most of the region [51]. Average summer temperatures are rising and summer droughts are more frequent due to CC [52]. This has led to an increase in the frequency of wildfires in the region [53]. In the future, soil organic matter may evolve faster due to higher temperatures and therefore its capacity to retain water and nutrients will be reduced, exacerbating the effect of drought and leading to changes in the species composition of ecosystems. This adds to the deterioration of ecosystems due to rural population decline and is further exacerbated by CC.
The aforementioned trends increase the need for GI planning that will help to mitigate the impact of CC on natural values. The zoning of GI is rather complicated in the region due to land cover fragmentation, fostered by smallholdings and a varied terrain composed of hills and small valleys.

Selection of Ecosystem Services
Delineation of multifunctional buffer zones in GI is based on the capacity of an area to provide ESS involved in the reduction of CC impacts on biodiversity. In order to identify these ESS, the potential impacts of CC in the study area were identified by reviewing the IPCC reports.
According to the 5th assessment report of the IPCC [1], in the coming decades, CC will have the following consequences in Europe: More frequent extreme precipitation events in northern and continental Europe, which will increase erosive processes and floods. -Increased drought and wildfire risk in the south of Europe. -Leaching of soil nutrients and pollutants in northern and continental Europe due to an increase in winter precipitation, which will lead to water pollution, poorer water quality, and reduced soil fertility. -Increased incidence of forest and agriculture pests and diseases. -Changes in ecosystems and in the composition of animal and plant populations. In southern Europe, these changes may lead to a decrease in biodiversity and to homogenization of populations. As a consequence, ESS such as pollination, pest control, and biodiversity may be affected.
The ESS involved in reducing the previously identified impacts of CC were selected by consulting version 5.1 of the Common International Classification of Ecosystem Services (CICES 5.1) developed by the European Environment Agency [54]. We used this classification as it is widely applied worldwide.
According to CICES 5.1, the following ESS will contribute most to counteracting the consequences of CC in the study area: Long-term carbon storage: Capacity of ecosystems to capture and store carbon for more than 150 years.

Mapping of ESS Provision Potential
The methods used for mapping the ESS depend on the data available. In most cases, the ESS were mapped by using land-cover as a proxy. Additional biophysical variables were used when available. All variables considered were combined by using a weighted sum model to estimate the ESS provision potential. When enough data became available, more complex regression models were used [29].

Materials
As previously mentioned, in most cases, land cover was used as the main variable to map the ESS. The most accurate and up-to-date data on land cover for the study region is the Spanish Information System on Land Occupation (SIOSE) for 2014. The scale of these data (1:25,000) is of high enough resolution for mapping the ESS at the case study scale.
The SIOSE includes numerous mixed land cover categories composed of the percentages of several land cover types. In order to be able to use the data, these mixed categories were assigned the land cover type occupying the highest proportion of the surface area. In those cases where two types of land cover occupied the same proportion of land, a mixed category was created. Many land cover categories were also combined into a single category in order to simplify the data and make them easier to handle. The following SIOSE categories were finally considered for mapping: bare rock, bodies of continental water, bodies of salt water, cliffs, artificial land cover, coniferous forest, crops and fields, crops, deciduous forest, eucalyptus plantations, eucalyptus and coniferous, wetlands, sport facilities, shrub land, shrubs and trees, shrubs and rocks, mixed forests, crops and shrubs, crops and urban, crops and trees, fields, beaches, afforested land, transportation infrastructure, landfill, vine yards and woody crops, mineral extraction sites, burnt areas, urban areas, urban green areas.
Additional data used for ESS mapping included the following:

Methods
We present below in more detail the methods developed in this study to map the provision potential of the following ESS: Regulation of hydrological flow (

Regulation of the Hydrological Cycle
The potential of ecosystems to regulate hydrological flow was estimated by considering the capacity of the land cover and soil to retain and infiltrate water [58]. The capacity of land cover to retain water was determined by assigning a weight between 0 and 1 to each SIOSE land cover type, where 1 indicates a high capacity to retain water and 0 indicates no capacity to retain water ( Table 1). The weights were estimated on the basis of those used in Burkhard et al. [59] to map flood protection and other ESS in a region of Germany by assigning weights to the land cover types included in the CORINE (European Union program for COORdination of INformation on the Environment) Land Cover database. SIOSE feeds data to the CORINE land cover database and therefore there is a correspondence between SIOSE and CORINE land cover categories, thus enabling use of the aforementioned weights as reference values. The capacity of soil to retain water was estimated considering the Topographic Wetness Index (TWI) and the percentage of sand in the soil [60]. The TWI was obtained from the DTM. The map of soil sand percentage was resampled to a resolution of 25 × 25 m using a bilinear interpolation method. Both the TWI map and the sand percentage map were normalized between 0 and 1 using a linear normalization method (Equation (1)).
Flood plains with a return period of 100 years were also taken into account due to their key role in regulating the hydrological cycle. Therefore, a binary raster map of resolution 25 × 25 m was produced where the flood plains were assigned a value of 1 and the other areas a value of 0.
All of the aforementioned maps were combined by using Equation (2) to obtain the potential for regulation of the hydrological cycle.
ES provision potential = normalized capacity of land cover to retain water × 0.35 Land cover and TWI were considered the two main factors that control water retention by the vegetation and the terrain. This is the reason why they have been given the higher weights. Percentage of sand only accounts for the drainage capacity of soils and thus runoff reduction; yet runoff will be mainly reduced by vegetation. In the case of flood plains, they overlap areas with a high TWI. therefore, these two factors were given a lower weight.

Fire Protection
This ecosystem service is related to the capacity of ecosystems to prevent the spread of forest fires, which is determined by the fire severity. The severity of a wild fire is conditioned by the ecosystem biomass (type, volume, and structure) and its humidity content [61].
The ecosystem capacity to prevent fire expansion due to biomass was estimated using SIOSE land cover as a proxy. Therefore, weights between 0 and 1 were assigned to land cover types, where 1 indicates land cover with the lowest content of biomass and thus the highest capacity to prevent fire. These weights were determined by consulting the studies of Fernandes, Proença et al., and González et al. [61][62][63] (Table 1).
Biomass humidity is determined considering precipitation and temperature [61]. Solar irradiance was used as a proxy for temperature as this variable is more representative of temperature at a fine scale [64] than the available temperature data. Maps of annual average precipitation were used to yield precipitation data. Both solar irradiance and precipitation maps were resampled to a 25 × 25 m resolution through bilinear interpolation, and they were normalized to between 0 and 1. The maps were normalized using Equation (1) for precipitation and Equation (3) for irradiation.
As soil moisture also influences the fuel humidity [65], the TWI was also taken into account in calculating the fire protection potential. Slope and wind also influence the speed of fire spread, which determines fire severity and thus the capacity of ecosystems to prevent fire expansion [66]. Slope is already considered in the TWI as it is one of the parameters used to calculate this index. However, wind data with a sufficient resolution for the scale of this study were not available. Therefore, only the TWI was considered and was normalized to values between 0 and 1 by using Equation (1).
All the variables were combined using Equation (4) to determine the potential of ecosystems for fire protection.
ES provision potential = normalized biomass × 0.4 + normalized TWI × 0.3 + normalized irradiance × 0.2 + normalized precipitation × 0.1 Biomass is the most determinant factor; this is why it has been given the highest weight. The moisture content of biomass is the second most important and is determined by water stored Sustainability 2020, 12, 10525 9 of 22 in soil and temperature. As it is difficult to obtain data on soil water content in the driest months it was decided to use the TWI as a proxy. Therefore, it was given the second highest factor. Irradiance was considered as a proxy of temperature at a very local scale. This is why it was given the third highest weight. Precipitation was used to account for humidity due to climatic conditions. However, there were only data at a global scale and it was not possible to estimate the water deficit in drier summers. This is why this factor has been given the lowest weight.

Pollination
The methodology used to map the provision potential of this ecosystem service is based on that reported in InVest [67], which can take into account several pollinator species and uses land cover maps as input data. In this study, we only considered bees as pollinator species, due to their importance in crop pollination. Thus, we considered a hypothetical generic bee species to represent the average behavior of the bee species that may be present in the area.
The potential of ecosystems to promote pollination is estimated by multiplying two indices that account for the accessibility of the pollinator to floral resources in an area (FR) and the suitability of land cover for providing nesting sites (HN). FR is calculated using Equation (5): where w sj is the weight that represents the relative production of floral resources for pollinator species s in season j, D (x,x') is the distance in a straight line between locations x and x' within the flight radius of the pollinator, α s is the average flight distance for the pollinators, RA(l,j) is the relative abundance of floral resources in land cover l during season j, and fa(s,j) is the relative feeding activity for the pollinator species s during season j.
As only a generic bee species was considered and seasonal variance in flower resources cannot be taken into account due to a lack of data, w sj and fa(s,j) were omitted from the equation. The average flight distance (α s ) for an average bee species is estimated to be equal to 1 km [68]. RA(l,j) was calculated by assigning weights between 0 and 1 according to land cover potential for providing flower resources. The weights (Table 1) were assigned according to Grundel et al. and Carré et al. [69,70]. The parameter related to nesting sites (HN) was calculated using Equation (6): where N(l,s) is an index that determines the capacity of land cover l to provide nesting for pollinator species s. If a pollinator species has several nesting preferences, N(l,s), the index assigned to cell x is that maximizing the preference for pollinator species s. As we considered a generic bee species, the species preference for one type of nest or other is not applicable. Therefore, land cover types were weighted according to their capacity to provide any kind of nest. Therefore, the weights for nesting suitability associated with each land cover type were used directly as the value of HN. These weights take values between 0 and 1 and were estimated taking into account the papers from Potts et al. and Carré et al. [70,71] (Table 1).

Maintenance of Breeding Populations/Habitats and Pest Control
Both ecosystem services were mapped in a similar way, using as a basis the 2012 Spanish Survey of Terrestrial Species. This survey uses a grid of 10 × 10 km where the number of individuals of each species in each cell is provided.
In the case of the ecosystem service related to maintenance of breeding populations, the abundance of species in each cell of the grid was used as an indicator of the ecosystem capacity to host breeding populations and provide breeding habitats. The data from the Spanish Survey of Terrestrial Species were downscaled to produce a map with a cell size of 25 × 25 m.
Many researchers use logistic regressions to downscale species abundance data to a 1 × 1 km grid [72,73]. Nevertheless, logistic regressions establish the relationship between the probability of occurrence of a binary variable (presence or absence of a species in a grid) and a series of independent variables. In this study, we aimed to downscale the 10 × 10 km data by establishing relationships between a number of independent variables and the number of species present in a cell. Thus, we used a Poisson regression model [74], which is highly suited to determining these type of relationships [75].
The independent variables were selected by considering the variables used by Olivero et al. and Keil et al. [72,73] and selecting those for which there were available data for the study area. Therefore, the following independent variables were considered and calculated for each cell in the species abundance map: Mean average annual precipitation, mean solar irradiance, mean slope, variance of the slope, mean elevation, variance of the elevation, mean distance from motorways, mean distance from rivers, mean distance from population settlements larger than 1000 inhabitants, mean population density, edge density (ED) of land cover, and percentage of land cover (PLAND). For calculation of ED and PLAND, the land cover categories were grouped into 10 land cover types according to their similarity ( Table 2). The ED and PLAND were then calculated for each of these 10 grouped land cover types by using FRAGSTATS [76]. The Poisson regression model was calibrated with a backward stepwise regression method using the "stepwiser" tool of R statistics software. The resulting model was then run in order to downscale the species abundance data, using the variables calculated for the cells of a grid with 1 × 1 km resolution as input. The resulting 1 × 1 km species abundance map was resampled to a cell size of 25 × 25 m using bicubic interpolation. This map represents the capacity of ecosystems to host breeding populations and habitats.
The same methodology was used to calculate the pest control potential. In this case, only the abundance of species that act as pest predators was considered. Civantos et al. [77] used a similar method. Therefore, the pest predators considered were those considered by Civantos et al. [77] and that are present in the study area.

Delimitation of Multifunctional GI Zones
The multifunctional buffer areas delimited using the current methodology are aimed at mitigating CC by means of carbon capture and storage and at reducing CC impacts in the core areas and the corridors. Therefore, delineation of these GI zones was based on the identification of areas with a high potential to provide several ESS involved in the mitigation and adaptation of CC and on the proximity of these areas to core areas and corridors of the GI.
The core areas of the GI correspond to the Natura 2000 network in Galicia. The corridors were delimited by the Institute of Agriculture Biodiversity and Rural Development of the University of Santiago de Compostela (IBADER) by interpreting aerial photographs and using the river network of the region as a basis.
The potential of an area to provide several ESS was determined by taking into account an average ESS provision potential map. This map was obtained by adding the normalized values of the ESS provision potential maps considered (Equation (1)) and dividing the result by the number of maps.
Multifunctional buffer zones were delimited by selecting 20% of the cells of the average ESS provision potential map with the highest values. The area of the patches formed by the selected cells was calculated and those patches of area less than 100 ha were discarded. Selection of patches with a minimum size ensures that the resulting multifunctional areas produce a large number of ESS [37] and can be properly managed.
Finally, from the remaining patches, those adjacent to a corridor or a core area were selected in order to ensure that the multifunctional buffer zones can reduce CC impacts in GI core areas and corridors. The holes inside the multifunctional buffer zones and the gaps between them and the core areas or corridors were added to the multifunctional areas when they met the following conditions: The holes and gaps are smaller than 5 ha or they are bigger than 5 ha but they do not correspond to artificial land cover.

Maps of ESS Provision Potential
Mapping the ESS provision potential (as described above) produced the maps shown in Figure 2. The maps of the provision potential for other ESS considered for delineating the multifunctional buffer areas were provided by García et al. [57]. For a more detailed analysis of the results of the mapping, the SIOSE land cover maps were overlaid on the ESS provision potential maps in order to calculate the average provision potential of each land cover type by using the QGIS zonal statistics tool. This analysis produced information about which are the types of land cover with the highest potential to produce several ESS and thus provide clues on which should be included in the planned GI.
The top five types of land cover in terms of provision potential for each ESS are shown in Table 3. Forest was the land cover type with the highest provision potential for filtering pollutants, erosion control, hydrological flow regulation, and carbon capture. In the case of carbon capture, the forests with the highest provision potential were productive stands of fast-growing species. On the other hand, fire protection followed the opposite behavior and the types of land cover with the highest provision potential were those with lower biomass, such as agricultural land and rocky terrain. Regarding the maintenance of breeding populations and pest control, the types of land cover with the highest provision potential were very similar and were related to coastal areas and wetlands. We can see some inconsistencies in the ESS long-term carbon storage, as bare rock had a high potential. This is due to the low resolution used for the cartography of bogs, which have a very high provision potential for this ESS and are usually located in depressions on mountain tops, where rocky outcrops are frequent. Afforested land also had a high ESS provision potential, which may be due to eucalyptus plantations encroaching on these bogs. Finally, the land cover types with the highest potential for provision of pollination-related ESS were those that provide more nesting sites for bees. These are namely: bare land where nests can be dug, rocks and trunks with cavities, and beehives (in urban areas). Sustainability 2020, 12, x FOR PEER REVIEW 13 of 22   The maps in Figure 2 were used to produce a suitability map for delineating the multifunctional buffer zones, so that they met the objective of maximum ESS provision. As explained in the previous section, this suitability (or average ESS provision potential) was obtained by adding all the normalized ESS provision potential maps and dividing the result by the number of ESS. As in the case of each individual ESS, the resulting average ESS provision potential map was overlaid with the land cover polygons of SIOSE 2014 to extract the average provision potential value for each land cover type. The top three land cover types with the highest average provision potential (Table 4) were native deciduous forest, followed by mixed forest and productive forest. Wetland also displayed a high average provision potential. The land cover type with the highest potential after forests and wetlands was shrubland.

Delineation of Multifunctional Buffer Zones
Multifunctional buffer zones were delineated according to the methodology explained in previous sections, producing the map presented in Figure 3.

Delineation of Multifunctional Buffer Zones
Multifunctional buffer zones were delineated according to the methodology explained in previous sections, producing the map presented in Figure 3. The multifunctional buffer zones delimited are quite irregular, and thus the edge they share with corridors and core areas of the GI is rather limited. This can be seen by calculating the ratio of the edge shared between multifunctional buffer areas and other GI by the total edge of the The multifunctional buffer zones delimited are quite irregular, and thus the edge they share with corridors and core areas of the GI is rather limited. This can be seen by calculating the ratio of the edge shared between multifunctional buffer areas and other GI by the total edge of the multifunctional buffer zones, which is 7.1%.
The irregularity and low adjacency were influenced by the heterogeneous landscape of the region, which is highly fragmented. To demonstrate this, the study area was split into uniform landscape areas considering the 12 Great Landscape Areas of Galicia (Figure 4) defined in the Galician Landscapes Catalogue [78]. The heterogeneity of the landscape in each Great Landscape Area was estimated by calculating the edge density (ED) of the SIOSE land cover types in each area with FRAGSTATS [76]. The irregularity of the delimited multifunctional buffer areas within each landscape area was then measured by calculating the mean fractal dimension (FRAC_MN) with FRAGSTATS [76] ( Table 5). The correlation between the ED and the FRAC_MN was calculated with the Spearman index, yielding a value of ρ = 0.69, indicating a high correlation.   From the map in Figure 3, it can be seen that despite the irregularity of the multifunctional areas, in most cases, the methodology successfully delimited zones close to core areas and corridors, where the ESS are required.
In order to test that the delimited multifunctional buffer zones meet the requirements of having a high potential for producing several ESS, the average ESS provision potential was calculated in the areas of the GI and outside the GI. Comparison of the average ESS provision potential in the GI areas and the rest of the region showed that the average ESS provision potential in delimited multifunctional buffer zones was higher (Table 6).  From the map in Figure 3, it can be seen that despite the irregularity of the multifunctional areas, in most cases, the methodology successfully delimited zones close to core areas and corridors, where the ESS are required.
In order to test that the delimited multifunctional buffer zones meet the requirements of having a high potential for producing several ESS, the average ESS provision potential was calculated in the areas of the GI and outside the GI. Comparison of the average ESS provision potential in the GI areas and the rest of the region showed that the average ESS provision potential in delimited multifunctional buffer zones was higher (Table 6). More detailed analysis of the average provision potential for the ESS considered in each GI area (Table 7) showed that most of the ESS have a higher average provision potential in the delimited multifunctional buffer zones than in the other areas of the GI. This is not the case for breeding populations, pest control, or long-term carbon storage. These 3 ESS were found to be closely correlated with biodiversity. For long-term carbon storage, the highest provision potential corresponded to bogs and wetland, both of which host high biodiversity values. As core areas and corridors were delimited by taking into account priority habitats and biodiversity values, it is not surprising that the average provision potential was higher for these ESS. Nonetheless, the average provision potential of these three ESS in multifunctional buffer areas was higher than for areas outside GI, except for long-term carbon storage, for which similar values were obtained. For fire protection, the average provision potential was similar in GI areas and areas outside GI. Nevertheless, provision of this ESS was related to low biomass land cover types, while the provision of other ESS was mainly related to forest land (see previous section). Therefore, there is a trade-off between provision of fire protection and provision of most of the other ESS. This may explain why the provision potential values for this ecosystem service in the delimited multifunctional buffer zones are not much higher than in other areas. Thus, the selected areas maximize the provision potential of most ESS related to forest areas (Table 8) and these areas have a low provision potential for fire protection. Given the aforementioned results, it is proven that the consideration of the average ESS provision potential as a suitability layer ensured that the delimited multifunctional areas have a high potential for providing several ESS.
Landscape fragmentation also represents a challenge in relation to ensuring that a certain amount of land is delimited as multifunctional zones. This is because fragmentation leads to the selection of many small and disperse suitable patches that do not meet the criteria of minimum surface or adjacency to core areas and corridors. When these patches were excluded, the assigned area was thus significantly reduced. As seen in Table 9, the percentage of multifunctional areas relative to the total area of the region was much lower than the initially selected 20% of cells with the highest average provision potential value for the ESS considered. Table 9. Area of the GI zones delineated and the percentage of area that they cover relative to the total.

Zone
Area (

Conclusions
The proposed methodology successfully delimited multifunctional zones in GI with a high potential to provide ESS that contribute to counteracting the consequences of CC in areas with high biodiversity and natural values. The use of land cover as a proxy for mapping the ESS provision potential facilitates identification of the types of land cover with the highest potential to provide ecosystem services. Furthermore, combining several ESS provision potential maps in a suitability map to delimit GI not only ensures the GI multifunctionality, but also further aids identification of the types of land cover with the highest potential to provide several ESS.
However, application of the methodology in Galicia has shown that complex landscape patterns produce irregular GI areas and faces challenges in relation to delineating the minimum surface area. On the other hand, by identifying the types of land cover with a high potential to provide several ESS, planners can obtain an idea of how to fine-tune the results of the delimitation. Therefore, the boundaries of the delimited areas can be adjusted to include land cover types with high potential to provide ESS and exclude those with a lower potential, thus ensuring that the GI is close to areas where there is a demand for the ESS it produces.
Trade-offs between ecosystem services occur, and although an area may have a high potential to produce several ESS, it may also have a low potential to produce others. In this case, using land cover as a proxy may hamper the identification of areas with a high potential to provide several ESS. For example, although forests have a high potential to capture carbon, filter pollutants, and regulate hydrological cycle, they do not have high potential for fire prevention. Nevertheless, irregular and mixed landscapes of forests and fields may have a high potential to prevent fires while producing other ESS more closely linked to forests. In addition, landscape patterns can also create synergies or trade-offs among ESS and either increase or reduce their provision potential or increase the provision of ESS in areas close to zones of demand [28,36]. Qiu and Turner [37] suggest delimiting areas that are large enough to comprise spatial heterogeneity and reduce trade-offs. Similarily, Walz et al. [79] suggest that landscape units should be used to delimit GI multifunctional areas.
Bearing all of the above in mind, future methodologies should explore the potential of including landscape structure in ESS mapping, mainly in those cases where the ecosystem service is not clearly linked to a particular land cover type [36]. In this regard, tier III methods, as defined in Maes et al. [29] (which rely on a regression model to link data from a statistical survey with spatial variables), will be more efficient for analyzing the relationship between landscape pattern or landscape units and ecosystem service provision potential. However, these methodologies require input of large amounts of data.
Another possible way of improving methodologies for the spatial planning of GI is to delineate landscape units and identify those with the highest potential to provide several ESS. This would reduce the influence of landscape heterogeneity on the results of the GI zoning and thus produce more regular delimited areas. Funding: This research was funded by Fundación Biodiversidad of the Spanish Ministry of Ecological Transition, project title "INVERCLIMA-Infraestructura Verde para la adaptación de la ordenación territorial al cambio climático".

Conflicts of Interest:
The authors declare no conflict of interest.