Assessing the effects of different land-use/land-cover input datasets on modelling and mapping terrestrial ecosystem services-Case study Terceira Island (Azores, Portugal)

Modelling ecosystem services (ES) has become a new standard for the quantification and assessment of various ES. Multiple ES model applications are available that spatially estimate ES supply on the basis of land-use/land-cover (LULC) input data. This paper assesses how different input LULC datasets affect the modelling and mapping of ES supply for a case study on Terceira Island, the Azores (Portugal), namely: (1) the EU-wide CORINE LULC, (2) the Azores Region official LULC map (COS.A 2018) and (3) a remote sensing-based LULC and vegetation map of Terceira Island using Sentinel-2 satellite imagery. The InVEST model suite was applied, modelling altogether six ES (Recreation/ Visitation, Pollination, Carbon Storage, Nutrient Delivery Ratio, Sediment Delivery Ratio and Seasonal Water Yield). Model outcomes of the three LULC datasets were compared in terms of similarity, performance and applicability for the user. For some InVEST modules, ‡ § | ¶,# ‡,¤ © Sieber I et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. such as Pollination and Recreation, the differences in the LULC datasets had limited influence on the model results. For InVEST modules, based on more complex calculations and processes, such as Nutrient Delivery Ratio, the output ES maps showed a skewed distribution of ES supply. Yet, model results showed significant differences for differences in all modules and all LULCs. Understanding how differences arise between the LULC input datasets and the respective effect on model results is imperative when computing model-based ES maps. The choice for selecting appropriate LULC data should depend on: 1) the research or policy/decision-making question guiding the modelling study, 2) the ecosystems to be mapped, but also on 3) the spatial resolution of the mapping and 4) data availability at the local level. Communication and transparency on model input data are needed, especially if ES maps are used for supporting land use planning and decisionmaking.


Introduction
Modelling ecosystem services (ES) allows us to predict the spatial distribution of different ES that sustain and support human life. Moreover, modelling ES enables us to assess changes in spatio-temporal distribution and the state of ecosystems and how they affect the flows of ES to people. Therefore, modelling ES has become an essential tool for mapping and assessing ES, which is heavily used, for instance, in the context of the European Union's (EU) Initiative Mapping and Assessing Ecosystems and their Services (MAES* ) which supports the implementation of the EU's Biodiversity Strategies 2020 and 2030.
The number of readily available model suites has grown in the last decades, with different options to map the actual or potential supply, demand or use of ES. In 2013, more than 17 decision-support tools for ES quantification and valuation models were identified for the assessment of ES (Bagstad et al. 2013) with a growing tendency (Olosutean 2015). Examples range from simple regression models, mechanistic and stochastic models to complex GIS-based ES models (Martínez-Harms and Balvanera 2012). Amongst the most popular open access models, the Integrated Valuation of Ecosystem Services and Tradeoffs model (InVEST) , ESTIMAP (Zulian et al. 2013), Resource Investment Optimization System (RIOS) or ARIES (Villa et al. 2014) are listed, allowing to model ES in either biophysical terms (e.g. Mg of carbon sequestered) or economic terms (e.g. net present value of that sequestered carbon) (Natural Capital Assessment ).
However, issues related to ES maps and models remain unaddressed, especially those related to uncertainty assessments (Schulp et al. 2014). Numerous studies have compared performances and focused on calibration and validation of models (Cong et al. 2020). 1 Many performance criteria have been established in support of this process and numerous methods to compare maps have been presented (Costanza 1989, Kuhnert et al. 2005, Hagen-Zanker and Martens 2008. Model ensembles have been tested for "goodness of fit" to identify uncertainties (Willcock et al. 2020). For these approaches, map similarity has been in the centre of attention, based on validation by comparing (reference) data from a certain point in time and alignment of model results with reality -assuming that model results are close to reality and, therefore, imply plausibility of the model itself (Hagen-Zanker and Martens 2008).
The InVEST model has been applied in numerous ES mapping applications and first calibration studies highlight the model sensitivities to specific input parameters (Hamel et al. 2015, Redhead et al. 2016, Sharps et al. 2017. Whilst model sensitivity to individual parameters has been well documented, an issue that has often been overlooked is how input land-use/land-cover (LULC) data impact the quality of the model outputs. This study tests the performance of three different LULC geodatasets. These included: 1.
2018 land cover map of the Azores Archipelago (COS.A 2018), the Azores Region official LULC map and 3.
Therefore, the rationale for this ES assessment is two-fold: • it presents an overview of modelled ES for Terceira Island, based on the InVEST model suite; • at the same time, it aims to assess the influence of different LULC input datasets for ES modelling and the respective consequences when using these different input LULC to compute ES maps for policy-and decision-making.

Material and methods
In this chapter, the case study area, the selected models and selected input datasets, as well as the model parameterisation, are described.

Case Study Area
Terceira Island is part of the Azores Archipelago, a European Outermost Region (OR) and an Autonomous Region of Portugal with political and administrative autonomy. The Region is an isolated oceanic archipelago located in the Northern Atlantic, approximately between 37° and 40ºN and 24° and 31ºW. Terceira Island is the third largest island of the archipelago and the third oldest (ca. 3.5 million years) with an area of approximately 400 km (maximum length and width of 21 km and 14 km, respectively). It is part of the Central Group, together with the Islands of Faial, Pico, São Jorge and Graciosa (Fig. 1). Most of the Island surface (72%) has altitudes lower than 400 m and only 1% higher than 800 m, with the highest point of the Island located in Serra de Santa Bárbara (1023 m). Hydrography tends to be irregular, with few prevailing small watercourses with torrential regimes, associated with springs located on the northern hills of the massif Guilherme Moniz -Pico Alto that are prone to dry out during drought periods. Existing lagoons are also relatively small (PGRH-Açores 2015).
Terceira Island preserved some pristine areas at high elevation (Gaspar et al. 2011) and few natural areas at lower elevations have remained, especially in Praia da Vitória Municipality . However, only 6% of the original native forest exists an area of 23 km (Triantis et al. 2010 intensively managed pastures subject to intensive cattle grazing and high inputs of fertilisers ).
Additionally to the balanced representation of Azores-specific ecosystems, Terceria Island provides well-documented land use information (see Chapter 2.1). With three official LULC datasets, Terceira surpasses other islands of the Archipelago and, therefore, presents an ideal case study area for conducting an ES modelling study.
Location of the Macaronesia biogeographic region, the Azores archipelago and Terceira Island (based on ESRI Basemap).
The MAES process is comparably advanced in the Azores Archipelago. Many scientific ES assessments of the archipelago have been published , Madruga et al. 2016, Mendonça 2012, as a recent review of scientific studies showed (Sieber et al. 2018). The outputs of such ES modelling studies have been used to support local land management strategies and plans (Cruz et al. 2011, Geneletti et al. 2020

The InVEST Ecosystem Services model suite
For this study, a modelling approach, based on biophysical models, has been selected (Tier 3; see Grêt-Regamey et al. 2015), applying the InVEST Models. InVEST is a suite of free, open-access, geographic information system (GIS)-based, stand-alone-software models used to map and assess the goods and services provided by nature. InVEST contains several modules and estimates different ES or ES indicators to quantify and map ES, based on spatial, statistical and physical data. InVEST was developed as an ArcGIS Map toolbox in 2008; however, the current version runs as stand-alone software. Besides InVEST's open access availability, its relatively easy implementation is a key advantage. Disadvantages are that, in some cases, the model uses rather simple assumptions, whereas, in other cases, InVEST functions as a black box model. Hence, not all model processing steps can be followed easily by the user. Furthermore, validation of modelling results often proves difficult (Lautenbach et al. 2010, Ochoa and Urbina-Cardona 2017, Polasky et al. 2011 due to a lack of reference data. Furthermore, some InVEST modules (e.g. the pollination module) give index values as results which cannot be measured directly. Therefore, in many cases, only literature-based plausibility checks can be performed instead of a data-based validation.
InVEST is organised in 21 final ES modules and four supporting tools. For this case study, six ES were selected considering the data availability for Terceira Island: recreation, pollination, carbon storage, erosion control, water quality and flow retention (Table 1).
Even though InVEST is stand-alone software, additional GIS software is required to preprocess data, visualise results and perform any further analysis (e.g. data overlays). In this study, data processing was based on ArcGIS ArcMap 10.5. Model post-processing is described in Section 2.5.

Input LULC datasets
Modelling ES is often implemented on the basis of geospatial data, such as LULC data.
The identification of the ecosystems assessed in this study was dependent on the categories of LULC included in each geodataset used (CLC, COS.A 2018 and Sentinel 2based LULC map, all available as shapefiles).  (Table 2). CLC has a minimum mapping unit of 25 ha for areal phenomena and a minimum width of 100 m for linear phenomena (European Environment Agency 2019). Nevertheless, despite its usefulness at European continental and regional level, CLC has been relatively ineffective to address LULC evolution and change issues in European Small Islands/Outermost Regions. This is due to its very broad geographic scale, very large minimum spatial unit (25 ha Minimum Mapping Unit (MMU) and 115 ha average feature size for Terceira) compared to island size and due to its rather inadequate legend, which prevents an efficient characterisation of typical natural and semi-natural LULC units in these territories (Gil et al. 2012 The official Azores Region LULC Map for 2018 (COS.A 2018) was developed by the Azorean Regional Government (Cruz et al. 2007, Hernâni et al. 2018) and covers the entire Azores Archipelago. The COS.A 2018 is structured hierarchically in three levels, from five major LULC classes on the first level to 29 LULC classes on the third level ( Table 2). The LULC nomenclature largely follows the CLC classification, with slight changes. For Terceira, this LULC contains 25 land-cover types (N3), adding aquatic ecosystem types, such as rivers, riparian zones and inland marshes. Initially, the COS.A 2018 was developed using high/very high resolution satellite imagery data classification and an average feature size of 15 ha. The second phase of development consisted of cartographic generalisation and validation tasks in order to address specific spatial planning needs and requirements.
The Sentinel 2-based LULC map of Terceira Island was developed by the Azorean Biodiversity Group* ( Fernández-Palacios 2018). It is based on multispectral satellite imagery with medium/high spatial resolution, derived mostly from available Sentinel-2 data from 2017 and also includes archived Rapideye and Landsat-8 data to fill gaps and increase classification accuracy. It follows its own classification scheme, based on 14 locally-specific ecosystem types and habitats (Table 2), such as Erica, Pittosporum, Acacia , Pinus or Cryptomeria-dominated forest/scrubland patches. Individual features are most detailed in this dataset with an average feature size of 0.15 ha; however, the focus is clearly on habitat types. Urban LULC is reduced to two classes: built-up areas and bare soil.
All three LULC geodatasets vary in purpose LULC classes and MMU. These differences make them suitable for a comparison of the effects that input datasets can have on ES model outcomes. A visualisation of the three LULC with their individual land use classes is shown in Fig. 2. 3 Figure 2.
Overview of the three LULC geodatasets available for Terceira Island.

Model parameterisation
In addition to the LULC data, the InVEST modules require additional input data, for example, geospatial or statistical information. The amount of needed input data strongly depends on the complexity of the modules, more specifically, the model processes. These range from very simplistic modules, such as Visitation, to highly complex, multiparametrical modules, such as Nutrient Delivery or Seasonal Water Yield.
The majority of input data was obtained from a thorough literature review. As locallyspecific information on Terceira Island was often scarce, these data gaps were filled with the best available comparable data. For example, most of the data, such as information on Azores-specific soils, were not available online as geodata. In many available datasets, the Azores were not included or were represented as hardly visible, undistinguishable pixels (e.g. FAO World Soil Database (FAO et al. 2009or Fischer et al. 2008. Therefore, additional soil data were collected from literature and added to the GIS to compile the needed information (Tsui et al. 2013Strohbach andHaase 2012). Table 3 shows a summary of the biophysical and/or statistical input data used for the modelling with InVEST.  Table 3.

Input data per module
Input data required for the different InVEST modules.
Assessing the effects of different land-use/land-cover input datasets on ...

Analysis
In order to assess the relative importance and effects of each input LULC dataset for the outcomes of the model, it was necessary to statistically compare the individual model results. Such an analysis exceeds the scope of existing pairwise map comparison (Hagen-Zanker and Martens 2008, Willcock et al. 2020). InVEST provides all model result maps as either ArcGIS shapefile or grids. Therefore, a fully automated Python script was written, in which data were converted from any InVEST output format to grid format using cell assignment type "cell centre" at an equivalent continuous scale. The resulting raster grids had the same spatial resolution (100 m x 100 m) and were clipped to the Island's extent. In a last step, the one-dimensional arrays with the results of the statistical evaluation were reshaped and re-assembled to its original format, so that each value was assigned to the former cell in the raster, visualising the difference between the arrays in maps for variance and standard deviation per InVEST module per cell in ArcMap.

Results
The six InVEST model modules were applied for Terceira Island. Each module was run with the three input datasets. This process resulted in 18 different ES maps for CLC, COS.A and Sentinel-2-based LULC maps (Figure 3 and Figure 4). Despite the differences amongst the input LULC maps, some models show similar output maps, for example, recreation and pollination. Other model results, for example, carbon storage and nutrient delivery ratio, show differences in the model outputs, indicating differences in the spatial modelled ES supply derived from the differences amongst the input LULC datasets.

Recreation
The InVEST Recreation module quantifies recreation, based on Photo User Days (PUD) uploaded on the online photo-community Flickr* between 2007 and 2017, where 2017 is the last possible year for modelling. The ES "recreation" is supplied predominantly in urban areas, based on the PUD/ha per year, with hotspots in the major cities of Angra do Heroísmo, Praia da Vitória and around Biscoitos (vineyards and coastal site). An exception is located at the touristic spots of "Furnas do Enxofre" (a very popular volcanic cavity). The overview of the three model results showed very similar results ( Figure 3). CLC obtains large areas of no to very weak service supply (0-0.002; 0.003-0.02PUD ha), containing agricultural lands. Continuous urban fabric scored highest for all three input LULC with a total maximum of 22 PUD (0.2 PUD/ha). The city of Angra do Heroísmo, visible as a dark green area in the south of Terceira Island, is one example of this, followed by the urban area around the city of Praia da Vitória in the east. Photo hotspots were found at touristic spots (e.g. Furnas do Enxofre) and natural reserves ("Reserva Natural da Serra de Santa Bárbara e dos Mistérios Negros, Terra Brava e Criação das Lagoas"). Whilst certain areas clearly marked hotspots for recreation and reach between 4 (CLC), 12 (COS.A) and up to 22 PUDs per feature in Sentinel2, averaging the PUDs per hectare results in leaving, for example, moors and heathlands and forested areas highly under-represented (0.004, 0.002 PUD/ha per year). For the CLC input LULC with its large polygons and coarse resolution of the pixels, tourist hotspots appeared large on the map. With smaller feature size, the individual maxima became higher; however, tourist hotspots decreased in size.

Pollination
Pollinator abundance was modelled, based on average pollinator abundance, a dimensionless index considering 0 as no to very low abundance to 1 as maximum abundance. Results showed highest potential for pollinator abundance inland of the Island, with the strongest supply in moors and heathlands and coniferous forests edges (CLC) as seen in Figure 3. Both COS.A and Sentinel2-based LULC maps had a lower average pollinator abundance modelled for the same location around the natural reserves (e.g. "Reserva Florestal Natural Parcial do Biscoito da Ferraria" and "Reserva Florestal Parcial da Serra de S. Barbara e dos Misterios Negros") in the Island's centre, with higher potential. The modelled species abundances in the three model results ranged from an average of 0 in port areas and dump sites to average maxima in grasslands (0.316), moors and heathlands (0.329) and territories mainly occupied by agriculture, with significant areas of natural vegetation (0.21). Forested areas, heath and moorlands scored highest in the CLC dataset, but overall pollinator abundance patterns throughout the three LULC remained similar. Agricultural sites, urban areas and port areas scored lowest with all three LULC datasets, with no to very little pollinator abundance (0 -0.09).

Carbon storage
The carbon storage module calculates a carbon balance for above-ground, below-ground and soil carbon, including carbon stored in dead material. The model computed high capacities for the ES "carbon storage" (Fig. 3) for all three datasets in peat bogs (up to 658 Mg C per ha), followed by forested areas (ranging from 237 -328 Mg C per ha). Urban areas, airport and port areas scored lowest with all three LULC maps. Nonetheless, model results showed diversity amongst the three LULC model outputs. Under CLC, values were slightly higher than in the other two LULC datasets. These findings are in line with Vergílio et al. (2016), who have calculated similar average Carbon storage values for Pico Island, based on a different LULC map. An exception resulted from the Sentinel-2-based dataset, as a result of the limited differentiation of agricultural classes in this LULC datasetdifferentiations of agricultural activities are restricted to "arable land" and "other crops". Due to this difference, the class "Other crops" included vineyards, horticulture and pasture areas. Averaging these three classes lead to an underestimation of vineyards and horticulture areas, but caused an overestimation of pastures (modelled values ranging from 81-180 Mg C ha ), which seemed unrealistically high compared to values of 30 Mg C ha found in literature, for example, by Seó et al. (2017). Potentially, this could be solved by adjusting the input parameters according to the LULC surface, taking into account falsely modelled low values for vineyards and horticultural areas. Despite these differences, all three input LULC maps reflected the low carbon storage capacity of urban areas and present similar spatial distribution. Freshwater-related ES nutrient export, sediment export and water quality (flow retention) showed related results.

Nutrient Delivery Ratio (NDR)
The NDR module spatially depicts the outwash of nitrogen (N) from different LULC types. The ES "nutrient export" shows overall low potential for the Island. As shown in Fig. 4, the -1 -1 model outcomes showed a high degree of dissimilarity between the different input LULC maps. Based on CLC, an annual nutrient export, or N outwash, between 45 and 51 kg N per ha was calculated on intensively cultivated agricultural areas and vineyards on the slopes and hillfoots of the Island, which occurred on 18% of the Island's surface, according to this dataset. N export was modelled with lower loading rates for all other LULC classes. Least annual N export (average of -4.1 -+0.24 kg N y per spatial unit) was modelled for forest and peat bogs in CLC. The NDR, based on COS.A, showed less drastic N outwashes. Results showed strong nutrient outwash on 5% of the Island, where steep slopes coincided with humanly-altered land-use types (-45 and -51 kg N y ), whilst all other areas showed no to very weak nutrient export capacities (-4.1 -+0.24 kg N y ). A strong contrast was presented in the Sentinel-2-based model results. With only two agricultural LULC classes, the model results for these areas were less distinct. The large surface cover of "other crops", a mixture of horticulture, vineyards and pastures, showed annual N export values between -30 and -40 kg N y . The mosaic of different shades of green in the map in Figure 4 reflects the many fragmented habitat types and their N export and retention capacities. InVEST Model results for recreation (Visitation), pollinator abundance (Pollinaiton) and carbon storage (Carbon) for the three input LULC datasets CORINE, COS.A and Sentinel-2.

Sediment Delivery Ratio (SDR)
This module estimates the export of sediment particles (Fig. 4) per ha or watershed per year. It draws upon complex functions of altitude, namely flow accumulation and the Universal Soil Loss Equation (USLE) . The modelled average sediment export on Terceira was relatively low for all three LULC datasets (between 0-15 t ha y ) with few locally specific extremes of up to 55 t ha y in steep agricultural areas towards the centre of the Island. This could be explained by the comparably small scale of agricultural fields, the abundance of support practices, such as stone rows and the high amount of pastures for dairy and meat industry . Agricultural fields showed an estimated average sediment export value ranging from 4.48 t ha y ("Land principally occupied by agriculture" CLC, COS.A) to 7.11 t ha y ("Non-irrigated arable land", CLC; COS.A), findings that are in line with field experiments by Fontes et al. (2004). The model results computed individual maxima per pixel up to 500 t ha y (CLC, COS.A, Sentinel-2) on steep slopes, which is high, considering that losses more than 55 t ha y are assumed to have a high erosion risk in comparable studies. However, calculated per feature in ArcGIS, this averaged to 35-51 t ha y . Soil loss in urban and forested areas was calculated to be higher than 1 t ha y for all three LULC datasets. InVEST Model results for Nutrient export (NDR), Sediment export (SDR) and flow retention (SWY) for the three input LULC datasets CORINE, COS.A and Sentinel-2.

Flow retention (Seasonal Water Yield)
The Seasonal Water Yield module calculates the quick flow water recharge and allows to calculate a flow retention index (1-Qn/P). The ES "flow retention" computed very good retention capacities (0.81 -1) for the majority of the Island, with coastal areas showing slightly lower retention capacities (0.61-0.8) for all three input LULC datasets. The model results (Fig. 4) using CLC showed values ranging between 0.35 to 1.0 [-] for Terceira, with the lowest flow retention values computed for continuous urban fabric (0.35). Highest mean flow retention values were found in natural grasslands (0.95), moors and heathland or urban sports facilities (0.94). Flow retention was lowest using the COS.A dataset with low capacities for streamlines and lower capacities for urban areas: the focus on aquatic ecosystems and riparian zones in the LULC classes might have contributed to this. The highest flow retention was obtained using the Sentinel-2-based LULC dataset: diversity and small size of LULC patterns enhanced retention capacity. Overall, all datasets indicated good flow retention for Terceira. These findings are in accordance with the outputs of the SDR module, showing low average soil exports.

Statistical Results
A statistical analysis of the different datasets reveals trends for each set of model outputs.
The arithmetic mean for all modules and LULC datasets is visualised in boxplots (Fig. 5).
For some modules, such as Visitation, SDR or Flow Retention, the differences between the LULC seemed small. However, the individual histograms revealed that the magnitude of ES supply differed with input LULC data. Regarding pollination, for instance, potential pollinator abundance was modelled to be highest in CLC with maximum values of up to 0.4, with a standard deviation between 0.05 and 0.18. In both COS.A and Sentinel-2 datasets, this deviation of modelled potential abundance is much smaller with a larger number of outliers. The strongest differences were obtained for the Nutrient Export module. The arithmetic means differ largely, with annual N outwash of -7.8 kg ha y using the CLC dataset, -43 kg N ha y for the COS.A map, and -27.9 kg N ha y year when using the Sentinel-2-based LULC dataset.
Nevertheless, Kruskal-Wallis-H proves significant for all six modules and all three model runs. With p-values < 2.2e-16 < 0.05 = α, the null hypothesis can be rejected in all cases and we conclude that there are significant differences amongst the three model outputs and that there is a significant effect of input LULC datasets on model results. This difference is most distinct for the Nutrient Delivery Ratio module, with differences in nitrogen outwash rates of up to 23 kg ha y for the agricultural areas deriving from the high nitrogen loading rates from the CLC dataset. In addition, differences between the LULC maps for the Pollination module were distinct on Terceira Island, as results of the size and type of the different forest patches. For Flow retention, based on the SWY module, the deviation reached up to 0.46, an effect of the aquatic LULC classes present in the COS.A dataset. For the Visitation module, the differences mainly occurred in a few small spots, showing overall the highest consensus of the three input LULC datasets (Fig.  6). In Fig. 6, the distribution of the dissimilarities for the three output maps per InVEST -1 -1 -1 -1 -1 -1 -1 -1 module are shown, highlighting where the differences, based on input LULC, are spatially located on Terceira Island.

Discussion
The applications of the six different InVEST model modules on Terceira Island present an overview of the spatial distribution of ES, based on three different input datasets. The statistical analysis showed that the choice of input LULC data largely affects the model outcomes. The maps of the ES model outputs of CLC, COS.A and Sentinel-2 clearly demonstrated significant differences in terms of the modelled distribution of ES supply. Hence, the decision on the input LULC data largely affected the modelled distribution of ES in a case study region. In the following, the reasons for the differences, as well as some practical guidance on what factors to look for when choosing an appropriate LULC dataset, will be given.
For those InVEST modules that model ES, based on specific indicators, such as recreation (Photo User Days), input LULCs only determined the spatial extent of ES supply as seen in the Visitation module. Differences arise, based on the size of LULC features when visualised as maps. Due to the overall low number of PUDs throughout large parts of the Island, differences between the input LULCs are highly local and small, as shown in Fig. 6.
As the InVEST modules require LULC as the major input for their modelling, modules such as Carbon Storage, Visitation or Pollination, showed little differences between the three output maps -potential reason for this could be the similar type of land use. Whether the classification foresees an area to be covered by "broad leaved forest" (CLC and COS.A) or "Pittosporum" (Sentinel-2) or "transitional woodland-shrub" (CLC), "shrubland, bushland, heathlands" (COS.A) or "Shrub peatland" (Sentinel-2), the corresponding data, for example, Carbon Storage values, remain similar and therewith, differences in ES distribution remain rather small, where ES classes closely resemble each other. For example, the outcomes of the pollinator abundance reflect similar trends throughout all three LULC maps. For other modules drawing upon more complex calculations and combinations of input data, the differences in ES maps are substantial. For example, the results of the NDR model show higher deviations between the individual datasets (Fig. 5). However, the extent of this difference becomes only visible on a map (Fig. 6). With an increasing number of model parameters, the differences between the input LULCs became more distinct -here, the choice of input LULC thus becomes more decisive.
The differences in model results can -as was expected -be explained by the use of different input datasets. Differences in input data influence the model outcomes. The different categorisation of LULC classes and the MMU of each dataset are important for spatial accuracy. The InVEST models assume that the supply of ES changes linearly with the land-use change (Balvanera et al. 2017). With changes in input LULC data, for example, different land uses, the definition of LULC, its spatial resolution, precision and extent, the modelled distribution of ES also changes (e.g. resolution is coarsest in the CLC, followed by COS.A and Sentinel-2-based LULC map). As stated by Gil et al. (2012), an island-scaled LULC cartography with high spatial resolution and with an adequate focus on biodiversity and natural resources protection and management is needed to better address the LULC and derived ES characterisation, evolution and change issues in the Azores Islands. Nevertheless, for ES maps indicating broader trends over larger areas, a coarse resolution may be sufficient. To map a first spatially-explicit distribution of ES, such as recreation (visitation), carbon storage capacity or water quality, the CLC dataset proved sufficient. This is in line with Hamel et al. (2017) who proposed to execute the SDR model on a watershed scale.
Based on our results, the choice of input LULC datasets depends on different factors. The proper selection not only depends on data availability, but also on: 1) the research or policy/decision-making question guiding the modelling study, 2) the ecosystems to be mapped, but also on 3) the spatial resolution of the mapping and 4) data availability at the local level.
Following the MAES guidelines, the purpose of an ES modelling and mapping exercise is often linked to a guiding decision-making or policy question (Maes et al. 2012). Such policy questions often include knowledge requests, policy support questions, questions on resources and responsibilities, application questions and technical and methodological guidance questions (Maes et al. 2018) . These questions often bear implications for the scope, scale and resolution of the modelling exercise.
Another factor of importance for the selection of input LULC data is the ecosystem(s) in the focus of the mapping exercise. Here, the level of detail in the classification scheme adopted for each LULC dataset is important for the model output. For example, the carbon storage maps increase in detail with more specific LULC types. Especially for agricultural areas, the difference becomes visible. Both CLC and COS.A datasets contain >=4 subcategories, whereas the Sentinel-2-based LULC map only uses two out of 13 LULC types for agriculture. This is important, when ES maps are to be used for urban or agricultural areas.
The spatial resolution of the modelling exercise is closely linked to the investigated ecosystem. To model ES that highly depend on the structure and composition of the natural mosaic landscape, with its interactions and processes, this study recommends to select habitat-specific LULC, such as the Sentinel-2-based LULC map. Particularly ES, such as pollination, require spatial information on (small) patch sizes of different ecosystem types. It is recommended to use the highest possible spatial resolution for modelling pollinator abundance, even though the results of all three input LULC datasets show a high degree of similarity. A comparison of these results with the work of Picanço et al. (2017), who mapped spatial distribution of insect pollination in Terceira Island, shows similarities even though both studies apply different methodologies (sampling versus modelling approach).
Depending on the size of the study area, the application of different LULC datasets can be useful. For national ES accounting and large scale ES modelling in the European Union, CLC is recommended, as its large average feature size might be sufficient. For regional and local ES assessments, LULC data, such as COS.A and Sentinel-2-based maps, can be used, especially for smaller areas with detailed landscape structures or feature sizes, including small patches and locally-specific ecosystem types. Examples of this are the EU's Outermost Regions and Overseas Countries and Territories, for which CLC is available, but this is unable to capture the dominant local ecosystem types ranging from tropical rainforest to arctic steppe (Petit and Prudent 2008). In this case, the modeller should strive for maximum class types in, for example, aquatic, forest or urban classes for highest output accuracy (Fig. 7).
Lastly, data availability at local level can impact the choice of LULC. With abundant data on a local level, it is possible to model ES for locally-specific ES types, as the example of Sentinel-2 with its endemic ecosystem types shows. Where such local data are limited, COS.A or CLC can be suitable options to conduct the modelling, as reference data on European level can be used as a proxy. However, this comes at the cost of model robustness and introduces uncertainties to the model outputs. Such effects can be severe for small islands as this affects the modelled distribution of ES and hence, the quality of ES maps.
The approach taken in this paper, running ES models with different available LULCs, attempts to minimise uncertainties, based on LULC input datasets. Each model entails uncertainties. As in all computer-based model approaches, model outputs are strongly dependent on the precision, quality and quantity of model input data. Many of the InVEST models require only few data inputs -for example, the carbon storage or visitation modules -which constitute an advantage of InVEST. This makes it broadly applicable across a variety of social-ecological contexts (Balvanera et al. 2017). Often, simple models bear the risk of being incorrect, but may perform well in matching the general patterns of the data, especially when there is a high variance and uncertainty in the data (Costanza 1989). However, if InVEST models rely on such few data inputs, these data need to be as accurate as possible. Nevertheless, on such a local scale as for small (oceanic) islands (as most EU Outermost Regions and Overseas Countries and Territories), many of the required data are not available. As a result, surrogate or proxy data need to be used, which can have a high uncertainty when InVEST models are applied with coarse and simplified secondary data (Balvanera et al. 2017). In short, the more precise and accurate the input LULC data and the better suited the quality of the data for the individual assessment are, the more precise and reliable the modelled maps of ES distribution will be. As this study shows, running ES models with different available LULC datasets can help to add an additional layer to the ES maps -it can predict where the distribution of ES can be modelled with high certainty and where uncertainty in spatial distribution occurs, based on LULC. This knowledge of the degree of uncertainty in ES maps can enhance the quality of modelled ES maps.

Conclusions
This study shows that the choice of input LULC datasets can have a significant impact on the outcomes of ES maps computed with the InVEST model suite on small islands, such Decision tree for the selection of input LULC datasets for Terceira Island for the InVEST model suite.
as the EU Outermost Region of Terceira on the Azores Archipelago. Comparing three input LULC datasets and six InVEST Modules, significant differences are found between each input LULC. The choice for a particular LULC can either lead to visually enlarging or visually diminishing areas of ES supply on the ES maps. Furthermore, this choice can affect the magnitude of ES supply through the inclusion or omission of certain LULC types. While studies acknowledge sensitivities of the models to input parameters, our work highlights the implications of selecting proper LULC input data -a novel aspect. The use of different input LULC maps in the modelling process can enhance the accuracy of ES maps. Studies and researchers should not only include information on their input parameters, but also on the input LULC dataset with its different classes and feature sizes in order to ensure transparency of the maps for potential users. This is especially relevant if policy and decision-makers or land-use planners are to base their decisions on ES model results.