Modeling how land use legacy affects the provision of ecosystem services in Mediterranean southern Spain

Land use decisions induce legacies that affect the welfare of future generations. Here, we present a spatial modeling approach for quantifying how past land use decisions influence provision of multiple ecosystem services (ESs) based on different land use trajectories. We modeled the effect of past land use changes on water regulation, soil protection and habitat quality in southern Spain, one of the most transformed areas of the Mediterranean region. We demonstrate a measurable influence of antecedent land use changes on the capacity of a given land use to provide ESs, and that the effect size can vary among different services and land use trajectories. Our results suggest that afforestation programs may decrease habitat quality but not alter soil protection, depending on whether the previous land use was cropland or shrubland. Although it is well-established that land use legacies motivated by past land decisions are ubiquitous and crucially important for effective landscape management, the question of how the magnitude and spatial distribution of ES supply vary under different land use trajectories remains unknown. Our approach enables quantification of how land use legacy affects ecological processes that underpin ES capacities at a regional scale, which will allow land managers to develop more accurate landscape planning strategies for preserving ESs.

Introduction Increasingly, regional land use decisions such as the implementation of restoration programs or declaration of protected areas are made based on the biophysical assessment of landscape-scale ecosystem services (ES) [1]. For that, managers and decision makers have access to a variety of new tools for mapping ES to identify areas of the landscape that have the capacity to provide simultaneously multiple ES or 'bundles' [2]. These tools enable land managers to test ES provision under different scenarios, i.e. different configurations of land use that result with different land policies [3]. ES provision is calculated for each scenario, and comparison among different scenarios enables managers to identify which land use decision will preserve future supply of multiple ES. While the ES bundles and scenarios approach has proven to be a powerful tool to enable better regional-scale land use decisions [4], so far, these approaches ignore the critical role that land use legacies play in understanding ES provision, i.e. effect from prior land use that are still propagating through the ecosystem [5].
It is well-established that land use legacies are ubiquitous, and crucially important for effective landscape management because they affect ecological processes underpinning ES supply [6][7][8][9][10]. However, most landscape-scale assessments of ES are based on the relationship between the spatial patterns of ES and current attributes of land uses [11]. It has not yet been empirically tested how multiple ES are influenced by past land use history (e.g. 10), and unraveling the effects of prior land use change on current ES provision would enable more accurate landscape planning strategies for preserving future ES supply [1]. Recent studies have made innovative progress on legacy knowledge gaps. For instance, Locatelli et al [12] reviewed existing literature and introduced the concept of land use trajectories as a mean of 'pathways of land change' that influence ES over time for mountain systems. Martin et al [5] developed a novel method to measure land use legacy for a single ES (i.e. water quality) in lake ecosystems. However, recent literature has not yet addressed how different land use trajectories may influence multiple ES, nor have they introduced approaches that can be applied to diverse ecosystem types at a regional scale.
This study presents a spatial modeling approach for empirically evaluating how diverse land use legacies affect multiple ES supply at a regional scale. We present this approach as transformative, in that it can be integrated into standard ES modeling approaches (such as Integrated Valuation of Ecosystem Services and tradeoffs (InVEST) and others) and as an advance tool that land managers can integrate in decisionmaking. We conducted our study in Southeastern Spain (figure 1(a)), where land use legacies are particularly relevant because the region has experienced massive land use transformations after the 80s [13], and because there is active landscape restoration planning underway to preserve future ES supply [14]. To demonstrate our approach, we: (1) quantified and mapped the provision of three regulating services (i.e. water regulation, soil protection and habitat quality) based on current land use; (2) mapped the five main land use trajectories that occurred over the last 50 years (figure 1(b) and table 1), and (3) modeled how these land trajectories have affected current ES provision. Finally, we discuss the implications of land use legacies underpinning changes in ES, and conclude with potential applications for land management and restoration programs.

Methods and materials
Study area The Arid Southeast Spain (figure 1(a)) has experienced since 1956 one of the most significant land use change transformations in all Europe [14]. This area covers approximately 1220 000 ha, and comprises highbiodiversity, ecologically vulnerable Mediterranean arid ecosystems, and land use changes altering their capacity to provide ES [15]. In the last 60 years, landplanning strategies to promote economic development have motivated three major land use changes: (1) a transition from traditional agriculture toward intensive greenhouse horticulture; (2) rural abandonment as rural people migrate to urban areas; and (3) the implementation of a protected natural areas network [14]. As a result, this region has high diversity of land uses, in which cropland (e.g. almond-trees or olive groves), shrubland, and forest (mainly reforestation of pines) are dominant (with 43.15%, 38.0%, and 11.97%, respectively). Greenhouse horticulture (3.77%), watercourses (1.23%), grassland (1.13%), urban (0.46%), and bare soil (0.29%) cover the remaining landscape ( figure 1(a)).

Modeling approach
Martin et al [5] defined legacy effects as those effects from prior human disturbances that are still propagating through the ecosystem. In particular, historical human-induced land use changes may result in underpinning legacy processes that influence current ecosystem functioning and structure, biodiversity and ES. Thus, the modeling consisted in first exploring the current capacity of different land uses to provide ES, and then exploring how land use trajectories affect ES provision. Specifically, our modeling approach was based on the three principles shown in figure 2.
ES and land use trajectories mapping ES mapping techniques included APLIS model for water regulation [16,17], the universal Soil Loss equation (USLE) model for soil protection [17], and the InVEST model for habitat quality [18]. Resulting ES maps were obtained in raster format with a resolution of 100 m (figure 3 and supplementary material is available online at stacks.iop.org/ERL/13/ 114008/mmedia-SM 1 and SM 2 for data procedure). Current and past land use types were extracted from the land use vector map of Andalusian region for the year 2007 and 1956, respectively (Environmental Information Network of Andalusia, www. juntadeandalucia.es/medioambiente/site/rediam). We generalized on eight land use types based on the International Geosphere-Biosphere Programme land classification (IGBP), as follows: bare soil, cropland, forest (mostly evergreen needle-leaf forest), grassland, shrubland, watercourse, urban, and greenhouses ( figure 1(a)). Although greenhouses do not belong to the general IGBP classes, we included them in our models because it is a very common intensive agricultural practice in some parts of our study area. We employed those eight land use types to model water regulation, and the same except urban and greenhouses to model soil protection and habitat quality because their capacity to provide these ES is considered as null [18]. To map the five most prevalent land use trajectories in the study area from 1956-2007 [19] we used tranUSE, a free software to interpret land use changes based on trajectories defined by the user [20]. These trajectories were: rural abandonment, agricultural intensification, deforestation, afforestation, and no change (table 1 and figure 1(b)). These land use trajectories have been recognized for initiating legacy processes by affecting forest composition, Figure 1. (a) Study area location and spatial pattern of land use types. Since our study case was focused on the arid and semi arid regions of Southeast Spain, we excluded the high mountain areas which did not meet this criterion. (b) Land use trajectories. Our spatial modeling approach quantifies different ES capacities within a single land use type, and demonstrates the role of land legacy. Note. * The-no change-trajectory does not assume that no land use change occurred in-between. We note that the 10.57% of the area labeled as-no change-had at least one land use change between 1956-2007. This area covered 0.00073% of the whole study area.  vegetation pattern, soil structure, etc [6]. As an example, forests reverting from agriculture have been shown to have legacy effects on processes such as soil nutrient dynamics and biodiversity [21,22]. Deforestation has long-term effects on N content in soils [23] (table 2). Finally, we rasterized land use and land use trajectory maps to a 100×100 m pixel size to extract the predictor variables used in models.
In summary, the three ES mapped (i.e. water regulation (APPLIS model), soil protection (USLE model) and habitat quality (InVEST model)) were used as response variables in the LU-models and LUxT-models (see 'Modeling of ESs and land use legacy' subsection). Likewise, both land use type and land use trajectory were used as predictors.

Modeling of ESs and land use legacy
Mixed-effect models were built (package lme4 and function lmer in R, www.R-project.org) to estimate (1) the current level of ES provision across land use types (hereafter LU-models), and (2) the influence of land use trajectories on the level of ES provision of current land use (hereafter LUxT-models). We modeled three key ES: water regulation and soil protection, gamma distributed with log as link function, and habitat quality, logit transformed and normally distributed with identity as link function (see supplementary material SM 2). We were interested in making inferences about the mean of current land use, compared to the whole of the study area in terms of ES provision rather than in testing differences between particular land use types. For that, LU-models included varying-intercept and land use as random effect. Similarly, LUxT-models included varyingintercept, but they also incorporated the statistical interaction between land use and land use trajectory as random effect (see LUxT-Models below). In addition, we tested the significance of land use trajectory effect on ES provision across current land use by comparing LU-models and LUxT-models (both estimated by restricted maximum likelihood) in terms of deviance explained by performing a likelihood-ratio test.

LUxT-models
The goal of these models was to explore how land use trajectories may modify the current capacities of land use types to provide ES, which were inferred previously by the LU-models. In these models, the mean of ES provision by current land use combined with the land use trajectories was compared to the ES mean of the whole study area. The model equation was:

Results
LU-models (i.e. models that included only current land use as a predictor) showed variation in the effects of land use on ES (table 3). Forest reached the highest positive effect for the three regulating services, while greenhouses, bare soil, and watercourse showed negative effects on the ES supply. Among all land use effects, cropland showed a significant positive effect for water regulation (effect size=0.34) and habitat quality (effect size=0.65), but showed a negative effect for soil protection (effect size=−0.28). Both grassland and shrubland showed a strong negative effect on habitat quality (effect sizes=−0.43 and −0.78, respectively) and a positive effect for water regulation (effect sizes=0.15 and 0.10, respectively) and soil protection (effect sizes=0.18 and 0.51, respectively).
By incorporating the land use trajectories in the models, we found variation in the effects on ES provision with respect to the specific land use (table 4). For instance, the three trajectories leading to current forest: (a) forest to forest (i.e. no change); (b) agriculture to forest (i.e. rural abandonment); and (c) shrubland to forest (i.e. afforestation) manifested in different capacities of ES supply ( figure 4). Specifically, the provision of water regulation and habitat quality varied among the three trajectories, but the soil protection capacity of current forest cover remained consistent regardless of past land use. The ES provision of cropland also differed depending on past land use, in particular for water regulation and habitat quality. For example, under the agricultural intensification trajectory, the effects of both ES moved from positive to negative (effect sizes=−0.04 and −0.89, respectively). The variation in the effects on ES provision by shrubland was the highest. For instance, water regulation and soil protection were positively affected under the deforestation trajectory (table 4).
Overall, the deviance explained by LUxT-models was significantly higher than the deviance by LU-models across all ES provision. Please see table S1 in supplementary material SM 3 for more details.

Discussion
Measuring the capacity of different land use types to simultaneously provide multiple ES is crucial to understanding the trade-offs and synergies associated with land management decisions [3,11]. While research has been conducted to model the ability of different past and current land uses to provide ES (see for example, 7, 5), our analysis here is the first modeling the effect of land use trajectories on multiple ES concurrently, and provides a transformative approach to incorporate potential effects of land use legacy on spatially-explicit ES assessments over broad spatial scales.
Our results demonstrate a measurable influence of antecedent land use changes on the current capacity of land use to provide ES. In addition, we measure the degree to which this effect varies among different ES. For example, it is well-established that forests are one of the most important land cover types in terms of ES provision [38]. Our results confirm this to be the case in our study area for the ES that we measure here: water regulation, soil protection, and habitat quality [39]. In the first step of our modeling procedure we quantified the provision capacity of forest compared to other land use types in our study area, and in fact, forests had the highest rates of all three ES provision.
In the next step, we incorporated the land use trajectories of forested pixels, and our results showed that the current forest capacity to provide habitat quality also depends on such trajectories. In our study area, the afforestation trajectory represents pine plantations that were established for the purpose of recovering areas affected by intense mining activity in the 19th century and rural abandonment in the middle of the 20th century. Results indicated that those plantations provide much less habitat quality compared to oldgrowth forests (e.g. pine forest that have not undergone change, i.e. no-change trajectory), but both trajectories were equally effective at soil protection. These differentiated patterns of ES among land patches with the same current land use but that come from different land use trajectories are likely motivated by legacy processes that still continue to affect ecosystems and the ES they provide at present (table 2). Indeed, afforestation and the homogenization of tree species composition at a regional scale have been recognized for initiating land use legacies on ecosystems function (by altering spatial-temporal dynamics of ecosystem productivity), structure (availability of habitat elements, for example, stand structure in forests), and biodiversity (changes of species composition) [6,22]. Our findings are consistent with case-studies which demonstrate the important role of natural forests in providing water regulation, soil protection and habitat quality [38,39] compared to pine plantations [40][41][42]. Thus, our modeling approach has important implications for the assessment of the restoration programs derived from the UE Rural Development Policy. This policy aims to restore and preserve ecosystems related to agriculture and forestry which were affected by past land use decisions [14,43]. We found that the capacity of a land unit to provide habitat increased with the deforestation-to-cropland trajectory. Mediterranean farmlands can result in beneficial environments for generalist wildlife species that can exploit the new food resources available in human-dominated landscapes and thus reach higher occurrences than in more natural areas [44][45][46]. These rural agricultural environments can be particularly favorable in lowlands of the arid Southeast Spain, which have low-diversity forests and therefore fewer resources are available for wildlife compared to heterogeneous semi-natural habitats [47]. Deforestation can result in a greater spatial heterogeneity of land cover types and hence, better habitat quality at a landscape scale [45]. However, it is also important to highlight the importance of maintaining at least some forest in this landscape, and specially scattered forest fragments, to support biodiversity at board spatial scales [45,48]. ES assessments based on regional land use scenarios are commonly incorporated into decision-making, but they do not consider the effect of land use legacy. The question of how the magnitude and spatial distribution of ES supply vary under different land use trajectories is one of the key knowledge gaps in ES science. Our approach measures different ES capacities within a given land-use type (e.g. forest), and links these within-type differences to land-use legacy. Many land use maps include only a single forest class, but multi-temporal land use maps showing forest/nonforest are becoming commonly available. Past land use maps can be used to define trajectories that serve as a proxy for different forest types according to our modeling approach. Thus, more reliable ES maps can be provided that will enable decision-makers to more accurately incorporate natural capital and ES into policy and management [1,43]. Modeling approaches such as proposed here are valuable to anticipate the regional-scale impacts of current land use decisions on future ES supply [49]. Future research should test the accuracy of the proposed approach for different ES categories and in diverse study systems.

Acknowledgments
This work is supported by the European LIFE Project ADAPTAMED LIFE14 CCA/ES/000612, the NSF Idaho EPSCoR Program and by the National Science Foundation under award number IIA-1301792 and the GLOCHARID Project. The study was conducted in the Arid Iberian South East LTSER Platform-Spain (LTER_EU_ES_027), https://data.lter-europe. net/deims/site/lter_eu_es_027. The models were run by using the High Performance Computing cluster of the Andalusian Computing Center, https://cica.es/ servicios/supercomputacion/. The research reported in this paper contributes to the Programme on Ecosystem Change and Society (www.pecsscience.org).

ORCID iDs
Juan Miguel Requena-Mullor https:/ /orcid.org/ 0000-0002-5120-7947 Model results are on a log scale for water regulation and soil protection, and on a logit scale for habitat quality. LU-models: models that included only land use variables; LUxTmodels: models that included both land use variables and land use trajectories.