Comprehensive decision-strategy space exploration for efficient territorial planning strategies

Multi-Criteria Decision Analysis (MCDA) is a well-known decision support tool that can be used in a wide variety of contexts. It is particularly useful for territorial planning in situations where several actors with different, and sometimes contradictory, point of views have to take a decision regarding land use development. While the impact of the weights used to represent the relative importance of criteria has been widely studied in the recent literature, the impact of order weights determination have rarely been investigated. This paper presents a spatial sensitivity analysis to assess the impact of order weights determination in Multi-Criteria Analysis by Ordered Weighted Averaging. We propose a methodology based on an efficient exploration of the decision-strategy space defined by the level of risk and trade-off in the decision process. We illustrate our approach with a land use planning process in the South of France. The objective is to find suitable areas for urban development while preserving green areas and their associated ecosystem services. The ecosystem service approach has indeed the potential to widen the scope of traditional landscape-ecological planning by including ecosystem-based benefits, including social and economic benefits, green infrastructures and biophysical parameters in urban and territorial planning. We show that in this particular case the decision-strategy space can be divided into four clusters. Each of them is associated with a map summarizing the average spatial suitability distribution used to identify potential areas for urban development. We finally demonstrate the pertinence of a spatial variance within-cluster analysis to disentangle the relationship between risk and trade-off values.


INTRODUCTION
Multi-criteria decision analysis (MCDA) is a wellknown decision support method that allows for the evaluation of decision alternatives based on a given set of criteria. MCDA serves as a general framework for supporting complex decision-making situations with multiple and often conflicting objectives that stakeholders groups and/or decision-makers value differently [1]. By its generic nature, MCDA is an "umbrella term to describe a collection of formal approaches which seek to take explicit account of multiple criteria in helping individuals or groups explore decisions that matter" [2]. MCDA can be used in a wide variety of contexts to deal with complex issues involving conflicting interests. In particular, it can be easily applied in GIS environment to analyze and produce land-use suitability maps [3]. Examples of applications include the mapping of wilderness [4], ecosystem services assessment [5], disease susceptibility mapping [6,7], or soil fertility evaluation [8] to name but a few. Within this context, each criterion is represented by a map and the method is applied on a pixel basis to evaluate its suitability in the light of the specific problem raised.
As it is the case in any modeling framework, MCDA's outputs are subject to uncertainty generated at different levels of the process. Indeed, it exists several sources of uncertainty regarding MCDA's model, inputs and parameters that might impact the final suitability map on which the decision will be based. If we consider MCDA using Ordered Weighted Averag- * Contributed equally to this work. † Corresponding author maxime.lenormand@irstea.fr.
ing (OWA) in a GIS context, four main sources of uncertainty can be identified in the aggregation process: the data used to build the criteria, the method used to generate the geographical layers (scale, standardization...), the methods used to determine the criteria weights and the order weights. While uncertainty of criteria weights has been already widely investigated in the recent literature [7,9], the sensitivity of MCDA outputs to order weight selection is a topic that is often neglected. The rigorous selection of order weight in MCDA is of major importance for discussing the outputs with policy makers [4]. Most of the studies focus on a few typical sets of order weights without clear method to select them or justify the choice of these weights. This is probably due to the lack of knowledge and process to build a rigorous experimental design to properly conduct a sensitivity analysis.
In this work, we tackle this problem by proposing an approach built around a method that allows for an automatic generation of order weights according to a given level of risk and trade-offs [10]. The main advantage of this method lies in the formalization of the relationship between risk (likelihood that the decision made for a given location is wrong) and trade-off (degree of compensation between criteria) providing a clear definition of decision-strategy space. Thus, a full exploration of the decision-strategy space becomes possible, fulfilling a major gap regarding the sensitivity of MCDA to order weights selection. With an experimental design defined properly followed by a clustering analysis, the decision-strategy space can be segmented into several clusters of order weights leading to the production of similar suitability maps. The resulting clusters are further analyzed to extract an average suitability maps summarizing the information contains in each cluster. We can finally perform a spatial variance within-cluster analysis to measure the degree of uncertainty associated with the average suitability maps. In what follows, we will describe in detail the methodology before presenting the results highlighting the importance of order weights selection in practical applications.
We illustrate, here, our approach through an ecosystem services (ES) valuation for an urban development perspective in South of France, particularly concerned by the growing demands for urbanization. Allowing to evaluate a high diversity of benefits and impacts, the concept of ES has gained importance for urban and coastal planning [11,12] and successfully combined with an MCDA framework [5] to undertake urban governance and related policies. Particularly, a promising use of ES as a decision-tool seems to reveal and anticipate trade-offs between the different ES and urbanization [5]. This is why we selected this concrete application to demonstrate the capacity of our approach for identifying potential spatial decisions regarding urban development in the light of different ES trade-offs.
The study area and the data are first introduced in the next section. We then describe in detail the methodology before presenting the results highlighting the importance of order weights selection in practical applications. The results and the potential limitations of our approach are finally discussed in the last section.

STUDY AREA, DATA AND CRITERIA
The territory of Thau, situated on the French Mediterranean Coast about 20km West of Montpellier, is an important populated area of 60,000 ha and 117,000 inhabitants. With a central lagoon of 7500 ha and 35 km of coastlines, it is a land dominated by water. The territory is remarkably rich in terms of biodiversity, natural areas and landscapes, from the coastal and agricultural plains to the wetlands, the wooden reliefs and the garrigue. The region is highly attractive and subject to a heavy demographic pressure. To conciliate urban development, traditional economic activities such as wine production or shellfish farming and protection of the ecosystems is the main challenge of the territory. To do so, several instruments of water management and territorial planning exist and an engineering structure, the SMBT (Thau Bassin Joint Association, Syndicat Mixte du Bassin de Thau in French), has been established. This contributed to build up a shared political vision and a common governance over the whole territory [13]. Our main challenge to support territorial planning needs in the region was to evaluate land suitability in order to identify the most suitable areas for spatial urban development. In agreement with the SMBT, which was our main partner for interacting with local stakeholders, we favored an approach based on the assess- We evaluated a high number of ecosystem services based on local expert knowledge. We followed the methodological framework of the capacity matrix, developed by [14] and improved by [15] and [16], which aims at assessing different land use and land cover (LU/LC) types capacities to provide ecosystem services. In practice, we applied the methodological framework presented by [17], using available LU/LC data ( Figure 1) and related biogeophysical information (see the Appendix for more details). In close collaboration with the SMBT, we determined a set of both, relevant services based on the Common International Classification of Ecosystem Services [18] and LU/LC types drawn on an accurate land cover map over the territory [19]. Here we leverage expertise from a participatory workshop with local experts and contacted individuals supported by the SMBT to complete the matrix [16]. In order to assess the actual supply of ES for the study region, we applied the matrix method, as aforementioned, a qualitative valuation approach that links spatially explicit LU/LC types to a defined set of ES. Each expert proposes a score to assess the capacity of each LU/LC type ( Figure 1) to provide to each ES a value of importance. The capacity matrix method reduces the complexity of human-environmental systems, allowing reproducibility. It has proven useful in addressing the urgency-uncertainty dilemma for ES assessments, i.e., to provide best available knowledge where detailed modeling is not feasible or where data gaps obstruct more explicit approaches [15].
Nine ES were considered at the end, on the bases of the importance provided by the local experts and stakeholders from the Thau region. These ES are listed in Table 1. Two provisioning services focus on food from crop production and gathering. The four regulating services are related to main territorial challenges in terms of conservation of natural habitats for biodiversity (i.e. nesting sites, breeding grounds,...), including maintaining a good biochemical water status, protection against erosion and flooding. Cultural services also ranked high, reflecting the importance that people provide to the benefits from the surrounding nature, through contemplating beautiful landscape or open air activities. Finally, risks were also considered, in particular, the risk of wildfires that is pondered a main disservice of Mediterranean ecosystems.
Each ES were spatialized using the high resolution (5×5m 2 ) LU/LC map ( Figure 1) and the average score given to the different LU/LC types by the experts. In order to refine the services assessment that is essentially based on the LU/LC type, we added external biophysical data to adjust the evaluation. We added for instance information about soil quality to the cultivated crops ES or a fire hazard index using physical informations (slope, exposure, wind, etc) to the disservice of wildfires. More details about this process are available in the Appendix. Note that we also added a tenth layer to measure the importance of habitat patches for landscape connectivity [20] using the GUI-DOS software [21].
It is worth noting that the ten standardized geographical layers described above have been normalized to obtain a suitability for urban development inversely proportional to the score given by the expert. At the end of the process, for a given criteria, the pixel value are ranging from 0 when it is unsuitable for urban land use to 1 in the opposite case. All the criterion maps are available in Figure S1. All the information was stored in a matrix Z which value z ij represents the suitability of pixel i for urban land use according to criteria j.
Finally, the relative importance of each criteria (i.e. criterion weights) is the number of experts that considered the ES important for the territory of Thau divided by the total number of experts (Table 1). An exception for the last criterion (the connectivity) that has been set to one after discussion with experts of the SMBT.

DECISION-STRATEGY SPACE'S EXPLORATION
This section describes the methods used to efficiently explore the decision-strategy space in MCDA using OWA operators. This methodology is summarized in Figure 2.

OWA operator
In this work, we rely on an Ordered Weighted Averaging (OWA) multi-criteria operator [22] to combine all the criteria together in order to obtain a spatial explicit representation, that was a final map, to allow decision makers and local expertes to base planning decisions. OWA operators can indeed be used in a GIS context [23] in which every pixel is considered independently and defined by a set of criteria values.
Two types of weights can be considered while using OWA operators: the criterion weights representing the relative importance of the criteria in the decision process, and the order weights characterizing the level of risk and trade-off taken in the decision. Criterion weights are represented by a vector V = (v j ) 1≤j≤n , where n stands for the number of criteria (n = 10 in our case). The j th criterion weight corresponds to the relative importance given to the j th criteria. In this particular application, the criterion weights are given by a panel of experts (Table 1) and reflects the importance of the criteria regarding urban land use. Order weights are represented by a vector W = (w j ) 1≤j≤n . These weights are applied to the ranked criteria on a pixel basis. For a given pixel, the first order weight w 1 is assigned to the criterion with the lowest value, the second order weight w 2 to the second lowest criterion... While criterion weights are usually based on expert knowledge, order weights are used to adjust the level of risk and trade-off in the aggregation process. Note that in both cases, the vectors of weights sums to 1, ∑ n j=1 v j = ∑ n j=1 w j = 1. Formally, the OWA operator is applied to every pixel i with the following formula.
where z i(j) is the j th lowest element of the collection of criteria z ij for the pixel i and v (j) is the j th criterion weight reordered accordingly.

Risk and trade-off
As mentioned above, order weights are deeply linked with the concepts of risk and trade-off in the decision-strategy process. Different order weights combinations (w 1 , ..., w n ) ∈ [0, 1] n lead to different levels of risk and trade-off values, and we will see later that the reciprocal is not so clear.
For a given pixel i, a combination of order weights that favors low z ij values, W = (1, 0, ..., 0) for example, represents a risk-aversion position. On the contrary, a combination of order weights that gives full weight to high z ij values, W = (0, 0, ..., 1) for example, represents a more risk-taking attitude. In the latter situation the risk is to give a high suitability to a pixel only based on the criteria with the highest suitability without taking into account the information given by the other criteria. As the name suggests, the trade-off represents the degree of accommodation between criteria. It can be seen as a measure of dispersion over the order weights. The maximum level of trade-off is reached when w j = 1 n for all criteria and a minimum (W = (1, 0, ..., 0)) and maximum (W = (0, 0, ..., 1)) risk implies an absence of trade-off. Therefore, risk and trade-off are not independent and they are completely determined by the order weights distribution, the risk by its skewness and the trade-off by its dispersion [23].
The relationship between order weights, risk and trade-off forms a decision-strategy space usually represented by a triangle to highlight the inconsistency of certain couple of risk and trade-off values. If we express the level of risk and trade-off as a couple of values (r, t) ∈ [0, 1] 2 , the three vertices of the triangle are represented by the three main configurations: low risk with no trade-off (r = 0, t = 0), high risk with no trade-off (r = 1, t = 0) and medium risk with full trade-off (r = 0.5, t = 1).

Experimental design and cluster analysis
Except for some trivial cases, like the ones described above, the limits of the decision-strategy space remain unclear, making an efficient exploration very complicated. This is due to the fact that, in practice, the generation of order weights according to a certain level of risk and trade-off is not formally established. [10] recently proposed a methodology to generate order weights using truncated distributions. The method allows for an automatic generation of order weights according to a certain level of risk and trade-off based on the two first moments of the order weights distribution. More specifically, each suitable couple of (r, t) values is associated with a truncated normal distribution with an average r and a standard deviation proportional to t. The probability density function is then discretized to obtain n order weights reflecting the predetermined level of risk r and trade-off t. It is important to note that for certain risk and trade-off values, no suitable truncated normal distribution can be found. This allows for a formal delimitation of the decision-strategy space, which takes, in this case, the form of a parabola whose equation is y = 4x(1 − x) (Figure 2). See [10] for more details.
In this work, we draw upon these recent advances in order weights determination to conduct a sensitivity analysis to assess the impact of the level of risk and trade-off in MCDA. The guiding idea is to efficiently explore the decision-strategy space in order to identify clusters of risk and trade-off values that leads to similar final suitability maps. As illustrated in Figure 2, we develop on a five steps approach to reach this goal. As it is usually the case in MCDA, we start with the selection and weighting of criteria with a panel of experts selected for their knowledge and heterogeneous point of view regarding the problem (urban land use allocation in our case). We then design an experimental plan to automatically generate order weights vector W covering the decision-strategy space. In this study, we generate 1, 000 vectors of order weights associated with 1, 000 risk and trade-off values. A suitability map is built for each vector of order weights using the OWA operator introduced in Equation 1. Then, the suitability maps are compared by measuring the euclidean distance between pixel values for each pair of suitability maps. We finally obtain a dissimilarity matrix used to perform a clustering analysis and identify clusters of risk and trade-off values. An average suitability map and its associated standard deviation is finally computed for each cluster.

RESULTS
In order to get a better grasp of the situation we plot in Figure 3 suitability maps obtained with typical vectors of order weights. We observe that the low risk attitude (W = (1, 0, ..., 0)) is very conservative and identify very few suitable areas (most of the  pixel show a suitability lower than 0.5). The area identified as not suitable are mostly composed of natural vegetation. Nevertheless, this extreme position toward risk attitude allows for the identification of areas that are suitable to urban development whatever the considered criterion. There are indeed several pixels with a high suitable score (higher than 0.5) that are mainly close to built-up area. Conversely, an extreme risk taking attitude (W = (0, 0, ..., 1)) tend to make the whole studied area suitable to urban development. This is due to the fact that for a given pixel one can always find a criterion suitable to urban development. These two configurations are not very nuanced since, on a pixel basis, the information is only coming from one criteria. Results obtained with the weighted linear combination that corresponds to an intermediate risk with full trade-off (top of the triangle) take better account of the heterogeneity of criteria. Note that in this particular configuration, the order weights have no impact on the process. We observe, then, that this full trade-off configuration leads globally to similar spatial features than the ones obtained in the low risk case, but with a level of suitability that is globally higher. However, beyond visual similarities between relative spatial suitability distributions, it is worth noting that the relative suitability of certain areas as compared with others is not the same in the two cases. Indeed, pixels with the same level of suitability in the low risk attitude case can have very different levels of suitability in the full trade-off configuration (see Figure S2 in Appendix for more details).  This first result provides informative evidence of the importance of comparing the impact of risk and tradeoff values on suitability maps in MCDA. This is why it is crucial to efficiently explore the decision-strategy space for a more accurate interpretation of the results obtained in MCDA. Following the method described in the previous section, we build a dissimilarity matrix between 1, 000 suitability maps associated with 1, 000 couple of risk and trade-off values drawn at random within the parabolic decision-strategy space. In order to identify potential cluster of risk and trade-off values leading to similar suitability maps, we apply an ascending hierarchical clustering (AHC) algorithm on this dissimilarity matrix using the Ward's agglomeration method. We choose the number of cluster based on the evolution of the ratio between the within-group variance and the total variance us a function of the number of clusters (Figure 4). We detect four clusters or groups of risk and trade-offs values that share similarities in the suitability maps that they produced. We then compute the average suitability maps associated with each of the selected clusters. The results are shown in Figure 5. We first observe that the segmentation of the parabolic decision-strategy space is mainly driven by the level of risk. Although there is a difference in trend across pixels, the global level of suitability increases monotonously with the risk. The segmentation obtained with high levels of hierarchical aggregation reflects therefore the global structure induced by the reordering of criteria. However, it must be noticed that when considering a higher number of clusters a trade-off based segmentation highlights differences among pixels for a given global level of suitability (see Figure S3 in Appendix for more details). The average suitability maps associated with the four clusters enable us to have a more nuanced view of the situation in the Thau region that with the typical sets of order weights (Figure 3). The low risk cluster provided a very conservative situation but bring into light the most suitable areas (the only ones with a correct suitability) that have a potential to be urbanized first within a sustainable planning framework. The low-medium risk cluster releases a large proportion of the agricultural land-use that is still suitable for urbanization. Conversely it allowed to identify the most important agricultural areas. The map of the medium-high risk cluster identifies almost all the agricultural areas suitable for urbanization. Ecosystems that provide a maximum of ES such as forests, garrigues or salt-meadows are the only areas that are considered not suitable using the approach presented. Eventually, the high risk cluster provided a map quite similar to the extreme risk situation, but allows to detect some critical ecosystems (mostly forests) of major importance in the region.
Although the total within-group variance obtained with four clusters is relatively low (Figure 4), it is also important to evaluate how this variability is spatially distributed. Notably, this variance based analysis will allow us to assess the impact of the trade-off for a given level of risk. We observe in Figure 6 that the spatial distribution of standard deviation of the suitability varies from one cluster to another. The uncertainty seems to be correlated with the suitability. However, the nature of the relationship changes according to the type of cluster. The uncertainty tends to increase with the suitability for the low and medium risk cluster, while the opposite behaviour is observed for the high risk cluster.

DISCUSSION
In this paper, we presented an approach to assess the impact of order weights and their associated level of risk and trade-off in multi-criteria decision analysis. The developed approach aims at efficiently explore the decision-strategy space formed by the relationship between risk and trade-off. Although the results might vary according to the formula used, a vector of order weights, whatever its size, can be summarized by its level of risk and trade-off. However, the opposite is not always true, making it complicated to automatically generate order weights according to a certain level of risk and trade-off and conduct rigorous sensitivity analysis. In this paper, we relied on the  recent advances made in order weight determination proposed in [10]. The method presented in this paper proposes to automatically generate order weights using truncated distributions according to a certain level of risk and trade-off based on the two first moments of the order weights distribution. This methodological framework enables us to generate a suitability map for each couple of suitable risk and trade-off values through an OWA operator, opening the possibility to perform systematic sensitivity analysis in MCDA.
We illustrated the potential of our approach on a case study of land use management in the territory of Thau (France). We selected this example to propose practical solutions to the common problem that is experiencing the south of France regarding urban sprawl, but also for its relevance as a concrete application of ecosystems' valuations using MCDA. The idea here is to prioritize suitable areas for urban development in the Thau territory, while preserving optimally "green" areas and their associated ES to promote a Standard deviation Figure 6. Spatial distribution of the within-cluster variability.
sustainable regional management. Indeed, one of the biggest challenges we are facing nowadays is to find suitable new areas for development in a sustainable way, without degrading existing valuable land uses. Spatial planning for urban development is directly related to decision-making and government strategies, besides outcomes of planning will clearly differ depending on stakeholders and territorial policies. The ecosystem service approach has the potential to widen the scope of traditional landscape-ecological planning to include ecosystem-based benefits, including social and economic benefits, green infrastructures and biophysical parameters to directly strengthen the role of landscape-ecological planning in urban and territorial planning [24,25]. The results obtained in this study demonstrates that a complete exploration of the decision-strategy space enable us to get a better understanding of how the different ecosystem services combine to produce the final suitability map.
In particular, the clustering analysis led to an in-formative segmentation of the decision-strategy space that helped us to really understand the role of order weights and particularly their impact on the trade-off between criteria. Our approach has the great advantage to be generic and relatively easy to implement in any multi-criteria evaluations. Although further analysis on more case studies are certainly needed, we believe that the approach developed in this study could provide an operational and collaborative tool to facilitate decision making processes to improve land use planning. This is particularly important as the final suitability map obtained in such process might be used to take important and irreversible decisions regarding land use management. Therefore, it is crucial to provide robust and quality-assessed materials to facilitate the dialogue between stakeholders and decision-makers.

ACKNOWLEDGMENTS
This study was partly supported by the IMAGINE project (ERANET BIODIVERSA 3) and the French National Research Agency (project NetCost, ANR-17-CE03-0003 grant). A special thank goes to Pierre Maurel, Julian Le Viol and the members of the SMBT of Thau. We acknowledge Marine Le Louarn for providing us with the landscape connectivity map.

DATA AVAILABILITY
Code and data are available at www. maximelenormand.com/Codes AUTHOR CONTRIBUTIONS ML coordinated the study, designed the study and wrote the paper. OB and MS wrote the code, processed and analyzed the data. SL coordinated the study. All authors read, commented and validated the final version of the manuscript.

Biophysical data
As mention in the main text, we added external biophysical data (Table 1) to adjust the ES valuation. For each pixel, the experts scores associated with the ES has been multiply by the factor extracted from the biophysical data as described below.
• The soil quality is an index made by the French National Institute for Agronomic Research and based on two main constraints (salinity and the available water capacity) and secondary constraints (soil sealing, hydromorphy, stoniness and pH). It has 16 categories so we assigned a multiplicative factor 1 to the most fertile and decreased of 0.05 for each category to end at 0.25 for salted ground.
• Habitats of ecological importance are defined areas by naturalists studies such as ZNIEFF I (Natural zone of ecological interest in terms of fauna and flora) or Natura 2000 (EU wide network of nature protection areas). We applied factor 1 when the pixel was located inside these protected area and 0.75 otherwise.
• Flow accumulation follow a hydrologic model of Multiple Flow Direction Algorithm. The value are converted to a 0-1 continuous scale. To avoid a crushed scale, we used 98% of the maximal value to set 1. Values above are also set to 1.
• We used a Slope Length Factor used by the Universal Soil Loss Equation (USLE) to adjust the mass stabilization and control of erosion ES. Here also, the value are converted to a 0-1 continuous scale. To avoid a crushed scale, we used 98% of the maximal value to set 1. Values above are also set to 1.
• The flooding hazard came from the official and mandatory risk prevention plans for inundations. It has three categories (high, medium, none) that corresponds respectively to a factor 1, 0.5 or 0.
• Distance to the main road was computed on SAGA GIS, taking into account of the relief. It is a proxy of accessibility, we used the threshold of 300m as the maximum distance for a recreational walk (5 min walk) to defined three categories: ≤ 300m, 300m-1km, ≥ 1km. Since almost the whole territory is suitable for hiking (no high elevation and slopes) we stopped the scale at 0.5 for a distance higher than 1km (around 15min walk) and computed a decreasing and continuous scale from 1 to 0.5 between 300m and 1km.
• The fire hazard has been derived from forest stands combustibility and spreading factors such as slope and speed propagation in a study carried out by the National Forest Office (ONF) and the Sea and Territory Departmental Direction (DDTM). For 6 values (very high, high, medium, low, very low, none) we started at a factor 1 for a very high risk and decreased of 0.1 for each category to end at 0.5 (unlike the flooding hazard, we considered that there was always a fire risk).
All the criterion maps are available in Figure S1.