Spatial Planning of Marine Protected Areas (MPAs) in the Southern Caspian Sea: Comparison of Multi-Criteria Evaluation (MCE) and Simulated Annealing Algorithm

Abstract

Protected areas are referred to around the world as the basis of conservation strategies. Designation of marine protected areas (MPAs) is to preserve marine biodiversity and protect species, habitats in the seas, and oceans. The simulated annealing algorithm (SAA) with other algorithms (swap iterative improvement, normal followed by two step, two step iterative improvement, and normal iterative improvement) in MARXAN conservation solutions software and the multi-criteria evaluation (MCE) method were used to locate MPAs in the Southern Caspian Sea. Then, four methods were examined for site selection that include: (1) Simulated annealing algorithm, (2) MCE with zonal land suitability (ZLS), (3) MCE with compactness and contiguity, and (4) combined method of multi-criteria evaluation with spatial constraints and a simulated annealing algorithm (improved MCE). In the MCE method, we applied different weighted scenarios to locate MPAs. The criteria for determining the desired regions of MPAs included 12 factors gathered in three groups, including: (1) Ecological criteria (distribution of fish Huso huso, Acipenser persicus, Acipenser stellatus, Rutilus frisii kutum, and Alosa braschnikowi; location of coastal protected areas, distance from coastal rivers (Coastline), distance from estuaries and deltas); (2) Physical criteria (distance from the coast, shore sensitive areas); and (3) Socio-economic criteria (distance from densely populated coastal cities, distance from industries near the coast). The results of comparing the algorithms in MARXAN 4.0.6 software showed that the simulated annealing algorithm has a better ratio of border-length/area than other algorithms. Also, the combined method of MCE (improved MCE) selects the best protection patches in terms of location, taking into account the seascape ecology metrics (e.g., patch compactness, edge density, normalized entropy, area metric for patches). Moreover, the results of the comparison of four methods for proposing MPAs based on seascape metrics showed that the combined method of MCE considers a protection network with more contiguity and compactness than the simulated annealing algorithm. The use of seascape ecology can help to preserve and create larger and denser patches in the arrangement of protective areas, because such a selection of protective areas is nature-inspired and can be more bold and appropriate in the course of conservation planning.

1. Introduction

Today, protected areas are referred to around the world as the basis of conservation strategies [1,2,3]. Designation of marine protected areas (MPAs) is to preserve marine biodiversity and protect species and habitats in the seas and oceans [4]. The term marine protected area (MPA) is defined in accordance with the International Union for Conservation of Nature (IUCN) as follows: “A clearly defined geographical space, recognised, dedicated and managed, through legal or other effective means, to achieve the long-term conservation of nature with associated ecosystem services and cultural values” (https://www.iucn.org/content/when-a-marine-protected-area-really-a-marine-protected-area) (accessed on 30 June 2022) [5].

The primary purposes of establishing such areas are the protection and rehabilitation of the marine environment, management of human activities and assessment of its effects on the marine environment, as well as sustainable exploitation of the world’s marine resources [6]. Currently, 8.16% of the world’s oceans are protected [7]. On the other hand, the EU Biodiversity Strategy for 2030 [8] requires the legal protection of a minimum of 30% of the EU’s marine area, of which 10% should be strictly protected [9]. In principle, the location of protected areas should be such that the site selected is representative of biodiversity [10,11,12,13], which in the past was usually done using expert opinions and some basic criteria. For this reason, many protected areas (due to incorrect location) have not been successful in conserving biodiversity [14].

Countries have the jurisdiction to create MPAs (like coral triangle as conglomeration in Southeast Asia) in the limits of their declared exclusive economic zones, which account for 39% of the world’s ocean area and include many hotspots of human activity and industry. Today, 6.3% of the world’s ocean within national jurisdictions is in implemented and fully/highly protected areas [15]. The top 10 fully/highly protected countries are Palau (78% or 477,418 km²), the United Kingdom (39% or 2,282,456 km²), Mauritius (29% or 637,916 km²), the United States (24% or 2,879,357 km²), Panama (21% or 68,771 km²), Ecuador (13% or 138,665 km²), Chile (12% or 452,140 km²), Kiribati (12% or 397,569 km²), Argentina (11% or 315,919 km²), and Australia (9.6% or 859,875 km²) [15]. In total, 75% of the global marine protected area is located in 42 very large (≥100,000 km²) implemented MPAs. The largest one is the Ross Sea region of 1,606,529 km² [15]. The Caspian Sea is the largest lake in the world, which is of great importance for five peripheral countries due to its special geographical location, habitat of unique species, fisheries-related industries (including fish resources), and offshore oil and gas resources [16,17,18]. Offshore oil and gas production, exploration, and transportation represent the main threat for the Caspian Sea marine environment due to serious oil pollution of seawater [19,20,21,22,23]. In the southern basin of the Caspian Sea, 24% of the fish are exclusive species of fish, 6% are highly endangered (sturgeon), and 21% need protection [24]. Population growth, agricultural development, oil and gas resources, maritime transport, the tourism industry, lack of an optimal wastewater treatment system, and the arrival of invasive species including Mnemiopsis leidyi (A. Agassiz, 1865) in the Southern Caspian led to water pollution, reduced seawater quality and biodiversity, destruction of vital habitats, and generally a decrease in the quality of life of local residents [25]. Therefore, in order to protect habitats and assess the possible effects of natural and human factors, the optimal determination and establishment of protected areas is essential, and can play an important role in preserving biodiversity, endangered species, and good practices of fisheries management on the southern shores of the Caspian Sea. Establishment of MPAs in the Caspian Sea is of high importance for sustainable development of the Caspian Sea region [26].

The Caspian Sea region has a big number of coastal natural state reserves and parks, but until recently, it had no MPAs. While there is no MPA in the southern Caspian Sea as of yet, it contains critically endangered species such as the Beluga. In February 2018, Azerbaijan announced plans for creation of the first MPA on the Caspian Sea by upgrading and expanding the Gizilaghaj State Reserve to become a national park which will include the first marine protected area. It was established in September 2018. The MPA will seek to protect six significant marine species on the brink of extinction, including the Beluga sturgeon and the Caspian salmon. The State Reserve located in south-eastern Azerbaijan was created in 1929 for the purpose of protecting migrant, swamp, and wild birds, and was declared a Wetland of International Importance by the Ramsar Convention on Wetlands in 2001. The reserve lies just south of the mouth of the Kura and Aras rivers, which provide the spawning grounds for different sturgeon species [27]. On 29 October 2021 the Secretariat of the Important Marine Mammal Areas (IMMA) announced that the Caspian seal—one of the red-listed species—in three areas of the Caspian Sea (Breeding Area, Moulting and Haul Out Areas, and Transitory Migration and Feeding Area) was officially awarded the IMMA status by the International Union for Conservation of Nature task force on marine mammal protected areas [28]. Establishment of these protected areas shows the progress in regional cooperation and a higher standard of marine environmental protection.

There are various approaches to propose spatial marine conservation planning, which can be considered as a multi-criteria evaluation (MCE) approach [29], including special software packages like MARXAN conservation solutions software [14,30,31]. The MCE approach includes techniques used to solve spatial planning problems [32] and simplify complex decision making that measures optimality using several criteria. MCE is a decision support tool in which the decision is a choice between several options (such as land allocation) and the basis of any decision is to know its criteria [33].

MARXAN software [30] is a conservation prioritization tool for generating protected area network design scenarios in marine conservation [34]. MARXAN’s goal is to solve the problem of selecting a set of proposed areas meeting quantitative conservation targets (e.g., X% of species Y or habitat Z) at minimal socio-economic cost (e.g., forgone opportunities). Many researchers have used the MCE method [35,36,37,38,39,40,41] and MARXAN software [9,30,42,43,44] to investigate the significance of MCE in the marine context.

Engelhard et al. [45] investigated MPAs in Moreton Bay in Eastern Australia. In the study, the indicators of landscape ecology indices (connectivity) were used for the relationship between marine habitats in marine protected areas. This study showed that the identification and prioritization of conservation spots without considering the indicators of landscape are not sufficiently represented in a reserve network.

Weeks [46] studied seascape connectivity to prioritize marine protected areas. In this study, MARXAN software algorithms were investigated for the relationship between marine habitats. This study showed that by using landscape indicators such as connectivity, it is possible to establish a relationship between breeding areas and birth areas of species and increase habitat stability.

Cajica et al. [47] investigated three indicators of seascape (habitat β-diversity, shape complexity, and connectivity) to prioritize conservation spots in Mexico’s Cozumel Reefs National Park (CRNP). The results of this study showed that seascape metrics can be used to prioritize protection spots and select the appropriate protection network.

Virtanen et al. [48] examined the spatial relationship of marine connectivity in spatial conservation. In this study, Zonation conservation planning software [49] developed in the Finnish Environment Institute (SYKE) was used to establish communication between marine protected areas. The results of this study showed that seascape indices, including the relationship between conservation spots, can contribute to the population dynamics and species of a habitat in marine protected areas and reduce the fragmentation of habitats in protected areas.

Javed et al. [44] used the satellite tracking data and the MARXAN software prioritization tool to propose the most important areas for flamingo protection in the coastal area of Abu Dhabi, and the results showed that this software shows the best solution with a more realistic area for proper protection.

Li et al. [50] used the MARXAN software to integrate inter-basin communication connections in the planning of lotic water protection in the South China Plain. The results showed that the existing protection system covers approximately 20% of freshwater wetlands.

Timonet and Abecasis [51] investigated an integrated approach to the analysis of ecological and socio-economic data for the MPAs of Portugal using the MARXAN software. The solutions provided by MARXAN, and subsequently improved with MinPatch, revealed that there is a need to considerably increase the area under protection. However, adequate protection could be achieved with a number of MPAs (between 5 and 14) similar to the number of already existing MPAs (n = 8).

According to the review of sources above, identification of marine protection areas has been done using a simulated annealing algorithm and multi-criteria evaluation for different ecosystems in the world, but so far, the effect of compactness and contiguity has not been studied in the multi-criteria evaluation method to select marine protection areas. Compactness means that cells with similar land-use are placed together. Meanwhile, contiguity means that patches with similar land-use are in the vicinity. Therefore, the innovation of the present study is to provide an improved multi-criteria evaluation model to determine suitable conservation areas for the Southern Caspian Sea. To propose marine protected areas in the southern basin of the Caspian Sea, MARXAN software algorithms, the zonal land suitability (ZLS) method, a multi-criteria evaluation method with spatial constraints (compactness and contiguity effect), and a combined method of multi-criteria evaluation with MARXAN software taking into account the seascape indicators—abbreviated as improved MCE—have been used.

2. Materials and Methods

2.1. Study Area

The Caspian Sea is the largest lake in the world with importance for five peripheral countries due to its special geographical location, habitat of unique species, fisheries-related industries (including fish resources), and offshore oil and gas resources. Its area is about 380,000 km², and it is politically, economically, and environmentally significant in the region and globally [16,17,19,52]. It has a maximum depth of 1025 m, average salinity of 10–13 PSU, and water temperature ranging between freezing point in the Northern Caspian and 28.6 °C in the Southern Caspian [52,53]. The study area is the southern basin of the Caspian Sea, which is located between 36.57° and 38.44° N, and 48.87° and 54.04° E (Figure 1). This area has an area of about 18,700 km² and is the main area for fishing whitefish in the coastal zone of Iran [54].

2.2. Methodology

Selection of criteria for MPA detection depends on the protection objectives, available data, available resources to analyze data, and collection of new data. The most important criteria for determining marine protected areas based on a review of natural resources for the study area in the present study include 12 factors gathered in 3 groups, including distance from the coast, distance from protected areas near the coast, distance from coastal rivers (coastline), distance from estuaries and deltas, distance from coastal sensitive areas, distance from densely populated coastal cities, distance from industries near the coast, habitats of Beluga (Huso huso, (Linnaeus, 1758)), Persian sturgeon (Acipenser persicus, (Borodin, 1897)), and Starry sturgeon (Acipenser stellatus, (Pallas, 1771)), Caspian Kutum (Rutilus frisii kutum, (Kamensky, 1901)), and Caspian Shad (Alosa braschnikowi, (Borodin, 1904)). The studied indicators, data, and data sources are shown in

Table 1. Criteria for proposing marine protected areas for the southern basin of the Caspian Sea.

Criteria for marine protected areas	Criteria Group and Final Weight	Criteria and Final Weight *	Sub-Criteria				Source
	Ecological (0.534)	Caspian Shad distribution (0.15)	Depth (all species); Chl a ** (Caspian Kutum, Caspian Shad); PAR (all species); SLA (Caspian Kutum, Caspian Shad); Water Temperature (all); Macrobenthos density, Bivalve density, Average organic matter, Bed material, Bed slope (sturgeon species)	Index	Sensor	Unit
				Water surface temperature	MODIS	°C	http://oceancolor.gsfc.nasa.gov (access on 18 December 2021)
				Chl a	MODIS	mg/m3	http://oceancolor.gsfc.nasa.gov
		Caspian Kutum distribution (0.15)
				Photosynthetic active radiation	MODIS	E/m2/d	http://oceancolor.gsfc.nasa.gov
		Beluga distribution (0.17)		SLA ***	Jason-1	cm	http://www.aviso.altimetry.fr (access on 18 December 2021)
				Depth	_	m	Based on the hydrographic data of the surveying organization of Iran
		Starry sturgeon distribution (0.17)		Macrobenthos density	-	n/m2	[55]
				Bivalve density	-	n/m2
		Persian sturgeon distribution (0.17)
				Average organic matter	-	℅
				Bed material	-	-
				Bed slope	-	℅
		Protection of areas with high ecological value (0.09)					Wetlands and protected areas on the Caspian Sea coast (Gomishan Wetland, Miankaleh Wetland, Bojagh National Park, Anzali Wetland, Lundville Protected Area)
		Distance from coastal rivers (0.05)					[56]
		Distance from the estuary, delta (0.05)					[56]
	Physical (0.233)	Distance from the beach (0.4)					Google Earth
		Sensitive coastal areas (0.6)	Coastal geomorphology	Sloping cliffs and high walls (1) Medium cliffs with serrated margins (2) Short cliffs, alluvial plain (3) Sandy beach, estuaries and coastal wetlands (4) Sandy and river deltas (5)			[56]
			Height (inland part of the beach) and depth	0 to 20% (1) 21 to 40% (2) 41 to 60% (3) 61 to 80% (4) 81 to 100% (5)			Prepared from hydrographic, elevation and elevation and depth points of the country’s surveying organization
			Layer of coastal structures (ports)	-			Google Earth
			Water level (water level)	0 to 20% (1) 21 to 40% (2) 41 to 60% (3) 61 to 80% (4) 81 to 100% (5)			http://www.aviso.altimetry.fr (access on 18 December 2021)
			Continental shelf	-			Based on depth line −50
			Population density of coastal cities	0 to 20% (1) 21 to 40% (2) 41 to 60% (3) 61 to 80% (4) 81 to 100% (5)			Statistical Calendar of Iran (2012)
	Socio-economic (0.233)	Distribution of residential areas (0.5)					Google Earth, Statistical Calendar of Iran (2012)
		Distribution of human activities (0.5)					WebGIS data

* Weights defined using the analytic hierarchy process method for every factor. ** Chlorophyll a is a measure of the amount of algae growing in a water body. *** A SLA—sea-level anomaly derived from satellite altimetry data.

.

Based on Table 1, two general phases of the analysis were performed in this research: (1) Species distribution modelling, and (2) Establishing the position and boundaries of MPAs.

2.2.1. Species Distribution Modelling Using MCE

Three species of sturgeon (Persian sturgeon, Starry sturgeon, and Beluga) represent endangered species, and Caspian kutum and Caspian shad were selected as commercial species regarding to the ecological criteria presented in Table 1. First, the database (ecological criteria of Table 1) is formed, and then the determined criteria are used for determining the suitability level of the criteria using fuzzy membership functions. Then, different weighting scenarios are determined based on the analytical hierarchical process. Finally, we model species distribution using MCE.

In the MCE model, in order to standardize (normalize) the decision rule criteria using fuzzy membership functions, the weighting of the criteria was done based on the analytical hierarchical process. Then, the accuracy of the model was assessed based on the descriptions below. Fuzzification is used to standardize variables. Since different variables have different units, fuzzification is used to make the evaluation scale the same for them. The equations used for defining fuzzy membership functions are shown in

Table 2. Fuzzy relationships of the decision rule criteria.

Type of Function	Shaped Fuzzy Functions	Equations Calculating the Numerical Values of Functions Per Pixel *
Linear function (L)		$X_i=\left(1-\frac{(R_i-c) }{\left(d-c\right)}\right)\times standardized \_range$	[57]
		$X_i=\left(\frac{R_i-a }{b-a}\right)\times standardized \_range$
	Rate each pixel before applying fuzzy and standardized range Ri $\mu =\cos \alpha \times 255$
Sigmoid function (S) **		$\alpha =\frac{\left(x-c\right)}{\left(d-c\right)}\times \frac{\pi }{2}$	[57]
		$\alpha =\left(\frac{\left(x-c\right)}{\left(d-c\right)}\right)\times \frac{\pi }{2}$
		c is first control point, d is second control point, pi = 3.14159
J-shaped function		$\mu =\frac{1}{1+\frac{{\left(x-c\right)}^2}{\left(d-c\right)}}\times 255$	[57]
		$\mu =\frac{1}{1+\frac{{\left(x-b\right)}^2}{\left(b-a\right)}}\times 255$

* a, b, c, and d are the control functions of fuzzy functions. In fact, they determine the threshold limits in standardization. ** Increasing sigmoid function for all criteria of all species except depth criteria for Caspian kutum and Caspian shad with decreasing sigmoid function.

.

The analytic hierarchy process, presented by Saaty [58] as one of the most efficient decision-making techniques, was used in weighting the criteria. This method is based on pair-wise comparisons. In the analytic hierarchy, matrices with CR values higher than 0.1 must be reevaluated, and values lower than 0.1 must be included in the process of determining the weight of the criteria [58]. To get more accurate criteria weights, different weighted scales were created to observe and measure the effect of each criterion on the modeling results. In fact, these five scenarios (

Table 3. Weight scenarios modeling to achieve maximum accuracy.

(A): Distribution of species of three sturgeon species
Scenario	Macrobenthos Density	PAR	Bivalve Density	Depth	SST	Average Organic Matter (TOM)	Bed Material	Bed Slope	CR
1	0.125	0.125	0.125	0.125	0.125	0.125	0.125	0.125	0
2	0.0728	0.1013	0.1239	0.0413	0.0413	0.3717	0.1239	0.1239	0
3	0.0831	0.1815	0.1515	0.0999	0.0605	0.1815	0.1815	0.0605	0
4	0.2111	0.0704	0.0704	0.1725	0.2111	0.0704	0.124	0.0704	0
5	0.0704	0.2111	0.1725	0.0704	0.0704	0.124	0.0704	0.2111	0
(B): Distribution of Caspian kutum and Caspian shad
Scenario	Depth		CLa		PAR	SLA	SST		CR
1	0.2		0.2		0.2	0.2	0.2		0
2	0.5201		0.2292		0.856	0.856	0.0797		0.02
3	0.1304		0.3913		0.3913	0.435	0.0435		0.04
4	0.1304		0.0435		0.0435	0.3913	0.3913		0.01
5	0.1008		0.2145		0.1217	0.2815	0.2815		0.04

) indicate the importance of each criterion in the final produced map and the extent of their effect on the spatial distribution of the species. The modeling of the weighted scenarios is shown in Table 3. These weights are based on the scenarios of the ordered weighted average (OWA) evaluation method, in which different weights include the same weights to all variables, or increase the weighting for some of them in different scenarios [40].

MCE is one of the most common methods for multi-criteria spatial evaluation. In this method, the weight of the criteria and the constraint on the criteria are applied. First, all the criteria are multiplied by the weight of the factor (compensation weight). (Equation (1)) [57]:

$$S={\displaystyle \sum }{W}_i\times X_i$$

(1)

where S indicates the suitability rate, W_i represents the weight of factor I, and X_i is the criterion score in factor i.

An accuracy assessment of MCE was made based on the receiver operating characteristic curve (ROC) index. The area under the ROC curve is used to compare a continuous map with a Boolean image (an image at the range of 0 and 1 which is related to the ground truth map) based on calculating false positive and true positive rates. A value of 1 represents a full spatial agreement, and a value of 0.5 shows a random matching agreement [57].

The ROC (Equation (2)) is extracted from the contingency tables based on comparing the simulated map with the reference map. In the present study, the reference map was the average of fishing data from the fisheries organization of Iran for Caspian kutum and Caspian shad, as well as from the Caspian Environment Project for three sturgeon species in waters of Mazandaran Province.

$$AreaUnderCurve={\displaystyle \sum _{i=1}^n[x_{i+1}}-X_i]\times [y_i+\left(y_{i+1}-y_i\right)/2]$$

(2)

After species distribution modelling, overlap of the species’ ecological niche was investigated. Ecological niche overlap between species x and y is calculated as follows:

$$Range\; overlap\; index =\left(N_{x,y}/\min \left[N_x,N_y\right]\right)$$

(3)

where N_x,y represent the number of pixels in which both species exist based on the model prediction, and N_x,y represent the predicted pixels of species x and species y, respectively.

In this method, Schoner D (Equation (4)) and Statistical Index I (Equation (5)) are two indices that are related to each other, and the closer they are to 1, the more niches overlap of the two species [59].

$$D\left(p_x,p_y\right)=1-\frac{1}{2} {{\displaystyle \sum }}_i\left|p_{x,i}-p_{y,i}\right|$$

(4)

$$I\left(p_x,p_y\right)=1-\frac{1}{2} \sqrt{{{\displaystyle \sum }}_i{\left(p_{x,i}-p_{y,i}\right)}^2}$$

(5)

2.2.2. Establishing the Position and Boundaries of MPAs

Criteria for establishing the position and boundaries of MPAs include 12 factors gathered in 3 groups (Table 1) including: (1) Ecological criteria (distribution of Huso huso, Acipenser persicus, Acipenser stellatus, Rutilus frisii kutum, and Alosa braschnikowi, distance from protected areas near the coast, distance from coastal rivers (coastline), distance from estuaries and deltas); (2) Physical criteria (distance from the coast, coastal sensitive areas); and (3) Socio-economic criteria (distance from densely populated coastal cities, distance from industries near the coast). Species distribution modelling was done according to methodology described in Section 2.2.1. Criteria related to distance parameters were calculated by using the distance function using ArcGIS 10.3 software, which calculates the shortest distance for each factor.

For coastal sensitive areas, we used the coastal vulnerability index (CVI). The CVI is the most widely used index for determining the sensitivity of beaches [60]. Both variables and stresses are involved in calculating this index. It is used to determine the relative and cumulative sensitivity of the coastal zone. To calculate this index, 6 indicators (Table 1) are used, and their ranges or scores are ranked on a scale from 1 to 5, where 5 indicates the highest sensitivity and 1 indicates the lowest sensitivity of that indicator (a_i). Then, the geometric mean of 6 indicators was calculated based on Equation (6).

$$CVI={\left(a_1\times a_2\times a_3\times a_4\times a_5\times a_6\right)}^{1/6}$$

(6)

After calculation of all effective criteria, establishing the position and boundaries of MPAs were performed based on two approaches—simulated annealing algorithm (SAA) [61,62] (in MARXAN software) and multi-criteria evaluation (MCE). MARXAN software has different types of artificial intelligence algorithms, which include swap iterative improvement, normal followed by two step, two step iterative improvement, and normal iterative improvement and simulated annealing algorithm (SAA). The capabilities of this software include the selection of protection units, production of conservation maps in raster and vector formats, and the evaluation of final maps based on the ratio of boundary length to area (P/A).

In the present study, 3 methods were used for the MPA site selection:

(1)

MCE with zonal land suitability. In this method, the goal is only to select areas with high suitability (Equation (7)) [63].

$$\mathrm{Sz}=\frac{{\displaystyle \sum }\left(\mathrm{Li}\right)\mathrm{z} }{\mathrm{nz}}$$

(7)

where Sz = zonal land suitability, (Li)z = local suitability of the cells I belonging to the zone z, nz = number of cells forming z zones. The minimum area for MPA is set as 10,000 hectares.

(2)

MCE with spatial constraints (compactness and contiguity). In this method, in addition to the high degree of stability, high contiguity and compactness of patches (seascape indices) are also considered. This criterion is expressed as Equation (8) (standardized equation of Bribiesca [64]):

$$C_{SC}={\displaystyle \sum }\left(\frac{CP-C_{min}}{C_{max}-C_{min}}\right)$$

(8)

where C_SC is standardized compactness, CP is ratio of P/A, C_min is minimum compactness—or in other words maximum dispersion, and C_max is maximum compactness or minimum dispersion. The number of patches in a seascape is an index for the intensity of contiguity. Accordingly, in this study, the calculation of the number of patches and its standardization according to Equation (9) has been used as a contiguity criterion.

$$C_G={\displaystyle \sum }\left(\frac{C_{G{T}_{max}}-C_{GT} }{C_{G{T}_{max}}-C_{G{T}_{min}}}\right)$$

(9)

where C_GT represents the sum of number of land use patches. To standardize this criterion, the minimum number of patches C_Gtmin was considered to be 1, and the maximum number of patches C_Gtmax was equal to the total number of target pixels in each land use [56].

(3)

Improved MCE, which is a combination of MCE with spatial constraints and SAA.

Finally, the last step of the present study was to compare the output maps of the MPA position and boundaries based on the following seascape metrics [65]:

$$\mathrm{E}=-({{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{k}}\mathrm{Pi}+\ln \left(\mathrm{Pi}\right))/\ln \left(\mathrm{n}\right)$$

(10)

where E is the normalized entropy (frequency of each patch (heterogeneous) among the total patches of the seascape); n is the number of image classes; Pi is the ratio of class I to total classes in image; I is the class index in the neighborhood unit; and k is the total number of classes in the neighborhood.

$$\mathrm{Ed}=\frac{\mathrm{E}}{\mathrm{A}} \left(10000\right)$$

(11)

where Ed is edge density (fragmentation index). Ed equals the sum of the lengths (m) of all edge segments in the seascape, divided by the total seascape area A (m²), multiplied by 10,000 (to convert to hectares).

$$\mathrm{C}=\mathrm{SQRT}\left(\mathrm{Ap}/\mathrm{Ac}\right)$$

(12)

where C is the patch compactness (patch shape relative to a complete circle); Ap is the patch area under calculation; and Ac is the area of a circle that has an area equal to the patch being calculated. In addition to three above mentioned metrics, area metric for patches (ha) was also calculated.

All the above processes are summarized in the Figure 2. In this paper, data pre-processing and preparation were done in ArcGIS 10.3 software. The MCE method and simulated annealing algorithm were done in Terrset 18.31 and MARXAN 4.0.6 software, respectively.

3. Results

3.1. Species Distribution Modelling Using MCE

MCE modelling of the species distribution allowed to construct maps of the habitat suitability for different species in the southern basin of the Caspian Sea, which has the highest probability (Figure 3).

3.2. Establishing the Position and Boundaries of MPAs

3.2.1. Comparison of the Algorithms Available in MARXAN Software and Selection of the Best Algorithm

Figure 4 shows the results of the Scenario 1 analysis where we investigated a relationship between border length and area of selected MPAs at different boundary length modifiers (BLM) that determine the optimum ratio of boundary length to area (P/A) of MPA. In Scenario 2, we investigated different conservation targets, of 40%, 50%, and 60% (target means X% of species Y or habitat Z for conservation), and compared the appropriate maps of optimal location and configuration of the resulting MPAs (Figure 5).

In the third scenario, we compared different algorithms with a BLM = 60 and a target of 30%, based on the diagram of border length vs. area of MPAs (Figure 6).

Table 4. Spatial characteristics of the MPAs calculated with different algorithms.

Algorithm	Area (ha)	Border Length (km)	P/A
Swap iterative improvement	299,191	837	0.27
Normal followed by two steps	291,906	802	0.27
Two steps iterative improvement	302,897	794	0.26
Normal iterative improvement	396,316	913	0.23
Simulated annealing algorithm	289,813	822	0.28

shows the spatial characteristics of the MPAs calculated with different algorithms. Figure 7 shows the spatial location and configuration of the MPAs calculated by different algorithms. For example, swap iterative improvement has a worse result, as fragmentation in MPAs is high in comparison with other algorithms. On the other hand, the best algorithm is simulated annealing, because the fragmentation of the resulted MPAs is the lowest.

3.2.2. Multi-Criteria Evaluation (MCE) Method

The GIS layers built from 12 parameters were calculated with appropriate weights indicated in Table 1. The highest weight was given to ecological criteria, and among the ecological criteria, the habitats of endangered species were set up as more important. To reduce the weighting error and the appropriate effect of the used parameters, multi-criteria evaluation was performed for each group of parameters. As a result, a map of optimal location and configuration of the MPAs along the coast of Iran resulting from the MCE method based on a 0–255 fuzzy range was plotted (Figure 8).

3.3. Site Selection for MPAs

Figure 9 shows the results of site selection for potential MPAs obtained with the following three different methods: (a) the zonal land suitability method, (b) MCE with compactness and contiguity, and (c) combined method of the improved MCE.

3.4. Comparison of Different Methods Based on Seascape Metrics and Selection of the Best Method

To compare different methods for the MPAs site selection, four seascape metrics were calculated and compared: (a) normalized entropy, (b) edge density, (c) area (ha), and (d) compactness (Figure 10).

4. Discussion

Based on results shown in Figure 3, the highest suitability region for Caspian kutum is located in the middle and western parts of the coastal zone of Iran, which is in agreement with [54,66]. The spatial distribution of Caspian Shad in the central and eastern part of the coastal zone corresponds to the results obtained by [13]. The most preferred areas for Persian sturgeon, Starry sturgeon, and Beluga are located in the eastern, western, and central parts of the coastal zone, respectively, which is in agreement with [67,68]. We found that there are high overlaps of environmental niches based on the obtained statistics among all three species of sturgeon, which may be related to the early period of their life and the fact that they have similar nutrition [69,70]. The calculated Coastal Vulnerability Index has shown that the eastern shores of Iran—including Gorgan Bay, Miankaleh Shelter, and the Gomishan Wetland—are the most vulnerable, which is consistent with studies by [71,72].

Figure 6 and Figure 7 and Table 4 suggest that the more compact patches are more important for conservation from an ecological and management point of view. Accordingly, the simulated annealing algorithm has better compactness (low fragmentation) and less area than the other algorithms. This is consistent with the results obtained by [9,42,43,44,61,73,74]. The simulated annealing algorithm also has suitable contiguity among other methods.

The results showed that the zonal land suitability method (Figure 9a) focuses on the preference of the region for conservation and does not pay attention to the spatial distribution of protective patches, which is consistent with [38,39,40,41,63,75,76]. In the multi-criteria evaluation method with spatial constraints, in addition to the patch area, the factors of compactness and contiguity (Figure 9b) are also taken into account. In the improved MCE method (Figure 9c), a cost layer (cost map was prepared from the reverse of the prioritization map of simulated annealing algorithm) was added, which led to a higher integration of several MPAs selected.

The amount of entropy increases with increasing numbers of heterogeneous patches in the seascape [77]. The results of our study showed that the highest value of this index is related to the zonal land suitability method, and the lowest is related to the combined method of the Improved MCE (Figure 10a). Increasing edge density means increasing fragmentation of the MPAs, which means reducing the biodiversity. According to Figure 10b, the highest edge density of 249.3 m/ha is related to the zonal land suitability method, and the lowest density of 85.09 m/ha is related to the combined method of the improved MCE. Calculation of the MPAs “area” index with different methods (Figure 10c) showed that the highest area is related to the simulated annealing algorithm, and the lowest patch area is related to the multi-criteria evaluation method with spatial constraints. It should also be noted that the MPA area index alone is not a good criterion for choosing the appropriate method, because the increase in the area sometimes is accompanied by an increase in length of the MPAs border, which leads to separation of areas and discontinuity of the protection network. This issue can be seen when using the simulated annealing algorithm (Figure 10c). Therefore, the compactness index should be investigated [9,43,74,78]. Based on the obtained results (Figure 10d), the lowest compactness index is related to the zonal land suitability method, and the highest value of this index is related to the combined method of the improved MCE.

Overall, landscape ecology (seascape in our study) is the study of the interaction of spatial patterns with ecological processes that has led to the development of new models and theories of spatial relationships (a new set of spatial pattern data and their dynamics) as well methods of investigation, which is rarely considered in ecological discussions. The results of comparing the above mentioned indices and their effect on the selection of the MPAs have shown that areas with low margin density, higher compactness, and lower entropy have a higher value in terms of protection, which is in agreement with studies [45,46,47,48]. In the present study, the combined method of the improved MCE is superior to other studied methods in terms of its ability to select a more compact protection network.

5. Conclusions

There are various approaches to propose spatial marine conservation planning and organization of MPAs, which can be considered as a multi-criteria evaluation approach. In the present research, identification of marine protection areas was done using the simulated annealing algorithm and the multi-criteria evaluation methods. For the first time, the effect of compactness and contiguity was studied and taken into account in the multi-criteria evaluation method to select marine protection areas. So, the innovation of the present study was to prepare an improved multi-criteria evaluation model to determine suitable conservation areas for the Iranian coastal zone of the Southern Caspian Sea. To propose site selection of marine protected areas, MARXAN software algorithms, the zonal land suitability method, the multi-criteria evaluation method with spatial constraints (compactness and contiguity effect), and the combined method of multi-criteria evaluation with MARXAN software taking into account the seascape indicators—abbreviated as improved MCE—were used. The results of comparing the algorithms in MARXAN software showed that the simulated annealing algorithm had a better ratio of border length/area than other algorithms. Also, the improved MCE selects the best protection areas in terms of location, taking into account the seascape ecology metrics (e.g., patch compactness, edge density, normalized entropy, area metric for patches). Moreover, the results of the comparison of four methods to select a location for MPAs based on seascape metrics showed that the combined method of MCE has revealed a protection network with more contiguity and compactness than the simulated annealing algorithm.

The improved MCE method, which is a subset of multi-criteria evaluation methods, can facilitate management of marine protected areas by considering landscape indices. Choosing areas with more compactness leads to easier identification of the codes of protected areas, and on the other hand, choosing fewer large protective areas is better than choosing more spots with no habitat connection, which leads to be more feasible and practical in the management of protected areas. Moreover, protection grading of the areas is in the next stage, in which zones with high sensitivity (no fishing in these areas) and zones with recreational fishing are also determined. According to the application of compactness and contiguity indicators, the above method can propose MPA sites that do not meet all the criteria. Also, this method reduces the effect of some criteria by applying some landscape indices. On the other hand, it will increase the habitat connectivity of the species and also facilitate the management of marine protected areas by the government.

In this regard, it is suggested that the developed methodology can be used to locate MPAs in further studies. It is also recommended to consider more criteria such as Caspian seal (Pusa caspica, (Gmelin, 1788)) habitat suitability, etc. The obtained results in determining the appropriate MPAs in the southern basin of the Caspian Sea are recommended for decision makers; for example, regarding the best location and configuration of the MPAs along the Iranian coast of the Caspian Sea, which is shown in Figure 9c.

MPAs play an important role in marine ecosystem services, and a healthy ecosystem has positive impacts on human welfare. According to the Millennium Ecosystem Assessment [79], ecosystem services are grouped into four categories: provisioning, regulating, supporting, and cultural. MPAs contribute to all four categories of ecosystem services in a variety of ways. The “provisioning services” may include marine products, biochemical resources, and genetic resources. The “regulating services” may include regulation of biodiversity, habitats of different species, fishery services, etc. The “supporting services” of coastal ecosystems include nutrient cycling, biologically mediated habitats, and primary production. The “cultural services” of coastal ecosystems include inspirational aspects, recreation and tourism, marine science, and education. Understanding of ecosystem services requires a strong foundation in the behavior of the marine environment, ecosystems, and ecology, which describes the underlying principles, processes, and interactions of marine organisms and the environment. Identification and establishment of the MPAs can be regarded as a useful tool for scientists to understand how relationships are interwoven among marine organisms and the marine environment via a variety of processes. The issues related to marine ecosystem services development in the Caspian Sea Region are not discussed in all five riparian countries. Only recently, several papers appeared dealing with marine aquaculture [80] and spatio-temporal variability of wind energy potential in the Caspian Sea [81]; another research project aimed to investigate Caspian seals [82]. We hope that the study presented in this paper will stimulate development of ecosystem services related to the establishment of marine protection areas in the Iranian coastal zone of the Caspian Sea. Unfortunately, it is too early to discuss the MPA-related ecosystem services in the whole Caspian Sea, because today, there is only one MPA in the Caspian. Today, the primary task is the scientifically-based selection of marine protected areas in the Caspian Sea and their legal establishment in the areas of responsibility of all Caspian countries. In this regard, according to the methodology, it is suggested that in other studies, the methodology of the present research can be used to locate marine protected areas and also consider more criteria such as Caspian seal (Pusa caspica) habitat suitability. In addition, the use of the output of the current results in determining the appropriate marine protected areas in the southern basin of the Caspian Sea is recommended for decision makers; for example, regarding the best location and configuration of the MPAs along the Iranian coast of the Caspian Sea.