Land use optimization through bridging multiobjective optimization and multicriteria decision‐making models (case study: Tilabad Watershed, Golestan Province, Iran)

This study aims to present an efficient methodology for land use optimization based on minimization of runoff and sediment and maximization of economic benefits, occupational opportunities, and land use suitability in the Tilabad watershed in northeast of Iran. The land use map of the area was prepared using the Landsat satellite images and field surveys. The amounts of runoff and sediment were estimated via SWAT model. The TOPSIS multicriteria decision‐making (MCDM) approach was applied on the results of the multiobjective optimization (MOO) based on non‐dominated sorting genetic algorithm II (NSGA II) to choose the final optimal solution among the Pareto solutions front generated by MOO. The results indicated that the area of agriculture and rangelands should decrease, and the area of forests should increase to achieve the defined objectives. Overall, results indicated that integration of MOO and MCDM provides an efficient procedure for land use optimization in a complex watershed.


| INTRODUCTION
Analysis and evaluation of land use patterns are the main objectives of natural resources management (Shaygan et al., 2014). These patterns are constantly changing due to increased manipulation by human-beings aiming to fulfill their various demands (Chakrabarty, 2001;Li & Liu, 2008). Rapid population growth and urbanization are the main driving forces of these changes. The consequences of these changes are reflected as hydrological changes and (Deng et al., 2015), ecosystem services perturbation, and land degradation (Grimm et al., 2002;Sunandar et al., 2014). Therefore, efficient land use management is necessary for managing these changes and creating sustainable development in society (Sharmin, 2017). Land use management refers to the process of allocating different standard land classes, such as industry, agriculture, forests, and so forth, to different land units consistent with the objectives of managers and decision-makers (Stewart et al., 2004). Multiobjective planning of lands is considered to be an efficient approach for sustainable use and protection of the environment. It provides an opportunity for watershed managers and decision-makers to make the best decision among different land-use alternatives (Min et al., 2015;Verburg et al., 2013;H. Zhang et al., 2016). Spatial planning is, in fact, a decision-making process related to the allocation of an activity among a set of activities or uses. It is mostly focused on a specific unit of land in which not only the extent of use but also the location of use is identified (Cao et al., 2012;Haque & Asami, 2014;Min et al., 2015;Stewart et al., 2004;Zhou, 2015). The existing problems in natural arenas are formulated via optimization methods and based on the relationship between different factors which leads to the best profit for residents of the watershed (Shively & Coxhead, 2004). Decision making in land use planning and protection of natural resources is a complex issue due to different and often contradictory objectives and extensive scope of the answers (Datta et al., 2007;Y. Liu et al., 2015). Application of computerized tools for computation and presentation of multidimensional and spatial data, such as geographical information system (GIS), is a great help for relieving these complexities. To increase the efficiency of GIS in solving complex spatial problems, intelligent algorithms such as multiobjective evolutionary optimization algorithms are used (Herzig, 2008). The main challenge in solving multiobjective optimization (MOO) problems is that a set of decisions may not optimize all the objectives simultaneously. In these cases, the optimization process results in a set of optimized solutions instead of a unique solution. Based on the nature of the problem, one of the mathematical methods, such as linear, dynamic, or nonlinear programming, can be used for optimization. The application of different optimization methods has expanded in recent years which resulted in comprehensive and logical management measures and tools.
The research carried out in the field of land use optimization are divided into three categories of linear programming models (LP) (Benli & Kodal, 2003;Hongrui et al., 2010;Min et al., 2015;Nikkami et al., 2009;Owji, Nikkami, Mahdian, et al., 2012;Sadeghi et al., 2009;Sokouti & Nikkami, 2017;Sunandar et al., 2014;Vafakhah & Mohseni Saravi, 2011), Cellular Automata (CA) models (Li & Yeh, 2002), and models based on intelligent algorithms, and evolutionary and innovative methods (Cao et al.,2012;Datta et al., 2007;García et al., 2017;Lazoglou et al., 2016;Y. Liu et al., 2013Y. Liu et al., , 2015Mohammadi et al., 2015;Shaygan et al., 2014;Yuan et al., 2014;H. Zhang et al., 2016;H. H. Zhang et al., 2010). In LP models, the optimal land use structure is identified rapidly based on the accurate mathematical definition of limits and one objective function (Aerts et al., 2003). However, LP models cannot be used for spatial optimization of land use (Santé-Riveira et al., 2008). To fill this gap, the CA models have been developed based on the rules of conversion and change of land use in local scales (Li & Yeh, 2002). Various and contradictory objectives cannot be easily included in CA models. Parallel to increasing the optimization space, the search for problem-solving becomes more complex and the calculations are intensified and time-consuming (Xiao, 2008), particularly for spatial analyses. Thus, various intelligent algorithms of spatial land use optimization such as genetic algorithm, simulated annealing algorithm, particle swarm optimization algorithm, and colony algorithm have been expanded as a compromise between timing and optimality of different final solutions (X. Liu et al., 2012). These methods provide a new approach for solving the problem of land use allocation and present multiple solutions to decision-makers (Cao et al., 2011). The nondominated sorting genetic algorithm II (NSGA II) evolutionary algorithm, as an optimization and genetic algorithm (Horn et al., 1999), is considered to be one of the most efficient and widely used multiobjective evolutionary algorithms. This algorithm does not require any additional parameters and has a high calculation speed (Deb et al., 2002), so it has been widely used in land use optimization and spatial planning (Cao et al., 2011;Datta et al., 2007;García et al.,2017;Hojjati et al., 2018;Lazoglou et al., 2016;Y. Liu et al., 2013;Mohammadi et al., 2015;Noor et al., 2016;Shaygan et al., 2014). However, this algorithm generates a set of nondominant solutions from which the selection of a unique solution becomes questionable for decisionmakers. Therefore, in this study, we used a multicriteria decision making (MCDM) technique to help decision-makers to select the final unique solution based on the running preferences.
Since watershed management aims utilization and protection of watershed resources in a sustainable manner, the implementation of optimal solutions for land use planning considering all the technical, social, economic, and legal potentials, limits, and preferences is imperative. Since a watershed is a complex spatial system, it is inevitable to use appropriate approaches and tools that can cope with the complexities of the system, multiplicity and the contradiction of objectives and preferences of various stakeholders and decision-makers in prescribing an optimal management solution. In this study a set of computer tools and applications including the ArcGIS spatial data analysis program, the SWAT (Soil and Water Assessment Tool) hydrological model, the Matlab based NSGA-II multiobjective genetic algorithm and the Microsoft Excel-based TOPSIS (Technique for Order of Preference by Similarity to Ideal Solution) MCDM technique has been used to optimize land use allocation in the Tilabad watershed in the northeast of Iran, where had undergone continuous land-use changes without clear and logical planning during the past decades. TOPSIS is a simple ranking method in conception and application and has features and advantages such as the simultaneous application of quantitative and qualitative criteria, weighting of criteria and ranking of decision variables (Opricovic & Tzeng, 2004).

| Study area
The Tilabad watershed, with a total area of 1540 km 2 , is located between geographic coordinates of 55°09′ 12″ to 55°50′ 40″ E longitude and 36°45′ 15″ to 37°15′ 10″ N latitude in the Golestan Province, northeast of Iran ( Figure 1). The weighted mean altitude of the watershed is 1146.32 m and the mean slope is 24.93%. Soils are mostly silty clay loam. The most common land uses include agriculture, broadleaf forests, rangelands, and bushes with area coverage of 45.23%, 28.45%, and 24.41%, respectively. The mean annual temperature, relative humidity, and precipitation are 14.16°C, 68.19%, and 597.13 mm, respectively. Despite considerable forest areas and relatively high precipitation, this watershed experiences frequent flooding events and a high rate of land degradation and sediment load (Sheikh et al., 2019). It is believed that land use planning and optimization can alleviate these threats. However, as far as our knowledge, there is no reported study regarding land use optimization and expected outcomes across this watershed.
F I G U R E 1 Location of the Tilabad watershed in Iran and Golestan province 2.2 | MOO Whenever an optimization problem includes more than one objective function, finding one or multiple optimal solutions for this problem is called MOO (Dias & De Vasconcelos, 2002). The decision space in single-objective problems is a sortable space; in other words, propositions which are defined by > or < are definable in a sortable area. However, decision space is not sortable in multiobjective problems (Shaygan et al., 2014). After identification of the problem in an MOO problem, a set of discovered decision variables are used in the optimization of objective functions concerning the limits of the problem. In these cases, finding a solution which optimizes all the objectives is nearly impossible and a set of solutions concerning all the values of the objective function is required (Coello Coello et al., 2007). In other words, the main objective is finding Pareto optimal solutions (Deb et al., 2002). The general form of multiobjective problems is as follows: An x solution is a vector of variable n: Each decision variable is subjected to a value between the lower limit of x i (L) and the upper . Moreover, f is a vector of the objective function m.  (3) g j (x) and h k (x) are constraints of a function, j stands for the number of limits of an inequality, and K stands for limit of a function. The mapping of the decision space into the objective space is indicated in Figure 2.
Finding the best solution for MOO problems, such as multiobjective land-use optimization, requires specific algorithms such as NSGA-II algorithm (Yuan et al., 2014). NSGA-II evolutionary algorithm is a kind of genetic and optimization algorithm which is widely applied for solving land-use problems (Cao et al., 2012). The NSGA-II algorithm was first presented by Deb et al., 2002. It utilizes elite attitudes and can search in an expansive space of decision and objective variables (Maringanti et al., 2009).

| MCDM
The MOO models provide many optimal solutions, which are equally good from the perspective of the given objectives. These solutions, known as Pareto optimal front and as nondominated solutions, provide multiple choices of equal priority for implementation (Wang & Rangaiah, 2017). Therefore, the selection of an optimal solution from the set of nondominated solutions requires a tool or methodology to avoid the arbitrary selection of the preferred solution. The TOPSIS is the most common MCDM method employed to select a preferred solution among a set of solutions (Wang & Rangaiah, 2017). TOPSIS is a simple ranking method in conception and application, which was originally developed by Hwang and Yoon (1981) as an MCDM method. According to this technique, the chosen optimal solution should have the smallest Euclidean distance from the ideal (also known as positive ideal) solution and also the largest Euclidean distance from the negative ideal solution. It has features and advantages such as the simultaneous application of quantitative and qualitative criteria, weighting of criteria, and ranking of decision variables or alternatives. The steps to follow in this method includes constructing a decision matrix (m × n), normalizing the decision matrix (n ij ), determining the weights of criteria and constructing normalized decision matrix (v ij ), determining the positive ideal solution (v j + ) and negative ideal solution (v j − ) according to each criteria, calculating distances of each alternative from the positive (d i + ) and the negative (d i − ) ideal solution using n-dimensional Euclidean distance, and calculating the relative closeness of each alternative to positive ideal solution (CL i + ) using the following equations: For positive elements of criteria } , For negative elements of criteria}, Decision variable space and objective function space in multiobjective optimization (Coello Coello et al., 2007) ∑ where, r ij is a normalized matrix, a ij shows entries of the decision matrix. Therefore, in this study after identifying the Pareto optimal solutions per each land optimization scenario, an MS Excel-based MCDM program (Wang & Rangaiah, 2017) has been employed to reach the final optimal solution according to the TOPSIS technique and the Shannon entropy weighting method. Shannon entropy is a well-known objective mathematical weighting method and its main advantage is that it works with quantitative data that are more reliable. It is especially useful in obtaining the weights when obtaining reliable subjective weights based on the decisionmakers preferences is difficult (Hosseinzadeh Lotfi & Fallahnejad, 2010).

| Data
In the present study, digital elevation model (DEM), drainage network map, land use maps, meteorological data, and soil data were used as base maps and data required for simulation and optimization models. Furthermore, the simulated runoff and sediment by the SWAT model were used as objective functions in the optimization process.

| Land use map
The land use map was prepared from OLI (Operational Land Imager) image of Landsat satellite in 2015, through processing it in the ENV14.8 software using maximum likelihood supervised classification method. The map included seven land-use classes, including waterbodies, coniferous forest, dense deciduous forest, sparse deciduous forest, agriculture, rangeland, and residential areas ( Figure 3). As seen, croplands (43%) and rangelands (29.9%) are dominant land-use types, whereas rangelands cover the uplands in south and croplands, although scattered throughout the watershed, but are mainly located in downstream near the watershed outlet. The deciduous forests (25%) are located in the middle parts of the watershed.

| Physical criteria
Runoff depth and sediment rate were used as the physical objectives and criteria for MOO and MCDM. To determine runoff and sediment values for different land-use types and scenarios, the SWAT model was used with 30 m resolution SRTM DEM and drainage network map. Based on soil, land use, and slope maps, the watershed was divided into 430 hydrological response units (HRU) (Figure 4). SWAT was run in monthly time steps for the period of 2000-2008 using the observed data in the hydro-meteorological stations within and around the watershed (Figure 4).
The initial 2 years of 1998 and 1999 were considered as warming-up period. The model was calibrated automatically using SUFI-2 software under the SWAT-CUP software package and observed discharge data of Arazkuseh hydrometric station at the watershed outlet. After optimization of parameters set in the calibration period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008), the model was run for the validation period of 2009-2014. The degree of uncertainty was measured via P-factor and R-factor. Moreover, the performance of the model in simulating flow discharge and sediment load was evaluated using efficiency criteria of coefficient of determination (R 2 ), Nash-Sutcliffe efficiency (E NS ), root mean square error (RMSE), and percent of bias (PBIAS) ( Table 1). Runoff and sediment values were calculated for each HRU and land-use scenarios ( Table 2). As seen, most of the model performance evaluation metrics indicate a relatively good performance in simulation of both discharge and sediment variables. Whereas, R 2 is always above 0.6 and also NS is above 0.6 except for the validation set of the sediment variable. P-factor, which is an indicator of the uncertainty of the model prediction, is always above 0.9. Moreover, the differences between mean of the observed and simulated values of discharge and sediment is less than 5%.

| Socioeconomic criteria
The ratio of benefit to cost and number of jobs created per hectare for land use types were obtained through field surveys, interview with local people and pertinent experts, and literature review as economic and social criteria, respectively (Table 3). Economic value per hectare of each main land-use type is equal to the total net production value of that land use, nonmarket services, and ecosystem services. To estimate economic value of the different land use types, the net present value was calculated with the discount rate of 20% over a time horizon of 10 years (minimum length of time that a forestation project requires intensive maintenance activities till it can self-sustain). The mean annual timber production of forests has been assumed 7.4 m 3 /ha and mean annual production of dry foliage by rangelands was estimated 312 kg/ha. For F I G U R E 3 Land use map of the study area for year 2015 calculation of nonmarket services, and ecosystem services by different land use types the reported values in the literature were consulted (De Groot et al., 2012;Hossein et al., 2017;Lv et al., 2010;Ninan & Inoue, 2013;Rastgar et al., 2013;Schaubroeck et al., 2016;Xue & Tisdell, 2001

| Land-use suitability criteria
To determine the suitability of an area for each land-use type, 24 indices were identified including slope, aspect, altitude, susceptibility to soil erosion, geology, susceptibility to landslide, distance from faults, current land use, vegetation density (NDVI), vegetation type, distance from the road, distance from the residential area, distance from spring, distance from the river, distance from well, maximum and minimum annual temperature, relative humidity, precipitation, soil texture, soil EC, soil depth, soil pH, and soil organic matter. Then, maps of these indices were prepared within ArcGIS software. The indices were assigned weights via analytical hierarchy process (AHP). After standardization of the indices (fuzzification) and identification of the limits for each index, land use suitability maps for forests, rangelands, agriculture, and residential areas were produced within the IDRISI software through using multicriteria evaluation (MCE) approach and weighted linear combination (WLC) method ( Figure 5) (Esa & Assen, 2017;Nahusenay & Kibebew, 2015;Salman-Mahini & Gholamalifard, 2006). Values of suitability maps vary in the range of 0-255, the higher the pixel value is the more suitable is for the considered land-use type ( Figure 5).

| Mathematical model of land use optimization
The first step in optimization is the formulation of an appropriate mathematical model, including decision variables, objective functions, and a set of constraints. Components of this mathematical model in this study are as follows:

| Decision variables
Land use types were used as decision variables.
LU members from 1 to 6 stand for coniferous forests, dense deciduous forests, sparse deciduous forests, agricultural lands, rangelands, and residential areas.
In these equations, E ij , M ij , R ij , S ij , and C ij indicate economic profit, social, runoff, sediment load, and land use suitability, respectively; and n ij = represents the number of HRUs.

| Constraints
The set of constraints includes: Where Slope i,j is slope of the unit (i, j), and Slope lu min and Slope lu max are slope constraints for lu land-use. Since minimum slope constraint for each land-use is zero, we will have: Where Altitude i j , is altitude of the unit (i, j), and Altitude lu min and Altitude lu max are altitude constraints for lu land-use. Since minimum altitude constraint for each land-use is zero, we will have: Land use conversion constraints: In spatial planning, the allocated area to different land-use types, and their convertibility are of great importance. On the other hand, the resulted chromosomes from the evolutionary method of land use types can be allocated by different areas. In this study, HRUs with dominant residential and dense deciduous forest land uses (according to inconvertibility of these land-use types) are removed from the optimization process. Therefore, the minimum residential and deciduous forest areas are fixed as the status quo. Table 4 indicates possible land-use convertibility matrix.
Conversion rate of land uses was limited to five scenarios of 10%, 40%, 30%, 20%, and 50% changes with respect to the current condition. Considering the inconvertibility of residential and dense forest land-use classes, only 339 out of 430 HRUs were used for optimization purpose. The mathematical formulation of the problem has been expressed as follows: where LnNow stands for the current land-use, and LnOptim stands for the optimal land-use. According to this constraint, 1 indicates equality of the optimal land-use with the current landuse. Otherwise, this value will be 0. Implementing this limit will cause minimum number of HRUs to be changed concerning the objectives.

| Implementation of the developed model
The whole study area was considered as a chromosome using the NSGA-II algorithm. Each HRU stands for a gene, and each gene has a vector of six elements indicating land-use types. If an HRU is allocated for specific land-use, its related element is 1, and other elements' values are 0. After generating the primary population, crossover and mutation operators were applied to evolve the population and obtain optimal solutions. For crossover operator, some chromosomes were selected as parents (in this study, crossover rate is 0.8; therefore, 16 parents were selected in each level based on the primary population of 20 chromosomes). Afterward, some genes from one of the chromosomes in each pair of chromosomes (170 genes) were selected and replaced with similar genes of the other chromosome. 10% of the genes (33 genes) were selected for mutation operator by considering all the dimensions of the problem and different scenarios with a specific percentage of land use changes for HRUs. According to the mutation ratio (0.3) and the total number of the population (20 chromosomes), six chromosomes were selected for mutation. The genetic algorithm has been used to generate primary chromosomes via a random approach (Beasley et al., 1993). However, due to the constraints of land-use change and the wide scope of the optimization problem (6 × 339 × 20), it takes time to generate the primary population. Therefore, some mutations were applied to some of the chromosomes of the current land uses to be used as primary chromosomes. The model was coded and processed within MATLAB software by taking the size of the primary population, crossover rate, and mutation rate as 20, 0.8, and 0.3, respectively, for all the runs. The number of generations was set to 20,000 and 20 solutions (optimal chromosome) were selected as Pareto set. Then the TOPSIS MCDM technique and the Shannon weighting methods were used to select the final optimum solution among the Pareto set. Furthermore, the percentage of changes in each objective function was estimated.

| RESULTS AND DISCUSSION
According to local authorities and experts' viewpoints, extensive unsustainable land-use changes are blamed as the major cause of many environmental degradation phenomena across the study area. Considering the current conditions of the study area, the scope of the dominant problems, and the expert's viewpoints, the optimization of land use patterns seems imperative to harness the increasing trends of land degradation and natural hazards risks (Dias & De Vasconcelos, 2002) in the Tilabad watershed. To this end, the first trends of land-use changes during three decades  have been calculated ( Figure 6). As seen in this figure, the extent of agricultural lands and residential areas has increased while the extent of forests and rangelands has declined during the past decades. Then, the NSGA-II genetic algorithm and the TOPSIS technique have been integrated to optimize the land-use allocation in the study watershed aiming to minimize the surface runoff generation and land degradation and maximizing land use suitability, economic profit, and job opportunities. The summary of the results has been presented in Table 5 and Figure 7.
The descriptive statistics of the Pareto solutions front including minimum, maximum, and mean values have been presented for each scenario and objective variable in this table. As seen, there is a considerable variation (the difference between the minimum and maximum) in the values of objective variables for the Pareto fronts which creates dubiousness in the decisionmaking process. Generally, by moving from one optimal Pareto solution to another, at least one objective function improves in lieu of worsening another objective function (Yee et al., 2003). Therefore, the MCDM-TOPSIS method has been used to prioritize the set of Pareto solutions and select the final optimal solution for the objective variables. All land-use change scenarios show a decline of surface runoff and sediment load. By increasing the level of scenarios, the absolute amounts of changes increase (Figure 8). However, the relative amounts of changes decrease as the level of land use optimization increases (Table 5).
The rate of changes in sediment load is higher than the surface runoff. Indeed, land-use optimization above the 10% level will not significantly decrease the surface runoff. As seen in Table 5 for the 10% land-use optimization scenario, surface runoff decreases as much as 20%, but by increasing the level of land use optimization from 10% to 50%, it decreases only about 5%. Therefore, it seems that hotspot areas regarding surface runoff generation and flood formation covers limited areas of the watershed. Among the five objective functions, the sediment load reduction has been highly influenced by the land use optimization scenarios, whereas it has declined from more than 16 ton/ha/year for the present condition (base scenario) to about 2 ton/ha/year for the 50% scenario, resulting into as much as 87% reduction (Table 5). The economic profit and land use suitability have continuously increased by increasing the level of land use optimization scenarios. The economic profit and land use suitability have respectively increased from 5% to 38% and from 7% to 31% by increasing the level of land use optimization scenario from 10% to 50%. However, despite using the maximization function for the job F I G U R E 6 Changes in area of land-use types in the Tilabad watershed during past decades opportunities objective, it has continuously decreased by increasing the level of land use optimization scenarios. In fact, by improving the land use condition, the number of job opportunities per hectare of land area will decrease. As seen in Figure 7 and Table 5, the job opportunities number across the whole watershed will decrease from 8795 man-days for the base scenario to 6433 in the 50% land optimization scenario. Therefore, by land-use optimization, land-based job opportunities will decrease and more manpower will be available for other socioeconomic sectors. The results of this study are in agreement with the reported results of Owji and Nikkami (2012) Distribution of the Pareto front elements for different paired objectives has been presented in Figure 9 for two scenarios of 10% and 50%. The final optimal solutions generated by the TOPSIS technique have been highlighted in the Pareto front. As seen in this figure, the final optimal solutions are stable for all the objective criteria. Therefore, it can be concluded that the combination of MOO and MCDM models has resulted into an efficient and robust methodology for land use optimization across the Tilabad watershed. Comparing the scatterplot of two scenarios indicates that by increasing the level of land-use optimization, the variation and length of the Pareto fronts increases. This means that dubiousness and uncertainty in decision making for land use optimization process increases by increasing the extent of land-use optimization. Thus, it is strongly recommended to apply an MCDM model to the results of the MOO model, particularly where the level of the interventions is high. According to Wang and Rangaiah (2017), among different techniques of the MCDM model, the TOPSIS technique performs better than other MCDM techniques for selection of optimal solution from the Pareto solutions front.
After selection of the final optimal solutions from the Pareto fronts corresponding for each land-use optimization scenario, the allocated area coverage and percentage of different land-use types have been calculated and presented in Table 6. Furthermore, the optimized land-use map for each scenario has been prepared and presented in Figure 10.
The optimization procedure carried out in this study is a kind of a complex and thorough computational optimization process regarding the number of decision variables, scenarios, models, and techniques. Traditional optimization techniques, such as linear programming, are not suitable for solving such complex computational optimization problems (Thompson et al., 1973). Other approaches, such as integer programming and mixed integer programming, are only able to solve problems with limited complexity and are inappropriate for highly complex problems such as land use allocation across a complex watershed system (Bevers & Hof, 1999). Therefore, only exploratory optimization techniques, such as Tabu search (TS), gradual annealing, genetic algorithm, and so forth, could be used for land use optimization of the Tilabad watershed considering its various potentials and constraints (Cao et al., 2011;Venema et al., 2005). Since the T A B L E 5 Descriptive statistics (min, max, and mean) of objective functions before and after optimization for the Pareto solutions front using the MOO model and the objective values for the final optimal solution selected using the MCDM model NSGA-II genetic algorithm is a widely used exploratory optimization algorithm in solving spatial planning problems (Cao et al., 2011;Datta et al., 2007;García et al., 2017;Hojjati et al., 2018;Lazoglou et al., 2016;Y. Liu et al., 2013Y. Liu et al., , 2015Mohammadi et al., 2015;Shaygan et al., 2014;Yoon et al., 2019). Despite good capabilities of the MOO models in solving complex spatial planning problems, the focus of them is on finding the nondominated solutions and the next step of selecting one of the nondominated solutions from the Pareto front is not considered (Wang & Rangaiah, 2017). The complex MOO models provide many nondominant solutions from which selection of a final solution turns into a decision choice problem. Because, by moving from one solution to the other, the values of the objective variables change considerably and make the decisionmaker to hesitate on reaching a certain final solution, particularly when the importance (weightage) of objective variables changes. Therefore, application of the MCDM model on the results of the MOO model help the decision-makers to reach an optimal solution considering the preferences dominant at the moment of the decision-making process.
According to the final optimal solution, the area coverage of all forest land use types should at least increase more than 10% for the 10% scenario to more than 70% for the 50% scenarios. On the other hand, the area coverage of agricultural lands and rangelands should decrease considerably. The area of agricultural lands should decrease from about 70,000 ha at the current condition to at least about 60,000 ha (more than 12% reduction) for the 10% scenario and at most about 40,000 ha (more than 40% reduction) for the 50% scenario. The predicted changes in the rangeland area are negligible for the 10% and 20% scenarios. But it is predicted to decrease by about 31% for the 50% scenario. This might be due to the fact that computation of the value of all objective variables for the rangeland has been based on the currently poor condition of the rangelands induced by overgrazing. Otherwise, it is usually expected to have an increase in the area coverage of rangelands due to several merits of good quality rangelands. Due to limited absolute area coverage of coniferous forests, the predicted rate of changes for this land-use type is highest among all land-use types. The dense deciduous and sparse deciduous forest will continuously increase by increasing the level of land use optimization scenarios. These results are in agreement with the reported results of Owji, Nikkami, Mahdian, et al. (2012) and Sunandar et al. (2014) on decreasing the area of agricultural land and rangeland and increasing forest land after the optimization process. These facts and figures indicate that land-use allocation of the Tilabad watershed is far from optimality. This finding reinvigorates the F I G U R E 8 Effects of land use optimization scenarios on various objective functions F I G U R E 9 Scatter plots of different objective functions for 20 patterns selected from the Pareto front of two scenarios (10% scenario has been presented with black Asterix symbols and 50% scenario has been presented with blue triangles), the optimal patterns selected using the TOPSIS method of MCDM approach are shown with filled symbols per each scenario. MCDM, multicriteria decision-making; TOPSIS, Technique for Order of Preference by Similarity to Ideal Solution

| CONCLUSIONS
This study aimed to make a comprehensive investigation on land use optimization with main objectives of reducing surface runoff and sediment load and improving the socioeconomic and land-use suitability conditions. Appropriate land use prioritization and allocation in terms of place and extent comprise the main aspect of the spatial land use planning process. The significance of the issue is highlighted when various and even contradicting objectives of different competing beneficiaries across a large watershed lead to complexity in solving the issue. The geographical information system and its capabilities in spatial data storage, management, and analysis offer great help but the gap in obtaining more systematic and confident methods simplifying the decision-making process remains. In this regard, spatial data should be combined with efficient algorithms and models to analyze decision-making problems and relieve the complexities. In the present study, a combination of MOO and MCDM models has been applied to generate a set of nondominant solutions and prioritize them. The NSGA-II F I G U R E 10 Land use maps after optimization according to 10%, 20%, 30%, 40%, and 50% change scenarios algorithm has high flexibility and efficient performance in obtaining non-dominated Pareto solutions. However, its major drawback is that it generates many locally optimal solutions, restricting the decision-makers in reaching a consensus on the final global optimal solution. Therefore, we applied the TOPSIS technique on the results of the NSGA-II algorithm to resolve this drawback. A combination of MOO and MCDM models has resulted into a robust procedure supporting the decision-makers in selecting the best land-use allocation solution based on various objectives including physical, social, economic, and environmental criteria and management scenarios.
ACKNOWLEDGMENT Gorgan University of Agricultural Sciences and Natural Resources is acknowledged for providing an opportunity to carry out this study as a Ph.D. Thesis (Grant No. 9219054102) in the study field of Watershed Management.