Combining biophysical optimization with economic preference analysis for agricultural land-use allocation

Agricultural production provides food, feed, and renewable energy, generates economic profits, and contributes to social welfare in many ways. However, intensive farming is one of the biggest threats to biodiversity. Although current market forces and regulations such as the European Union’s Common Agricultural Policy, seem to foster agricultural intensification, a socially and ecologically optimal land-use strategy should seek to reconcile agricultural production with biodiversity conservation. Research on spatial land-use allocation lacks studies that consider both aspects simultaneously. Therefore, we developed a method that finds landuse strategies with a maximum contribution to social welfare, taking into account the landscape’s biophysical potential. We applied a multiobjective optimization algorithm that identified landscape configurations that maximize agricultural production and biodiversity based on their contribution to social welfare. Social welfare was approximated by the profit contribution of agricultural production and society’s willingness to pay for biodiversity. The algorithm simultaneously evaluated the biophysical outcomes of different land uses using the Soil and Water Assessment Tool (SWAT) and a biodiversity model. The method was applied to an agricultural landscape in central Germany. The results show that, in this area, overall social welfare can be increased compared to the status quo if both social benefits from biodiversity and economic profits from agricultural production are considered in land-use allocation. Further, the resulting optimal solutions can create win-win situations between the two, usually conflicting, objectives. The integration of preference information into the biophysical optimization allows reducing the usually large set of Pareto-optimal solutions and thus facilitates further stakeholder-based analyses. Our explorative study provides an example of how socioeconomic data and biophysical models can be combined to support decision making and the development of land-use policies.


INTRODUCTION
Traditional land uses have been shaping the European countryside for millennia. They have created a diversity of cultural landscapes that has fostered biodiversity (Vos andMeekes 1999, Plieninger et al. 2006). Today, European agrobiodiversity is considered a valuable resource that needs to be protected (European Learning Network on Functional Agrobiodiversity 2012). However, in the mid-20th century, agricultural production became more intensive and has since led to a dramatic and still ongoing decline in biodiversity (Donald et al. 2001. Farmland birds, which are commonly used as an indicator for biodiversity in agricultural areas, have seen the highest decrease in population sizes compared to other bird species (Inger et al. 2015). In Germany, for instance, one reason for the decline in open-land bird species is the loss of grassland. Additionally, large fields and fast-growing, tallstatured crops hinder ground-breeding birds during the breeding season, and the use of pesticides reduces the available food (Hötker et al. 2014). Although agricultural production is considered one of the biggest threats to biodiversity (Behrman et al. 2015), it also provides food and economic benefits to humans. However, biodiversity is an important factor for the proper functioning of various ecosystem services (ESS), particularly for agricultural production (European Learning Network on Functional Agrobiodiversity 2012). One of the most prominent ESS examples is pollination (Hass et al. 2018). A high level of species diversity in agricultural systems also enhances resistance to pests, adaptability to changes in the system, and resilience. These functions are of special importance in the face of climate change and the resulting climate variability (Frison et al. 2011).
Because of the conflicting societal demands on land for agricultural production and biodiversity, it is of growing interest to find land-use strategies that enhance landscape multifunctionality (i.e., the provision of multiple ESS within a landscape) and heterogeneity (e.g., diversity of habitat types, arrangement of habitat patches; Hass et al. 2018, Hölting et al. 2019) so that landscapes are used and shaped in a way that reconciles ecological and socioeconomic objectives. So far, only a few land-use optimization studies simultaneously take into account biodiversity and socioeconomic aspects (see, for example, Polasky et al. 2008, Butsic and Kuemmerle 2015, Verhagen et al. 2018). Cavender-  emphasize that for managing socialecological systems, it is necessary to understand both the biophysical trade-offs between different ESS and how these ESS contribute to the well-being of different stakeholders. Analyzing trade-offs between ESS helps in finding ecologically meaningful land uses and land-use allocations. These analyses are often done by applying multiobjective optimization methods that lead to a whole set of "optimal" land-use allocation strategies (Kaim et al. 2018). To find applicable solutions for real-world implementation, stakeholders are increasingly involved in the decision-making process (Memmah et al. 2015, Lienhoop andSchröter-Schlaack 2018). Stakeholder inclusion can be achieved by considering their preferences at different stages of the optimization, for example, Fig. 1. Methodological framework for the analysis. The biophysical models for agricultural production and biodiversity are coupled with their respective willingness to pay (WTP) and interact with the optimization algorithm. From the resulting set of nondominated solutions (Pareto frontier), the solution with the highest sum of both WTP categories creates the highest contribution to social welfare, taking into account the biophysical potential of the landscape.
with previous weight setting, interactive and iterative adjusting of objective functions and constraints, and by applying methods from the field of multiattribute decision-making afterward (for more information, see Kaim et al. 2018). A similar approach has been put forward by Cavender- , who present an ecological-economic framework that combines utility functions of different stakeholders (indifference curves) with the biophysical Pareto frontier. They illustrate the method by giving an example that identifies land-use strategies for agricultural production and biodiversity that are most desirable for stakeholders (see also King et al. 2015).
Here, we develop an approach that combines methods from ecology and economics with multiobjective optimization. In contrast to Cavender- , we use preference information (willingness to pay for biodiversity and contribution margins for agricultural production) directly within the optimization process, instead of using utility functions after the optimization. This way, the applied algorithm identifies optimal land-use strategies that maximize social welfare (i.e., an aggregate monetary measure of preference satisfaction), given the available preference information, while simultaneously taking into account the biophysical potential of the landscape. As a proof of our concept, we demonstrate a simplified real-world application of the method to a river basin in central Germany.
We next define social welfare for the purposes of this study and explain how we approximate. We then provide insight into the optimization process and the calculation of a final optimal solution. Finally, we present the application of our method to a case study area and discuss the results.

METHODS
Our approach combines biophysical optimization with preference information, following the economic concept of social welfare. Social welfare is an aggregate function of individual welfares, i.e., expressions of "the value [each individual] attaches to his personal circumstances in a social state" (Dasgupta 2001:14). It is an aggregate measure of preference satisfaction, including preferences for ecosystem services and biodiversity (Millennium Ecosystem Assessment 2005).
Here, we focus on the contributions of agricultural production and biodiversity to social welfare. We approximate social welfare by combining aggregate willingness to pay (WTP) for biodiversity and contribution margins for agricultural production. Both factors are measures of opportunity costs of the two goods in question (biodiversity and agricultural production).
In the methodological framework of our analysis (Fig. 1), agricultural production and the level of biodiversity are both modeled by means of biophysical models and linked to the respective WTP. This information is then passed on to an optimization tool that evaluates different land-use configurations. The result is a set of mathematically optimal land-use configurations (Pareto frontier; see Methods: Optimization for a detailed definition). For each point on the frontier, i.e., each specific land-use strategy, we calculate its contribution to social welfare by summing the respective WTP for both objectives. Consequently, the land-use strategy with the highest sum is the solution with the maximum contribution to social welfare.

Modeling social welfare
We next describe the details of how the contributions of agricultural production and biodiversity to social welfare are modeled by combining the respective biophysical models with socioeconomic data (yellow and green boxes in Fig. 1).

Social value of agricultural production
We model agricultural production using the Soil and Water Assessment Tool (SWAT; Arnold and Fohrer 2005). For the simulation period from 1995 to 2009, SWAT calculates the annual amount of harvested biomass for each potential land-use type in each hydrological response unit (HRU; i.e., the smallest spatial model unit defined as a unique combination of soil type and land use per sub-basin). Information on model input data and performance is provided in Appendix 1. The social value of agricultural production (an approximation of the WTP for it) is assumed to equal the annual contribution margin (€/ha) averaged over each HRU and the entire simulation period. HRU-level contribution margins are estimated by multiplying the simulated biomass yield by crop-specific market prices less the variable costs, e.g., for field activities and fertilizer, as obtained from standard agronomic data of the Association for Technology and Structures in Agriculture (KTBL: http://www.ktbl.de/). Subsequently, this information is transferred to a look-up table (for code and input data see Jungandreas et al. 2020). It should be taken into account Ecology and Society 26(1): 9 https://www.ecologyandsociety.org/vol26/iss1/art9/ that, because of the focus on agricultural landscapes, we do not consider any contribution margins for forest.

Social value of biodiversity
To estimate the social value of biodiversity for a specific landscape, it is first necessary to derive the level of biodiversity in the area. We use a bird habitat suitability model and then combine its output, a bird habitat index, with the WTP for the respective biodiversity level. To describe the relationship between different values of this bird indicator with the respectively changing WTP, we developed WTP functions for biodiversity. These functions are also used in the optimization process to evaluate the social value of biodiversity for changing landscapes.
We apply a biodiversity model that has been developed by Anne Jungandreas in the context of the project "Towards multifunctional agricultural landscapes in Europe -TALE" (https://www.ufz.de/tale/). The model has been tested for the Middle Mulde River basin in Germany, of which the case study area is a sub-basin (see Methods: Application). More information about the model can be found on the TALE learning environment (http://tale.environmentalgeography.nl/wp2-indicators-and-models/) and in the presentation by Jungandreas et al. (http://tale. environmentalgeography.nl/wp-content/uploads/2020/03/ Gfo2018_JungandreasA.pdf). Furthermore, the model, as used here, is available on GitHub (see Jungandreas et al. 2020). The biodiversity model calculates an index based on the percent change in suitable habitat for different bird species compared to the initial land use of the respective area (in our case, the status quo). For each species, it runs a random forest model that is driven by spatially distributed data on land use and land cover, climate, and soil. We selected nine Red List bird species (Table 1) assuming that the habitat suitability for these species is positively linked to the habitat suitability of more common birds with similar requirements. For the status quo, the indicator always returns a value of 1 because it serves as a reference for the change in suitable habitat for different land uses. Thus, if the land use changed and the model returns an indicator value that is > 1 (e.g., 1.3), the amount of suitable habitat increased (by 30%) compared to the initial land use; if the indicator value is < 1 (e.g., 0.8), suitable habitat decreased (by 20%). We obtained data about WTP for biodiversity in Germany from the study by Hirschfeld et al. (2021). In 2013, they carried out a choice experiment as part of a population survey with citizens from across Germany. In total, 8800 questionnaires were completed. The study was considered representative; only salary and education level of respondents were slightly above average. In the survey, respondents could select from among scenarios that described how the landscape should develop within a 15 km radius of their area of residence. These scenarios included various landscape features such as "share of forest in the landscape", "size of fields and forests", "share of corn production in agricultural crop land", etc., that were considered cultural ecosystem services. Additionally, they included two biodiversity features: "biodiversity in forests" and "biodiversity on agricultural land".
The results were used to identify the individual WTP for the different features, of which we used only WTP for biodiversity on agricultural land. These monetary values were defined as a yearly financial contribution to a landscape fund. Changes in biodiversity that the respondents had to evaluate were expressed by means of a bird indicator developed by the Federal Agency for Nature Conservation (Bundesamt für Naturschutz -BfN; Dröschmeister andSukopp 2009, German Government 2017). This indicator uses birds as an umbrella species for overall biodiversity (plants and animals) for different landscape types. The biodiversity level is evaluated using a scoring scale: A score of 100 points serves as a reference value for the desired level of biodiversity. Thus, if an area gets a score ≥ 100 points, it is deemed particularly suitable for typical animal and plant species. In Hirschfeld et al. (2021), the status quo indicator value was assumed to be 65, according to the indicator values for agricultural landscapes by the BfN. The results of the choice experiment showed that people were willing to pay 22 € person −1 yr −1 and 53 € person −1 yr −1 to reach a score of 85 and 105 points, respectively.
Because Jungandreas' bird model is applied in the optimization (see Methods: Optimization), the bird indicator values need to be aligned to the biodiversity indicator values from Hirschfeld et al. (2021). Therefore, we develop functions that assign each biodiversity indicator value (based on Jungandreas) to the WTP for the respective biodiversity level. This alignment can be done because both bird indicators (from Hirschfeld et al. 2021 and Jungandreas) are based on suitable habitat and we therefore consider them comparable. Furthermore, we assume that the status quo biodiversity index of 65 used by Hirschfeld et al. (2021) is representative of the current status in the case study area. This assumption is based on a study from the State Office for Environmental Protection Saxony-Anhalt (Landesamt für Umweltschutz Sachsen-Anhalt 2015) that used the same bird indicator as Hirschfeld et al. (2021) and found similar results for the status quo. Saxony-Anhalt is geographically close to the case study area in the northern part of Saxony and has a similar landscape (see Methods: Application).
The first step is to translate the biodiversity indicator used by Hirschfeld et al. (2021), for example, 85, into the indicator developed by Jungandreas, which will be used throughout the optimization process. This transformation is done by applying the mathematical "rule of three": for the status quo, 65 points (Hirschfeld et al. 2021) equates to an index of 1 (Jungandreas). Taking this latter into account, the transformation of 85 points (Hirschfeld et al. 2021) is as stated in Eq. 1, with [HF] referring to the index by Hirschfeld et al. (2021) and [JA] to that by Jungandreas.
Ecology and Society 26(1): 9 https://www.ecologyandsociety.org/vol26/iss1/art9/ ⇔ ⇔ = = =^^ Additionally, the monetary values for WTP must be expressed in the same units as those that represent the contribution margin of agricultural production, i.e., € ha −1 yr −1 . WTP obtained by Hirschfeld et al. (2021) is in units of € person −1 yr −1 referring to the agricultural area within a 15 km radius. For simplicity, we assume that population and agricultural area are evenly distributed within Germany. Therefore, it is sufficient to multiply WTP by the number of people of legal age in Germany in 2013 (67.23 million; Destatis 2016, Hirschfeld et al. 2021, when the choice experiment took place, and divide it by the total amount of agricultural area in the same year (16.7 million ha; World Bank 2018). The respective calculation is given in Eq 2.

€ € ha ha
Eqs. 1 and 2 can be calculated analogously for a biodiversity score of 105. Linking the indices to the respective WTP, we obtain the results given in Table 2. These data form the basis for the development of the WTP functions. Because three data points are not enough to perform an interpolation, we define two different WTP functions to explore different options of how these data points could be related (Fig. 2). Furthermore, the functions should be economically reasonable: In economic theory, it is assumed that the marginal utility of a good, and thus the WTP for it, decreases for every additional unit that is provided until it converges at some point (Gossen's First Law;Gossen 1983). Therefore, both WTP functions follow the shape of a saturation curve. The level of biodiversity for the highest WTP (105 points) already slightly exceeds the goal of the BfN that has been set for the year 2030 (100 points). Additionally, no data about the WTP beyond this point are available. Consequently, we assume that saturation starts at a biodiversity index of 105/1.6. Both WTP functions are noncontinuous and defined piecewise as follows.
The aim of our study, i.e., to find the optimum combination of agricultural production and biodiversity conservation with the highest contribution to social welfare, can be formulated as a biobjective optimization problem: x∈{1,2,3,4} n ∈ s.t. land transition rule s minimum/maximum land cover i∈I x∈{1,2,3,4} n max max x∈{1,2,3,4} n max ℝ Ecology and Society 26(1): 9 https://www.ecologyandsociety.org/vol26/iss1/art9/ The landscape is divided into multiple decision units I = {1, ..., n} that determine the number of decision variables of the optimization problem and x is a vector of indicator variables x i stating the land use of patch i ϵ I. The decision variables take only discrete values that represent the land-use type of the respective unit, i.e., 1 for cropland, 2 for extensive grassland, 3 for intensive grassland, and 4 for forest. Therefore, Eq. 5 can be characterized as a combinatorial optimization problem. The first line of the problem formulation maximizes the sum of the contribution margins c for all patches x i, and the second line does the same for the WTP for biodiversity b. Furthermore, the optimization problem is constrained by specific land transition rules (i.e., which land use can be converted into what other land use) and values for the minimum and maximum land cover of each land-use type (for further details, see Methods: Application).
This type of problem is usually solved by applying local search algorithms such as genetic algorithms (Kaim et al. 2018). We use the optimization tool CoMOLA (constrained multi-objective optimization of land allocation; Strauch et al. 2019), which is based on the nondominated sorting genetic algorithm II (NSGA-II; Deb et al. 2002). For each optimization run, CoMOLA first creates an initial set of land-use maps (individuals). These solutions are passed to the biodiversity model and SWAT lookup table that return the respective WTP. Then, the NSGA-II evaluates the WTP of each solution, and those with the best performances (i.e., highest WTP for each objective) are merged by crossover and mutation, and thus, form the next generation. Individuals of the offspring population that do not satisfy the constraints are considered infeasible and are transformed into feasible individuals by a repair mutation algorithm . Again, the WTP for these solutions are calculated, evaluated, and used as parents for the next generation. This process continues until a certain stopping criterion, e.g., the maximum number of generations, is met. The algorithm then selects those solutions that are nondominated. In a maximization case, a feasible solution is nondominated or Pareto-optimal if there is no other feasible solution that performs better in at least one objective without decreasing another objective simultaneously (Coello Coello et al. 2007). Finally, the set of all nondominated solutions, i.e., land-use strategies, forms the Pareto frontier.
We select the solution with the highest contribution to social welfare from the Pareto frontier by maximizing the sum of each objective's WTP (see simplified illustration in Fig. 1 red box and Eq. 6).
Here, P is the set of nondominated solutions p, B the WTP for biodiversity, and C the WTP for agricultural production.

Application
We applied our approach to the Lossa River basin, a sub-basin of the TALE project's (https://www.ufz.de/tale/) case study area, the Middle Mulde River basin, which is located in central Germany (Fig. 3). The case study area was selected because it is a river basin (as required for running SWAT) and its size (14,076 ha) is not too large for running the optimization via CoMOLA. It is dominated by agricultural land use, primarily winter wheat (36.1%), rapeseed (22.6%), winter barley (15.6%), and corn (14 %). Also, forests and mainly intensively used grassland cover parts of the landscape (Figs. 3 and 4A). Furthermore, the area is habitat of different open-land bird species, of which we take into account nine Red List species in the biodiversity model.  The study area was divided into HRUs, which served as decision units for the optimization but were also needed to apply the SWAT model, which calculates the contribution margins for each HRU. For cropland, we used typical crop rotations representative of the status quo. As shown in Table 3, We distinguished four land-use types: cropland, extensive grassland, intensive grassland, and forest (Table 3). The land-use transition rules and maximum allowed land cover of each land-use type (Table 3) were used to constrain the optimization (Eq. 5) and were set according to the results of a stakeholder workshop with farmers, conservationists, and employees of state ministries. The procedure and results of that workshop are described in Karner et al. (2019). Parts of the case study area that were urban, wetlands, water bodies, or forests in the status quo were not allowed to change and were thus excluded from the optimization process. Similarly, to decrease the number of decision variables, HRUs < 5 ha were not considered in the optimization, though the associated contribution margins were added to the total value of agricultural production. Approximately 9657 ha (~69%) of the case study area was considered in the optimization. To avoid extreme land-use changes due to very large HRUs, we split them along roads and railways so that no decision unit was > 1 km², leading to a total of 243 HRUs. For each HRU, the optimization algorithm varied the land use according to the transition matrix (Table 3). Furthermore, the optimization was constrained by minimum fitness values for the two objectives, i.e., solutions with a value < 0 €/ha were excluded in both cases. With these settings, we ran the genetic algorithm with a population size of 100 individuals (i.e., land-use maps) over 200 generations (each generated by crossover and mutation from the previous generation). This procedure was done in seven independent repetitions (i.e., optimization runs), due to the stochastic nature of the algorithm, for each of the two WTP functions.

RESULTS
The overall Pareto frontiers for both biodiversity WTP functions illustrate that, with increasing biodiversity, the solutions of WTP 1 first dominate those of WTP 2 and then, after an inflection point in the middle of the plot, the solutions based on WTP 2 tend to dominate those based on WTP 1 (Fig. 5). This pattern could be explained by the shape of the WTP functions (Fig. 2), where WTP 1 first dominates WTP 2, and then after their intersection at a biodiversity index of 1.3 (88.57 €), WTP 2 dominates WTP 1. The plot also shows that the current land use can be improved by all solutions found by the optimization algorithm. For example, the solutions located at the edges of the Pareto frontiers are those with the highest WTP for agricultural production (Agrimax) or biodiversity (Biomax). Compared to the status quo, the Agrimax solutions indicate a possible gain in contribution margin by up to 10 €/ha, whereas the WTP for biodiversity only increases by about 8 €/ha (WTP 2) and 19 €/ha (WTP 1; Fig. 6). In contrast, the highest possible gain in WTP for biodiversity (Biomax) is, according to the definition of the WTP functions (Eqs. 3 and 4), 213.36 €/ha. In this case, the gain in contribution margin would only be approximately 3 €/ha for both WTP functions.

Fig. 5.
Pareto frontiers based on the two willingness to pay (WTP) functions, WTP 1 and WTP 2. Each of the seven optimization runs per WTP function returned a set of Paretooptimal solutions (light-colored, unfilled points). The nondominated solutions of these sets form the overall Pareto frontiers. Orange asterisks at the upper left indicate points with the highest WTP for agricultural production (WTP 1, light orange; WTP 2, dark orange) and green asterisks at the lower right indicate points with the highest WTP for biodiversity (WTP 1, light green; WTP 2, dark green). Blue asterisk = status quo. Fig. 6. Changes in willingness to pay (WTP; € ha −1 yr −1 ) and land use (%) compared to the status quo for the Pareto-optimal solutions with maximum WTP for agricultural production (Agrimax) and maximum WTP for biodiversity (Biomax) based on WTP functions 1 and 2, respectively.
By comparing WTP 1 and WTP 2, there is no major difference between the Agrimax or Biomax solutions in the land-use map. Generally, for both extreme solutions, cropland decreases in favor of intensive grassland (Fig. 6), which nearly reaches the maximum https://www.ecologyandsociety.org/vol26/iss1/art9/ land cover allowed by the constraints. Nevertheless, there is a higher loss in cropland for Biomax than for Agrimax, which is due to a gain in extensive grassland. In all cases, there is no significant change in forested area. Because the results for the two WTP functions are very similar, we next focus only on those based on WTP 2.
According to Eq. 6, we identified the Pareto-optimal solution with the highest contribution to social welfare (yellow star in Fig.  7; Fig. 4B) to be the solution with the highest WTP for biodiversity (Biomax), which has already been discussed. Instead of comparing spatial land-use changes for this solution only with the current (status quo) landscape, we decided to analyze the 15 best solutions (yellow to orange in Fig. 6). This procedure makes it possible to identify not only patches that have to change to achieve a particular solution, but also how frequently patches have been selected for a particular land use when the solution reached a high contribution to social welfare.

Fig. 7.
Overall Pareto frontier of the optimization based on the willingness to pay (WTP) function 2. Light grey unfilled points = the sets of nondominated solutions for all seven optimization runs. Colored points and scale = the social welfare of the overall nondominated solutions starting from 250.85 € ha−1 yr−1 (status quo, dark blue) to 466.84 € ha−1 yr−1 (best solution, yellow). Blue asterisk = current state of the landscape; yellow asterisk = optimal solution with highest social welfare.
As already noted from the comparison of land-use changes for the extreme solutions (Fig. 6), there was also no significant change for the 15 best solutions of WTP 2 in forest. Cropland was mainly transformed into intensive grassland because it has a higher contribution margin per area than does extensive grassland and, in many cases, cropland (Figs. 8 and 9). Furthermore, a considerable number of patches with intensive (extensive) grassland were mainly converted into extensive (intensive) grassland, and few of them into forest. Intensive grassland close to water bodies was not converted to another land use (Fig. 8). On the contrary, there was a significant gain in grassland, which occurred because the biodiversity model uses a distance parameter for water bodies and thus prefers grassland along these areas. Fig. 8. Gain and loss of intensive grassland patches for the 15 best solutions for willingness to pay function 2 compared to the status quo. Red = loss of intensive grassland; orange = no change; green = gain in intensive grassland. Counts indicate the number of times out of 15 that the patch was affected. Fig. 9. Gain and loss of extensive grassland patches compared to the status quo and prediction of suitable habitat for the Northern Wheatear across the 15 best solutions for willingness to pay function 2. Red = loss of extensive grassland; orange = no change; green = gain in extensive grassland. Purple dots = suitable habitat for Northern Wheatear. Counts associated with colors indicate the number of times out of 15 that the patch was affected (grassland) or how often the bird was predicted at that spot (Northern Wheatear).

DISCUSSION
The results of the optimization suggest land-use strategies that improve the social welfare in the Lossa River basin, as measured in terms of contribution margins from agricultural production and WTP for increases in biodiversity. Particularly, the social welfare gains from biodiversity can be substantial without compromising on economic gains from agricultural production. This analysis shows that optimization algorithms can help to identify desirable land-use options as compared to the status quo. Ecology and Society 26(1): 9 https://www.ecologyandsociety.org/vol26/iss1/art9/ They can do so not only by changing one objective in favor of another; in some cases, even win-win solutions may be attainable (Hanspach et al. 2017).
The main aim of our exploratory study was to test and implement the conceptual idea of combining multiobjective land-use optimization with preference information in a real-world context. With this in mind, our specific results should be regarded with caution and not be interpreted in terms of policy recommendations. Rather, they are stylized results showing a "proof of concept". We next discuss a number of challenges to provide a ground for future research in the field of socialecological landscape optimization.

Quality of underlying data and models
To a large extent, the optimization results are influenced by the quality of the underlying data and models. The bird model that we used considers nine Red List bird species. Their habitat requirements are relatively balanced among the different land-use types, though with a focus on open-land bird species because of the characteristics of the study area. We assumed that the selected species serve as a good indicator for the overall biodiversity level in the region, which implies that if there is suitable habitat for these threatened species, more common species will also occur in these habitats. However, we are aware of the possible bias of this approach and that there might be species that cannot be represented by this indicator (Bonn et al. 2002). However, similar bird-based indicators of biodiversity are quite common in the literature and are, for example, used in the German national strategy on biological diversity (Federal Ministry for the Environment, Nature Conservation and Nuclear Safety 2007). Furthermore, the bird model from Jungandreas can be applied to any bird species for which the necessary data are available so that future studies could select other species they consider representative for the study area.
The biodiversity model can also help to explain the optimization algorithm's allocation of intensive grassland close to water bodies. The combination of urban areas that include abandoned buildings, parks, and small gardens with cropland and, particularly, grassland, offers habitats for nesting and feeding for many bird species considered in our study (e.g., Western Jackdaw, Common Redstart, Crested Lark, and Barn Owl). In the Lossa River basin, many towns and villages are located close to water bodies, which might cause the allocation of grassland close to these areas. Some birds also prefer the proximity of water (e.g., Northern Lapwing, Kingfisher). Additionally, the grassland along Lossa Creek connects the different habitats, purely due to the geographical characteristics of the case study area because the bird model does not consider habitat connectivity. This idea should be taken into account when applying the model to other areas.
As could be expected, forest cover did not change significantly. This result is because of both land-use transition rules that prohibit the cutting of forest and the low potential for afforestation (maximum land-cover rules) in the mainly agricultural area, which is why we focused on open-land bird species.
Extensive grassland patches increase, especially in the southeast of the study area (Fig. 9). Contrary to expectation, this increase is not because of low soil quality and thus low contribution margins for agricultural production. Rather, the biodiversity model predicts this area as a highly suitable habitat for Northern Wheatear (Oenanthe oenanthe). The gain in large extensive grassland patches coincides with the prediction of this bird (i.e., dark and light purple dots in the dark and light green patches, respectively, in Fig. 9). The prediction of Northern Wheatear in the north and central south relate to the suitable habitats the bird already had before the optimization.
The high values for the bird indicator and thus the resulting high WTP for biodiversity for these solutions can be explained by the immense gain in suitable habitat for Northern Wheatear. Breaking down the bird indicator, which is about 1.6 for all best solutions, into the indicator values of each bird species, Northern Wheatear achieves an average index of 5.61. The indicators of all other species remained near the status quo (see Table A2.1 in Appendix 2).
Similarly, the simulation of biomass production with SWAT is based on simplified assumptions about farm management operations (e.g., crop rotations, fertilizer, and tillage practices). A different and possibly more detailed model may influence the results but would also require more detailed information input, which is usually not available at the landscape scale.

Detail of modeling vs. computational effort
If data availability allows, the models and the formulation of the optimization problem could be improved by adding more detail. For example, one of the outcomes of optimizing agricultural production and biodiversity is that additional patches of intensive grassland are located mainly close to water bodies, although that may negatively affect water quality. Considering other objectives such as maximizing water quality could thus lead to different optimization results. Also, the analysis is only partial in economic terms. It is based on the implicit assumption of ceteris paribus ("other things being equal"); ideally, a larger bundle of ecosystem services (and possibly other effect domains) would be considered in the optimization. However, such inclusion could be challenging given that the approach requires that all objectives can be described economically using the same units (i.e., € ha −1 yr −1 ). Furthermore, we used the CoMOLA optimization framework that can be coupled with various spatially explicit models. An economic model could be added instead of or additional to the biodiversity and SWAT model. However, it should always be taken into account that adding models, objective functions, or decision variables to the optimization can increase the computational effort substantially.

Quality of preference information
The shape of the WTP function can influence optimization results (Fig. 5). We tested two possible versions, a simple piecewise linear function and an economically more realistic 6th degree function ( Eqs. 3 and 4). Of course, other functional relationships between the three given points (Table 2) could be considered. In this context, it should also be noted that the saturation point has an influence on the point with highest social welfare from biodiversity and in our case also on the optimal solution because it sets a limitation for the maximum bird index that can be achieved in the landscape. For the optimization algorithm, solutions with a bird index of, e.g., 1.8 are equal to solutions with an index of https://www.ecologyandsociety.org/vol26/iss1/art9/ 1.6 because they coincide with the same WTP. Therefore, there is no incentive for the algorithm to find solutions with an even higher biodiversity level; instead, it would aim to improve the second objective, the agricultural contribution margin.
The quality of the preference information used in an optimization exercise such as ours plays an important role. In this respect, we used rather simple proxies. For policy relevance, regionally specific contribution margins and WTP for biodiversity would have been necessary; however, these were neither available nor necessary for our exploratory purposes. In this context, it should also be mentioned that the simplified assumptions made in Eq. 2 (population and agricultural land are evenly distributed in Germany) could be improved. Ideally, the calculation of WTP per unit area would consider information about the share of agricultural land and the number of people living within a 15 km radius of each respondent in the choice experiment. Because we did not have any information about the distribution of the survey participants and their WTP beyond this 15 km radius, our calculation in Eq. 2 is imprecise and only an approximation. These challenges should be taken into account in the assessment of WTP that might later be used on a different spatial scale.
Furthermore, it may be objected that monetary expressions of preferences are imperfect; their informational value is limited, especially in the context of biodiversity (Pascual et al. 2017). They provide only a preliminary indication of the contribution of ecosystem services and biodiversity to human well-being and should be complemented by other indicators in the decisionmaking process (Förster et al. 2019).
In our case, even if there were different values for WTP for biodiversity, the optimization results would probably not change significantly because of the Pareto-based two-objective optimization approach. The algorithm searches for solutions that maximize both WTP for agricultural production and WTP for biodiversity. In other words, changes in the scale of one objective would not change the topology of Pareto-optimal solutions. In the case of a weighted-sum approach, in contrast, the results would be completely different because these algorithms consider only a single-objective function that already aggregates both WTP values. For further discussion of Eq. 2 and a sensitivity analysis of the effect of WTP values on the results, we refer to Appendix 3.

Identification of winners and losers
A related issue is the implicit application of the Kaldor-Hicks criterion (Johansson 1993) that goes along with the type of preference information used in our study. It does not allow for the analysis of winners and losers. Generally, there are two possible ways to add such an analysis: ex ante or ex post. The exante variant was discussed by Cavender-  and implies the identification of differently affected groups (e.g., farmers vs. the general public) before preference data are collected and then uses them to see how much the losers lose. If preference data are collected directly in the study area, it may also be possible to conduct cluster analysis on the basis of demographic and socioeconomic control variables. In this case, however, it would be necessary to have objective functions that maximize WTP for biodiversity for each stakeholder type to distinguish among stakeholders in the Pareto frontier. This process could easily lead to a many-objective optimization problem (i.e., more than four objectives) with an increased complexity that has to be solved by algorithms other than NSGA-II. The alternative is to calculate only biophysically optimal land-use configurations and to include preference information in an ex-post evaluation, preferably involving a wide range of stakeholders.

CONCLUSION
We presented a method that combines mathematical optimization with methods from ecology and economics to identify socially optimal land-use strategies that take into account the biophysical potential of a landscape. We applied the method to a case study in the Lossa River basin in central Germany. The optimization could identify landscape configurations that would increase overall social welfare in the region by creating a win-win situation between the maximization of social benefits from biodiversity and agricultural production. The integration of preference information into the biophysical optimization allows reducing the usually large set of Pareto-optimal solutions to a single (or manageable number of) societally most beneficial solution(s) without having to apply subsequent utility functions  or methods such as weighting approaches (see, for example, Qi and Altinakar 2011). This selection of solutions can also be used as input for further stakeholder-based analyses and decisions. Though our approach does not allow for a purely biophysical trade-off analysis, a trade-off analysis with respect to the different contributions to social welfare of each objective is still possible. The method can help to answer questions of sustainable land-use planning by providing spatially explicit solutions. These solutions can also be used for the development of site-specific policy instruments, e.g., spatially explicit agrienvironmental schemes or management standards, which will help to balance ecological and social demands on landscapes.
Our method can be applied to other real-world case studies, even in broader or different contexts. The selection of the objectives is flexible as long as they can be linked to preference information measured in the same (e.g., monetary) unit so that calculating changes in overall social welfare is possible. For future studies that follow the optimization of biodiversity and agricultural production on a regional or local scale, we suggest paying particular attention to the quality of the data and models used because they have a strong influence on the final result. Another general issue that should be borne in mind when interpreting the results of multiobjective social-ecological optimization (and an area where we see the largest need for further innovative research) is that it is a purely "comparative-static analysis", i.e., it compares one state (X) to another (Y). It does not tell much about the road from X to Y, regarding practical feasibility given political, behavioral, economic, and other constraints. To study these contraints and how the socially optimal state Y may actually be achieved, one would need additional information and methods such as agent-based models (Brunner et al. 2016, Haslauer et al. 2016 -Hydrological data o Measured discharge (daily) and water quality (e.g. N and P species, bi-weekly to monthly) for different gauging stations and varying time periods o Point sources: average nutrient loads from waste water treatment plants

A1.2 -Calibration and performance of the SWAT model
For the Middle Mulde Basin, SWAT was calibrated against daily streamflow and monthly loads of sediment as well as nitrogen and phosphorus fractions at three gauges (Bad Düben, Golzern and Erlln), taking into account four performance metrics (Nash-Sutcliff-Efficiency -NSE, relative index of agreement -rd, Kling-Gupta-Efficiency -KGE and percentage bias -PBIAS). We used the R package 'SWATpasteR' (Schürz et al. 2017) for both sensitivity analysis and model calibration. After sensitivity testing of 130 model parameters which were sampled for 14,076 model simulations using the STAR method (Razavi and Gupta 2016), we carried out further 10,000 Latin-Hypercube sampled SWAT simulations to calibrate the most sensitive 29 model parameters. Table A1 summarizes the performance values of the best model simulation, which we defined as the simulation with the minimum total rank sum of all performance metrics for all variables and gauges.
According to the classification of Moriasi et al. (2015), model performance ranged from satisfactory to very good for daily streamflow and was very good at all gauges for monthly loads of Nitrate-N. For all other variables, model performance was rather unsatisfactory with a few exceptions at single gauges (Tab. A1.1, see also Fig. A1.1).
Moreover, we manually adjusted SWAT plant parameters in order to fit simulated yields of single crop types to the yields observed in period 1995 to 2009 (Fig. A1.1). Our results illustrate the difficulty in identifying a single parameter set that simultaneously satisfies multiple performance criteria for multiple variables at multiple gauges in the Middle Mulde Basin.
However, it nevertheless appears reasonable to use the model for land use change impact assessments regarding crop yield, streamflow, Nitrate-N and sediment loads. Figure A1.1: SWAT model performance for predicting daily streamflow, monthly water quality and mean annual crop yields the Middle Mulde Basin.