Ecosystem Services in a Protected Mountain Range of Portugal: Satellite-Based Products for State and Trend Analysis

Mountains are facing strong environmental pressures, which may jeopardize the supply of various ecosystem services. For sustainable land management, ecosystem services and their supporting functions should thus be evaluated and monitored. Satellite products have been receiving growing attention for monitoring ecosystem functioning, mainly due to their increasing temporal and spatial resolutions. Here, we aim to illustrate the high potential of satellite products, combined with ancillary in situ and statistical data, to monitor the current state and trend of ecosystem services in the Peneda-Gerês National Park, a protected mountain range in Portugal located in a transition climatic zone (Atlantic to Mediterranean). We focused on three ecosystem services belonging to three broad categories: provisioning (reared animals), regulating (of water flows), and cultural (conservation of an endemic and iconic species). These services were evaluated using a set of different satellite products, namely grassland cover, soil moisture, and ecosystem functional attributes. In situ and statistical data were also used to compute final indicators of ecosystem services. We found a decline in the provision of reared animals since year 2000, although the area of grasslands had remained stable. The regulation of water flows had been maintained, and a strong relationship with interannual precipitation pattern was noted. In the same period, conservation of the focal iconic species might have been affected by interannual fluctuations of suitable habitat areas, with a possible influence of wildfires and precipitation. We conclude that satellite products can efficiently provide information about the current state and trend in the supply of various categories of ecosystem services, especially when combined with in situ or statistical data in robust modeling frameworks.


Introduction
Mountains are among the most important ecosystems for the supply of ecosystem services (ES), which are valued by mountain inhabitants as well as by people living in lowlands and mountain visitors [1,2].For example, they are known as "water towers" because of their strong enrollment in the hydrological cycle, namely their important role in storage and gradual release of water [3].They are essential refuges for biodiversity and hold a remarkable diversity of landscapes and cultural heritage.However, many mountain areas across Europe are undergoing rapid socio-ecological transitions following the decline of traditional land management systems [4].Extensive land abandonment reflects the decline in importance of provisioning services (food, fodder, wood) from local ecosystems explored through traditional farming due to external inputs of commodities and demographic changes.Generally, this has resulted in increased landscape vulnerability to wildfires, invasive species, and other hazardous events, although positive effects can also be pointed out, such as the decrease in soil erosion and the reduction in agricultural diffuse pollution [5,6].Satellite-based products have been used to provide information about the area of farmland abandonment in Europe, estimated to be up to 7.6 Mha from 2001 to 2012, in part occurring in mountains [7].Future scenarios for mountains of northern Portugal at the transition between Atlantic and Mediterranean climates suggest a trend of declining agricultural activities to be replaced by forestry as the major source of income [8].
Mountains are often classified as protected areas, holding an essential role in the conservation of natural and cultural heritage and in the supply of ecosystem services.Tourism activities, especially those based on agro-tourism, eco-tourism, and wildlife (in particular symbolic species) in protected mountain areas, has been increasing in the last decades [9].Several studies have shown that existing networks of protected areas supply considerable quantities of ecosystem services.Castro et al. [10] reported that two semiarid mountains in Spain contribute to safeguard both biodiversity and ecosystem services, in particular carbon stocks and groundwater recharge.Ecosystem functional types (EFTs), computed from satellite products such as those derived from MODIS (Moderate Resolution Imaging Spectroradiometer), may support the analysis of current protected areas under the new conservation paradigm that considers ES conservation targets to achieve biodiversity conservation, as evaluated by Cabello et al. [11] in Portugal and Spain.Recently, an extensive list of ecosystem functions that can be evaluated by satellite remote sensing products was outlined by Pettorelli et al. [12] with the aim of fostering the efficient monitoring of biodiversity and ecosystem services.
Earth observations, namely satellite products, have been used to support ES analysis for land management and policy options [13,14].Recently, efforts are being developed to select essential ecosystem services variables (EESVs) under the GEO BON (Group on Earth Observations Biodiversity Observation Network) initiative [15].Earth observations (EO) are at the genesis of the concept of essential variables, which arise to characterize changes in the systems, and have focused primarily on climate, followed by the ocean and biodiversity [16], and finally on ecosystem services [15] and Sustainable Development Goals (SDG) [17].Traditionally, satellite-based products provide information about the biophysical components of ecosystem services supply, namely the processes and functions that sustain human well-being, such as productivity, evapotranspiration, and temperature regulation [14,18].The demand side of ES is often very difficult to monitor from space, although some attempts (with limitations) can be made using satellite products to estimate, for instance, population density and human settlements [15].In addition, the mapping and analysis of the full range of ecosystem services is currently beyond the capabilities of remote sensing alone and should thus rely on complementary information from other sources, such as statistics and in situ data [13].
In this context, the objective of this study is to illustrate how satellite products, combined with other sources of information (in situ measurements, statistical data), can be used to monitor the state and trends of several categories of ecosystem services (provisioning, regulating, and cultural).This study is intended as a contribution to the development of EESVs, previously outlined by Cord et al. [15], and a contribution to the characterization of ecosystem services supply in mountain areas using satellite products, either through direct interpretation of the products themselves or as input variables to inform ecological models.To our knowledge, this is the first study assessing one ecosystem service per category using different sources of satellite products in a protected mountain area.The provision of reared animals, the regulation of water flows, and the cultural service of supporting the conservation of a narrow endemic and iconic plant species (Iris boissieri) were evaluated following guidelines suggested by the CICES (Common International Classification of Ecosystem Services).Satellite-based products (grassland cover with Sentinel-2 and Landsat, Sentinel-1 and ESA CCI soil moisture, and Ecosystem Functional Attributes (EFAs) from MODIS) were computed and analyzed together with in situ (presence data on Iris boissieri, field measurements of soil moisture, grassland occurrence data) and statistical sources (precipitation, temperature, number of cattle heads, burnt area) to obtain ES estimates.This study partially builds on satellite products and ecological models previously developed by the authors under the H2020 project Ecopotential (www.ecopotential-project.eu/),extending their application to the assessment of ecosystem services.Finally, the monitoring of ecosystem services using satellite-based products is discussed under the scope of EESV development and sustainable management of mountain areas.

Study Area
The Peneda-Gerês National Park (PGNP) covers around 70,000 hectares and was founded in 1971 (Figure 1).The park includes a group of neighboring mountain ranges, with the highest elevations reaching 1545 m.It is located in a climatic transition zone (the climate is temperate/Atlantic sub-Mediterranean), with 1-2 dry months in summer and mean annual temperature of 13 • C. Annual rainfall is usually above 1500 mm, reaching 3000 mm in summits.As this is the region where it rains the most in Portugal, several artificial reservoirs were built in the main valleys for water storage and hydropower production.Granite is the dominant bedrock type and the predominant soils are leptosols and regosols, usually with low agricultural potential.The main land cover types in highlands are heathland and shrubland, rock outcrops, deciduous oak forests, and meadows.Small-scale agriculture, shrubland, and mixed forests of pine and eucalypt dominate the lowlands.evaluated following guidelines suggested by the CICES (Common International Classification of Ecosystem Services).Satellite-based products (grassland cover with Sentinel-2 and Landsat, Sentinel-1 and ESA CCI soil moisture, and Ecosystem Functional Attributes (EFAs) from MODIS) were computed and analyzed together with in situ (presence data on Iris boissieri, field measurements of soil moisture, grassland occurrence data) and statistical sources (precipitation, temperature, number of cattle heads, burnt area) to obtain ES estimates.This study partially builds on satellite products and ecological models previously developed by the authors under the H2020 project Ecopotential (www.ecopotential-project.eu/),extending their application to the assessment of ecosystem services.Finally, the monitoring of ecosystem services using satellite-based products is discussed under the scope of EESV development and sustainable management of mountain areas.

Study Area
The Peneda-Gerês National Park (PGNP) covers around 70,000 hectares and was founded in 1971 (Figure 1).The park includes a group of neighboring mountain ranges, with the highest elevations reaching 1545 m.It is located in a climatic transition zone (the climate is temperate/Atlantic sub-Mediterranean), with 1-2 dry months in summer and mean annual temperature of 13 °C.Annual rainfall is usually above 1500 mm, reaching 3000 mm in summits.As this is the region where it rains the most in Portugal, several artificial reservoirs were built in the main valleys for water storage and hydropower production.Granite is the dominant bedrock type and the predominant soils are leptosols and regosols, usually with low agricultural potential.The main land cover types in highlands are heathland and shrubland, rock outcrops, deciduous oak forests, and meadows.Small-scale agriculture, shrubland, and mixed forests of pine and eucalypt dominate the lowlands.This site is the only national park in the country, with over 45 years of experience in conservation and management.Since 1997, it is also part of a transnational park and of a UNESCO Biosphere Reserve (Gerês-Xurés).It is also included in the European Natura 2000 network.It holds one of the highest diversity of Annex I habitat types in the country, several Iberian endemic plant and animal species of high conservation value, including Iris boissieri, the Gerês lily, an endemic plant listed in Annex IV of the Habitat Directive and one of the Park's iconic species.Peneda-Gerês also shows the highest number of satellite-based EFTs of several protected areas in the Iberian Peninsula as well as the highest carbon gains [11].
The area of the park holds several protection levels, with some areas of native forest and low shrub at the top of the mountains having higher protection level and less human disturbance.There is high potential for rewilding in the highest elevations, where human activities (mainly pastoralism) were widespread in the past but occur only occasionally nowadays [19].Traditional farming has been decreasing in the last decades, with declining number of farmers and cattle heads from 1999 to 2009 (dates of the most recent farming censuses), although the cattle heads per farm holding has increased, indicating that subsistence agriculture is not so important as it once was.

Selection of Ecosystem Services and Indicators
The physical and biological characteristics of the mountains in PGNP give exceptional conditions for the provision of ecosystem services [20].Three ecosystem services representing the three main sections-provisioning, regulating, and cultural-were chosen because of their socio-economic importance for the park and for the region (Table 1).For a common understanding of the definition of the three ecosystem services, the CICES classification (v5.1) was used [21].
As a provisioning service, reared animals were chosen because of their importance for the local economy.This service refers to livestock raised in housing and/or grazed outdoors, an indicator already suggested by Maes et al. [22].In addition, cattle in this region include traditional breeds, which produce meat and derived products of high interest not only for the local people but also for consumers downstream.Cattle, such as cows, and horses graze almost in the wild in some areas of the park.The direct indicator for the supply of reared animals is the number of cattle heads from statistical sources.Here, satellite remote sensing was used to map grassland areas, an important ecological clause for growing cattle and to assess the potential supply of the service.
As a regulating service, water flow regulation was selected due to the importance that mountains hold for regulation of the water cycle.In addition, precipitation in these mountains can often reach 3000 mm, fostering natural hazards such as floods and landslides.In turn, even in this wet region, summers are dry, and drought or water shortage may occur in lowland farmland and urban areas.Soil moisture is an important indicator of water regulation and used in many environmental studies, from agricultural and hydrological to climate change [23].
As a cultural service, the protection of an endangered and iconic species (Iris boissieri, or Gerês wild lily) was selected because it represents the characteristics of living systems that have an option or bequest value, something in nature that we want future generations to enjoy or use.The recent version 5.1 of CICES includes much more detail in the categories of cultural services, describing in depth more benefits that can be appreciated by humans from ecosystems [21].Here, the model for habitat suitability of this iconic species as an indicator of a cultural service was based on ecosystem functioning attributes derived from MODIS.

Provision of Reared Animals
Livestock is the indicator for reared animals suggested by Maes et al. [22].Here, livestock was evaluated through data from farming census, namely the number of cattle heads for 1999 and 2009 by parish [24], indicating the service provision and trend evolution.The ecological clause used to evaluate reared animals was grassland cover evolution, obtained through the classification of Sentinel-2 and Landsat images for 2016 and 2002, respectively.In this case, satellite grassland cover is an indicator of potential supply, whereas the actual service indicator is derived from agrarian statistical sources (number of cattle heads).In fact, not always satellite productsprovide all the information needed to characterize the supply of ecosystem services [13].It should be noted a mismatch between dates for monitoring grasslands and livestock due to data availability.Considering that the extent of grassland areas does not change substantially in a few years, we assume that this mismatch did not have an important effect in the analysis.Finally, as an indicator of environmental pressure, the grazing livestock index (GLDI) was calculated as the ratio between livestock and grassland area.
Satellite-based grassland cover maps were derived from a supervised classification with the Random Forest classifier forced by the greenness tasseled component (TCT greenness) of Landsat-7 (2002) and Sentinel-2 (2016) satellite scenes.A set of 50 grassland sites common across scenes were used as training areas for classification.Accuracy assessment followed the good practices protocol of Olofsson [25] and used 150 validation sites based on a stratified random sampling with unequal sample that finally included sites with grass (n = 29) and sites with no grass (n = 121) from the grassland map of 2002.For 2016, corrections to the assignment of class (grass/no grass) were done using visual inspection of very high-resolution imagery (VHR) available in Google Earth [26].For each map, an error matrix giving accuracies unbiased stratified area estimates and uncertainty was obtained.
The grazing livestock density index (GLDI) is an indicator from Eurostat (www.ec.europe.ue) that measures the pressure of livestock farming on the environment.In essence, it is the number of livestock units (LSU) per hectare of fodder area, with LSU being a reference unit to aggregate different cattle species and ages.Following the coefficients from Eurostat, 1 LSU is equivalent to one dairy cow, and the remaining coefficients are lower according to the species and age, e.g., horses are equivalent to 0.8.In our study, LSU was calculated for years 1999 and 2009 per civil parish based on the number of cattle heads.Then, GLDI was calculated using LSU and grassland cover for 2002 and 2016.

Regulation of Water Flows
Soil moisture (SM) is an important variable in several hydrological processes, namely in runoff generation, drought development, and evapotranspiration [23].It was considered an essential climate variable (ECV) in 2010 by the Global Climate Observing System (GCOS), which considers relevant, feasible, and cost-effective variables obtained from satellite products for systematic observations of the climate system [27].The regulation of the flows of water in the environment strongly depends on soil moisture, which may directly influence natural hazards, such as floods and droughts, given a surplus or lack of soil moisture [28].
For the regulation of water flows in PGNP, two satellite soil moisture products at different temporal and spatial resolutions were analyzed: Sentinel-1 soil moisture, with a high spatial resolution (10 m), for the period 2015-2017 and ESA CCI soil moisture, a long-term series for evaluating trends since year 2000, though with coarser spatial resolution (25 km) [28].
Soil moisture is highly variable in space and time, therefore Sentinel-1 retrieval is an important tool to monitor hydrological processes and ecosystems services related to water flow, in particular when in situ observations are rare and in specific locations not representative of local conditions.Nevertheless, soil moisture based on Sentinel-1 retrievals can accurately represent the degree of wetness only for the first ~5 cm of soil layer, therefore only cells with agriculture and low shrubland as predominant land cover are eligible.This information at high scale resolution (10 m) may be very useful for farmers to know the degree of moisture in their fields to plan irrigation or for authorities in the calculation of fire risk based on 6 to 12 days retrievals.
The Copernicus Sentinel-1 mission is composed of two satellites-Sentinel-1A and Sentinel-1B-that carry identical synthetic aperture radar (SAR) systems, operating at a frequency of 5.4 GHz (5.5 cm), with a revisiting time of six days.In this study, a retrieval model based on a machine learning approach developed under the project Ecopotential [29] was used for PGNP.It relies on support vector regression (SVR) to construct an empirical model considering σ0 (Sentinel-1) combined with auxiliary ground information (land cover and topography), as well as coarse-resolution-modeled surface soil moisture from the Global Land Data Assimilation System [30], which acts as a first guess and constrains the model estimations.The model was optimized for global application using data from the International Soil Moisture Network [31].The average global estimation accuracy was shown to be at 0.08 m 3 m −3 , but the comparison with time series of in situ monitoring stations showed that it is around 0.05 m 3 m −3 for many stations.Soil moisture estimations with the abovementioned approach can be performed with the Python package PYSMM [32], which is freely available and is open source.
Sentinel-1 soil moisture retrievals were averaged for PGNP and converted to monthly values to smooth the line and to be compared with precipitation and temperature local data collected at SNIRH (National Information System for Water Resources).Analysis were carried out in R statistical software [33].
For validating Sentinel-1 soil moisture retrievals in PGNP, a field campaign was carried out during the summer of 2016 in three dates synchronized with the Sentinel-1 orbit.Twenty sites were established and measured repeatedly in each date, but only eight sites were used for final comparison because the others fell in areas where the retrieval was not calculated (forest, tall shrub).Each sample was collected in the top 20 cm with a soil sampler, stored in an aluminum soil sample box (120 cc), weighted with a digital scale, and stored in a cooler.In the lab, samples were dried during 48 h at 105 • C, weighted, and moisture content was estimated in percentage using the following equation: where WW is wet weight and DW is dry weight.The values were compared through the coefficient of determination (R 2 ) and the percent bias (PBIAS) for each pixel during the surveyed periods.As a result, average R 2 for the eight final pixels was 0.88, and average PBIAS was −13%, indicating that satellite Sentinel-1 was overestimating the in-field observed values for soil moisture.The pixel that performed better through time was number 16 (Figure 2).Soils are very thin in this region and local conditions, such as light rain in the morning and clear sky in the afternoon, may strongly influence the values of moisture at the top soil level.Even though the model was tuned with globally distributed in situ soil moisture from the International Soil Moisture Network (ISMN), model accuracy for local conditions of PGNP can be considered satisfactory.
As part of the Climate Change Initiative (CCI) program, the European Space Agency (ESA) released the first multidecadal, global satellite-observed soil moisture dataset, which combines single sensor active and passive microwave into three products: merged ACTIVE, merged PASSIVE, and COMBINED [28].Here, the COMBINED daily product version 03.3-the weighted average of microwave data for both ACTIVE and PASSIVE-was used, taking into account the error characteristics of the individual data sets [34].
As ancillary data to characterize environmental pressure, monthly precipitation was used from the E-OBS gridded dataset for Europe and extracted for the PGNP, with a spatial resolution of 25 km [35].In order to analyze trends in precipitation and soil moisture, the time series were first tested for breakpoints in the trend component using BFAST (Breaks For Additive Season and Trend) [36].After the "remainder" component was removed to reduce the influence of anomalies, the time series were tested for monotonic linear trends using seasonal Mann-Kendall tests [37].Both were performed in R statistical software [33].
Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 18 were tested for monotonic linear trends using seasonal Mann-Kendall tests [37].Both were performed in R statistical software [33].

Habitat Suitability for the Iconic Species Iris boissieri
Cultural services are primarily defined as the nonmaterial benefits people obtain from ecosystems [38].Recently, CICES described these services as the "outputs of ecosystems (biotic and abiotic) that affect physical and mental states of people" [21].These intangible benefits can be evaluated by the relative contribution of ecological structures and functions specified in the biophysical domain and by applicable social evaluation approaches [39].Focusing on the functional approach and given the importance of biodiversity conservation in PGNP-in particular the preservation of narrow endemic and iconic (flag) plant species, the Gerês wild lily (Iris boissieri)-the availability and trends of suitable habitat areas can be a proxy indicator for this cultural service considering that the presence of this species and its habitat indirectly contributes to moral well-being [21].We selected the Gerês lily not only for its current "critically endangered" conservation status (http://www.iucnredlist.org/details/162312/0)but also because the only information available so far about the ecology of this rare and narrowly distributed species is based on field surveys covering essentially the core of its distribution area.A comprehensive field survey carried out in 2007 across the National Park returned indications of a reduced and fragmented population size (area of occupancy well below 10 km²) and in continuous regression, possibly due to various threats, such as changes in land use (e.g., abandonment of pastoral systems), recreational activities, and changes in fire regime.Moreover, the species is listed in Annex IV of the EU Habitat Directive, and therefore under mandatory, regular monitoring of its conservation status by EU member states (according to Article 17 of the HD).
In order to encourage the preservation of ecological conditions necessary for sustaining populations of Iris boissieri, we developed a model to calculate Predicted Suitable Areas (PSA) at 1 km 2 cell size.The methodology was based on Arenas-Castro et al. [40], which showed that EFAs can perform as good or better than traditional habitat predictors (land cover and climate) for Iris boissieri.Model accuracy was measured as the area under the curve (AUC) of receiver operator characteristic (ROC) curves.AUC is a robust threshold-independent measure of a model's ability to discriminate presence from absence [42], ranging between 0 and 1 (models with values below 0.7 were considered poor, 0.7-0.9 were moderate, and >0.9 were good).The resulting models with the highest AUC median (among those satisfying the condition ≥0.7) were selected.We used AUC

Habitat Suitability for the Iconic Species Iris boissieri
Cultural services are primarily defined as the nonmaterial benefits people obtain from ecosystems [38].Recently, CICES described these services as the "outputs of ecosystems (biotic and abiotic) that affect physical and mental states of people" [21].These intangible benefits can be evaluated by the relative contribution of ecological structures and functions specified in the biophysical domain and by applicable social evaluation approaches [39].Focusing on the functional approach and given the importance of biodiversity conservation in PGNP-in particular the preservation of narrow endemic and iconic (flag) plant species, the Gerês wild lily (Iris boissieri)-the availability and trends of suitable habitat areas can be a proxy indicator for this cultural service considering that the presence of this species and its habitat indirectly contributes to moral well-being [21].We selected the Gerês lily not only for its current "critically endangered" conservation status (http://www.iucnredlist.org/details/162312/0)but also because the only information available so far about the ecology of this rare and narrowly distributed species is based on field surveys covering essentially the core of its distribution area.A comprehensive field survey carried out in 2007 across the National Park returned indications of a reduced and fragmented population size (area of occupancy well below 10 km 2 ) and in continuous regression, possibly due to various threats, such as changes in land use (e.g., abandonment of pastoral systems), recreational activities, and changes in fire regime.Moreover, the species is listed in Annex IV of the EU Habitat Directive, and therefore under mandatory, regular monitoring of its conservation status by EU member states (according to Article 17 of the HD).
In order to encourage the preservation of ecological conditions necessary for sustaining populations of Iris boissieri, we developed a model to calculate Predicted Suitable Areas (PSA) at 1 km 2 cell size.The methodology was based on Arenas-Castro et al. [40], which showed that EFAs can perform as good or better than traditional habitat predictors (land cover and climate) for Iris boissieri.In short, a correlative model based on presence-only (and pseudo-absences) for the test species and on interannual average values of EFAs was calibrated for the period 2001-2007, and annual projections from 2006 to 2001 (backcasting) and from 2008 to 2016 (forecasting) were carried out using biomod2 R package [41].Data on Iris boissieri (72 records) resulted from a local survey carried out in 2007.Six EFAs derived from MODIS (Terra MODIS products, collection 6; 2001-2016 time series) were used to predict PSA for each year, calculated from Enhanced Vegetation Index (EVI) (MOD13Q1), Land Surface Temperature (LST) (MOD11A2), and Albedo (MCD43A3).The two best scored EFAs were EVImn = EVI annual minimum and LSTdmxc = Cosine of the date of maximum LST.Model accuracy was measured as the area under the curve (AUC) of receiver operator characteristic (ROC) curves.AUC is a robust threshold-independent measure of a model's ability to discriminate presence from absence [42], ranging between 0 and 1 (models with values below 0.7 were considered poor, 0.7-0.9 were moderate, and >0.9 were good).The resulting models with the highest AUC median (among those satisfying the condition ≥0.7) were selected.We used AUC median values for excluding models with lower scores, for the binary transformations of model probability predictions, and for building the ensemble predictions.Differences between the mean predictive ability of the best models in terms of AUC were also compared using Tukey's HSD (Honest Significant Difference).Model evaluation was based on cross-validation, with the species datasets divided into 80% of the records for model calibration and 20% for model evaluation.To transform predicted probabilities into predicted suitable/unsuitable areas, we used a threshold minimizing the straight-line distance between the ROC plot and the upper-left corner of the unit square [43].We finally reclassified the binary PSA pixels considering their stability as suitable areas (i.e., number of years predicted as suitable between 2001 and 2007) to build a map with a range of suitable areas.PSA response curves were plotted for the EFAs with the highest contribution for the model.Finally, we explored the possible effect of fire occurrence (socio-environmental pressure) on annual PSA projections based on official fire data obtained from the Institute for Nature Conservation and Forests of Portugal.

Reared Animals
Livestock units and area of grasslands were the indicators for actual and potential ES supply, respectively.We found a decline in the number of livestock units from 1999 to 2009 (−18,000 units, −30%), in particular for goats, sheep, and cows (Figure 3A).In contrast, the number of horses increased in this period, most probably due to financial incentives [24].The area of grasslands has basically stabilized in the parishes of PGNP, with only a small decline (−34 ha) verified between 2002 and 2016.Therefore, it is likely that the actual supply of the reared animal service is decreasing although the potential supply side has stabilized.
This finding is consistent with the ongoing mountain agricultural abandonment process, widely documented across Europe [4].Migration and aging of populations has been leading to the collapse of traditional farming systems in PGNP [44].In fact, local communities in PGNP recognize a decline in food production mostly due to land abandonment, with negative effects on biodiversity, cultural heritage, and fire risk, as shown by local participatory studies [45,46].Accordingly, the environmental pressure expressed by the GLDI has decreased from 3.1 to 2.6 livestock units by hectare on average for all parishes in PGNP.At the local level, there is a general trend for decreasing GLDI; however, one parish (Gavieira) increased GLDI for more than seven livestock units per hectare in 2016, indicating a trend into localized and specialized livestock farming systems in the region, as also reported by Poças et al. [47].Hence, the agrarian census reports higher number of cattle heads per farm holding in 2009 compared to 1999.Some parishes with low GLDI may supply grass (or grazing land) to farm holdings in other parishes with more cattle.Although cattle are decreasing in the PGNP, some species such as horses and cows graze in a free-ranging regime throughout the year, with possible conflicts (competition for resources) with wildlife, namely the wild goat, as reported by Moço et al. [48].The major implication of the decrease in reared animals is the likely decrease in food supply in a region that holds potential (grasslands, local breeds) for growing cattle in a traditional and sustainable way.
Another important issue arising from agricultural land abandonment is the increase in fire risk due to the decline in pastoral activities and the resulting fuel accumulation across the landscape [49].Grazing systems were acknowledged by different stakeholders across the Mediterranean to mitigate fire [5].The increase in vegetation biomass due to shrub and forest encroachment in abandoned land, often by fast-growing species such as eucalypts and invasive acacia trees, is creating optimal conditions for fire occurrence and an increase of burnt area, as can be seen since 2003 in PGNP (Figure 6C).From the point of view of biodiversity, land abandonment may negatively influence plant diversity in most farmlands of PGNP, as reported by Lomba et al. [50], but little impact was noticed on bird species diversity, suggesting adaptation of these communities to farmland abandonment [51].

Regulation of Water Flows
Soil moisture is strongly dependent on precipitation and temperature.In PGNP, summers are relatively dry, but high temperatures determine low soil moisture values, as shown during the period July-September of 2016 (Figure 4).In turn, higher precipitation events influenced the peak in soil moisture between December 2015 and April 2016.Although the average value of February 2016 was slightly higher than the average of February 2015 (Figure 4B), the maps show stability in the darker tone (higher soil moisture) in February 2015 than in February 2016.However, in February 2016, there are darker tone spots, corresponding to very high soil moisture pixels that may influence the higher average value.
Continuous time series allow trends or sharp changes in ecosystem service potential supply to be analyzed.Precipitation showed a generally flat linear trend without significant breakpoints since 2000 (Figure 5).On the other hand, a major breakpoint was detected for soil moisture around January 2004 (with 95% confidence interval between October 2003 and May 2004), setting a shift in the linear trend from significantly decreasing to significantly increasing (Figure 5).The water regulation service (represented by soil moisture) has been increasing in PGNP since 2004, which is important for drought mitigation but also for the buffering of high precipitation periods, when compared to the winter of 2000/2001 in which major precipitation events had led to high soil moisture content and consequently led to natural hazards events in the area, such as major landslides.

Regulation of Water Flows
Soil moisture is strongly dependent on precipitation and temperature.In PGNP, summers are relatively dry, but high temperatures determine low soil moisture values, as shown during the period July-September of 2016 (Figure 4).In turn, higher precipitation events influenced the peak in soil moisture between December 2015 and April 2016.Although the average value of February 2016 was slightly higher than the average of February 2015 (Figure 4B), the maps show stability in the darker tone (higher soil moisture) in February 2015 than in February 2016.However, in February 2016, there are darker tone spots, corresponding to very high soil moisture pixels that may influence the higher average value.
Continuous time series allow trends or sharp changes in ecosystem service potential supply to be analyzed.Precipitation showed a generally flat linear trend without significant breakpoints since 2000 (Figure 5).On the other hand, a major breakpoint was detected for soil moisture around January 2004 (with 95% confidence interval between October 2003 and May 2004), setting a shift in the linear trend from significantly decreasing to significantly increasing (Figure 5).The water regulation service (represented by soil moisture) has been increasing in PGNP since 2004, which is important for drought mitigation but also for the buffering of high precipitation periods, when compared to the winter of 2000/2001 in which major precipitation events had led to high soil moisture content and consequently led to natural hazards events in the area, such as major landslides.The ESA CCI soil moisture product can be criticized as its spatial resolution might be too coarse for local evaluations.However, it was the only product available for analyzing a fairly extensive period of time.In fact, long-term ESA CCI SM has shown good performance for characterizing soil moisture, revealing the effectiveness of merging active and passive soil moisture products, as shown by Cui et al. [52].With future climate projections indicating an intensification of dry periods in the region with consequences for water regulation services [53], Sentinel-1 soil moisture with high temporal and very high spatial resolution may be a robust instrument to support operational planning for natural hazards such as drought and fire.The ESA CCI soil moisture product can be criticized as its spatial resolution might be too coarse for local evaluations.However, it was the only product available for analyzing a fairly extensive period of time.In fact, long-term ESA CCI SM has shown good performance for characterizing soil moisture, revealing the effectiveness of merging active and passive soil moisture products, as shown by Cui et al. [52].With future climate projections indicating an intensification of dry periods in the region with consequences for water regulation services [53], Sentinel-1 soil moisture with high temporal and very high spatial resolution may be a robust instrument to support operational planning for natural hazards such as drought and fire.[28] and precipitation from E-OBS [35] for the period 2000-2016 in PGNP.

Habitat Suitability for the Iconic Species Iris boissieri
Suitable habitat conditions for Iris boissieri are mostly located in the core areas of the National Park in the Peneda (north) and Gerês (south) mountain ranges (Figure 6A).The spatial projections of the model based on six satellite EFAs also allows the pinpointing of the populations that might be at risk due to recent shifts or fluctuations in habitat suitability.The fact that most recorded Iris presences coincide with areas predicted as higher stability classes (pixels that ranged from 4-7 of PAS reached 93.4% of positive validation with 2007 presence records) confirms the high predictive performance of the model [40].
Nevertheless, PSA projections show high inter annual variability, explained by variations of the two most important predictors-EVImn and LSTdmxc-but also linked to fire occurrence (Figure 6B,C).Specifically, higher EVImn values are usually related to lower PSA values (Figure 6B), indicating that the species prefers habitats with lower minimum productivity.In turn, higher negative values of the cosine of the date of maximum LST are linked to higher PSA values, indicating that the species prefers habitats reaching their maximum LST in summer (Figure 6B).
There were two periods with several consecutive "bad" years for PSA (2003-2005 and 2009-2013; Figure 6C).These low PSA periods seem to coincide with years recording large burnt areas, indicating that fire may negatively influence the suitable habitat conditions for the species,  [28] and precipitation from E-OBS [35] for the period 2000-2016 in PGNP.

Habitat Suitability for the Iconic Species Iris boissieri
Suitable habitat conditions for Iris boissieri are mostly located in the core areas of the National Park in the Peneda (north) and Gerês (south) mountain ranges (Figure 6A).The spatial projections of the model based on six satellite EFAs also allows the pinpointing of the populations that might be at risk due to recent shifts or fluctuations in habitat suitability.The fact that most recorded Iris presences coincide with areas predicted as higher stability classes (pixels that ranged from 4-7 of PAS reached 93.4% of positive validation with 2007 presence records) confirms the high predictive performance of the model [40].
Nevertheless, PSA projections show high inter annual variability, explained by variations of the two most important predictors-EVImn and LSTdmxc-but also linked to fire occurrence (Figure 6B,C).Specifically, higher EVImn values are usually related to lower PSA values (Figure 6B), indicating that the species prefers habitats with lower minimum productivity.In turn, higher negative values of the cosine of the date of maximum LST are linked to higher PSA values, indicating that the species prefers habitats reaching their maximum LST in summer (Figure 6B).There were two periods with several consecutive "bad" years for PSA (2003-2005 and 2009-2013; Figure 6C).These low PSA periods seem to coincide with years recording large burnt areas, indicating that fire may negatively influence the suitable habitat conditions for the species, probably associated to unfavorable years of EVImn and LSTdmxc, which are the best scored model predictors.Since fire was not directly included as a model predictor, these relationships may be indirect.In fact, fire may directly cause plant mortality and devastate suitable areas for the species to sprout; however, the co-effects of fire-related environmental variables (such as maximum temperatures, low humidity, or longer dry periods) may also negatively influence PSA.In turn, the highest PSA values, predicted for years 2008 and 2014 (Figure 6C), were associated with lower values of burnt areas.It should be noted, however, that model predictions were based on the presences of Iris in year 2007 only and thus the related uncertainty must be taken into consideration.
Overall, the species distribution model calculated with satellite-based predictors allowed the spatial and temporal monitoring of PSA as an adequate indicator of the stability of suitable habitats for Iris boissieri and consequently for the provision of a cultural service related to the protection of the habitat of an iconic and endangered species.For the time series under study (2001)(2002)(2003)(2004)(2005)(2006)(2007), the more stable PSA pixels (i.e., those that ranged between 4 and 7; see Figure 6A) included 93.4% of the presences recorded during the 2007 field campaign.Future campaigns based on model projections of habitat suitability will allow further validation of the model and of PSA as a suitable ecosystem service indicator.
Other approaches have used expert group evaluation to assess symbolic species for cultural services provision, for instance, in European Alps [9], with relevance because the human perceptions were considered.In the future, the combination of the two approaches might be a reasonable option for modeling cultural services through the habitat of iconic species.Here, satellite-based products to map the potential distribution, as well as the major socio-environmental pressures, are essential to target conservation measures in protected mountain areas.This allows the traditional way of using satellite products to evaluate cultural ecosystem services to be complemented, usually based on feeding a Geographic Information System (GIS), with conservation, tourism, and heritage features [13].

Satellite-Based Indicators for Monitoring Ecosystem Services in Mountains
For the three selected ecosystem services chosen as examples of the three broad categories of services (provisioning, regulating, and cultural), satellite products have shown high potential, together with in situ and statistical data, to support the assessment of ecosystem services in the current status and trend analysis in a mountain area, including the ecological clause (potential supply), service supply, and pressure indicators.Demand was not considered in this study as it is very difficult to monitor demand in all categories of ES, especially for regulating services.Even though examples of how satellite-based products can be used for demand monitoring were identified by Cord et al. [15], ES demand remains a challenge for assessment and monitoring using satellite products.In addition, satellite products cannot always inform us directly about the supply of a service, as it was shown with soil moisture for the regulation of water flows.Satellite products can complement statistical data (reared animals) or serve as input for models (habitat suitability of Iris boissieri, for cultural services).
Another important point is the accuracy of satellite products to evaluate ecosystem services at local scales, in particular for mountains, where cloud cover may influence the quality and the number of available scenes for optical products or topography may influence the radar signal [14].New products and more efficient models retrieving signals may bring along more accurate data, which should be calibrated and/or validated with in-field measurements, especially for local case studies [12].This is the case of Sentinel-1 soil moisture retrieval used in this study for water flow regulation.The product was conceived mainly for global assessments, as the majority are, but local validation still showed sufficient accuracy.A better retrieval could therefore be performed due to the input of more local environmental variables.Local-scale studies, with the use of satellite products to assess ecosystem services provision, are important to enlarge the capacity for local management and decision making.
In recent years, the limitations of the use of satellite products have been dissipating with growing spatial and temporal resolutions of the products, politics of open access data, user-friendly platforms to collect information often using cloud processing, and creating opportunities for democratic use of satellite products, among others advantages.However, significant challenges remain to monitor ecosystem services from space for which five priorities were outlined by Cord et al. [15] and supported by our research: (i) defining standardized and monitorable EESVs; (ii) integrating satellite and socio-economic data; (iii) ensuring open-access, maintenance and operability of satellite products; (iv) assessing spatial mismatches between service supply and demand, trade-offs across regions and global teleconnections; and (v) designing and maintaining opportunities for long-term collaboration across disciplines [15].

Conclusions
Three ecosystem services were evaluated in Peneda-Gerês National Park, a protected mountain range in northern Portugal.Provision of reared animals, regulation of water flows, and cultural services through habitat suitability for an endangered and iconic plant species were analyzed using satellite-based products together with in situ and statistical data.
Grassland cover based on Sentinel-2 (2016) and Landsat-7 (2002) images, together with statistical data on livestock, revealed that the provision of reared animals is decreasing in PGNP, although the grassland area has stabilized since 2002 after a decline since the 1980s.Accordingly, the grazing livestock density index has decreased since 2002 in PGNP; however, in one parish, it has increased, showing a tendency toward specialized livestock farming systems together with generalized land abandonment in the park.
Sentinel-1 soil moisture retrievals were used for mapping the current state of water regulation in PGNP at very high spatial resolution, and ESA CCI soil moisture with lower spatial resolution were used to evaluate the trend since 2000.The provision of the water regulation service remained stable, and there was a strong relationship with interannual precipitation pattern.The very high spatial resolution of Sentinel-1 soil moisture allows local targeting actions to cope with natural hazards, such as drought and fire, particularly under future consequences of climate change.
The cultural service of preserving the habitat of Iris boissieri, an iconic species of this National Park, was computed based on ecosystem functional attributes derived from MODIS.Predicted suitable areas for Iris boissieri were modeled and yearly projected since the year 2001, which showed interannual fluctuations of higher and lower PSA values with possible influence of burnt areas and precipitation.The model based on satellite predictors allowed the spatial and temporal monitoring of PSA for this iconic species as an adequate indicator of the stability of habitats for future conservation and consequently for the provision of a cultural service.
We believe that this study can be replicated in other mountain areas to assess the supply of these three ecosystem services considering that satellite products used here are global and freely available, and the codes to run the models are open.Notwithstanding, local data is important, especially to validate satellite products and to calibrate/validate the models.Overall, satellite-based products were shown to be efficient indicators to support the current state and trend analysis of ecosystem services in a protected mountain range, in particular for increasing the spatial and temporal accuracy of assessment and monitoring.This study contributes to innovative perspectives on the use of satellite-based products on ecosystem services assessments, with relevance for monitoring efforts and in particular for the ongoing definition of essential ecosystem service variables.

Figure 2 .
Figure 2. Comparison of volumetric soil moisture (SM) between Sentinel-1 retrievals and in-field measurements in PGNP for three dates where Sentinel-1 passed the region during the summer of 2016 (7th of June, 12th of July, and 21st of September).
In short, a correlative model based on presence-only (and pseudo-absences) for the test species and on interannual average values of EFAs was calibrated for the period 2001-2007, and annual projections from 2006 to 2001 (backcasting) and from 2008 to 2016 (forecasting) were carried out using biomod2 R package [41].Data on Iris boissieri (72 records) resulted from a local survey carried out in 2007.Six EFAs derived from MODIS (Terra MODIS products, collection 6; 2001-2016 time series) were used to predict PSA for each year, calculated from Enhanced Vegetation Index (EVI) (MOD13Q1), Land Surface Temperature (LST) (MOD11A2), and Albedo (MCD43A3).The two best scored EFAs were EVImn = EVI annual minimum and LSTdmxc = Cosine of the date of maximum LST.

Figure 2 .
Figure 2. Comparison of volumetric soil moisture (SM) between Sentinel-1 retrievals and in-field measurements in PGNP for three dates where Sentinel-1 passed the region during the summer of 2016 (7th of June, 12th of July, and 21st of September).

Figure 3 .
Figure 3. (A) Reared animals service provision demonstrated by the evolution in the number of cattle heads from 1999 to 2009 and the grassland area cover (ha) in 2002 and 2016 in the parishes of PGNP.(B) Grazing livestock density index (GLDI) in the parishes of PGNP, illustrating the level of cattle pressure in the grassland cover area in 2002 and 2016.

Figure 3 .
Figure 3. (A) Reared animals service provision demonstrated by the evolution in the number of cattle heads from 1999 to 2009 and the grassland area cover (ha) in 2002 and 2016 in the parishes of PGNP.(B) Grazing livestock density index (GLDI) in the parishes of PGNP, illustrating the level of cattle pressure in the grassland cover area in 2002 and 2016.

Figure 4 .
Figure 4. Sentinel-1 soil moisture (%) in PGNP.(A) Spatial distribution of soil moisture in two dates of the month of February 2015 and 2016, respectively.Cells with no data (white) refer to cells that the algorithm does not retrieve, such as forest and bare rock land cover classes.(B) Average monthly Sentinel-1 soil moisture for the years 2015, 2016, and 2017 with average monthly local precipitation (mm) and temperature (°C) from SNIRH (snirh.apambiente.pt).

Figure 4 .
Figure 4. Sentinel-1 soil moisture (%) in PGNP.(A) Spatial distribution of soil moisture in two dates of the month of February 2015 and 2016, respectively.Cells with no data (white) refer to cells that the algorithm does not retrieve, such as forest and bare rock land cover classes.(B) Average monthly Sentinel-1 soil moisture for the years 2015, 2016, and 2017 with average monthly local precipitation (mm) and temperature ( • C) from SNIRH (snirh.apambiente.pt).

Figure 6 .
Figure 6.(A) Distribution of predicted suitable areas (PSA) between 2001 and 2007 by classes for each grid unit (1 km 2 ) for Iris boissieri in PGNP.Pictures from Iris boissieri and suitable habitat gently provided by ADERE association (http://www.adere-pg.pt/?lg=2).(B) Response curves of the two main contributing predictors: EVImn (enhanced vegetation index annual minimum) and LSTdmxc (cosine of the date of maximum land surface temperature) with PSA.(C) Interannual variation of PSA (black line) and PSA predictions for 2007, year with observed presences (dashed red line), along with auxiliary data related to the total amount of burned area in the study site each year (grey bars).

Figure 6 .
Figure 6.(A) Distribution of predicted suitable areas (PSA) between 2001 and 2007 by classes for each grid unit (1 km 2 ) for Iris boissieri in PGNP.Pictures from Iris boissieri and suitable habitat gently provided by ADERE association (http://www.adere-pg.pt/?lg=2).(B) Response curves of the two main contributing predictors: EVImn (enhanced vegetation index annual minimum) and LSTdmxc (cosine of the date of maximum land surface temperature) with PSA.(C) Interannual variation of PSA (black line) and PSA predictions for 2007, year with observed presences (dashed red line), along with auxiliary data related to the total amount of burned area in the study site each year (grey bars).

Table 1 .
Ecosystem services evaluated in PGNP according to the CICES (Common International Classification of Ecosystem Services) classification (v5.1), and their respective indicators used in this study (bold), including satellite-based products.The Pressures column is not part of CICES but was included here to list the major socio-environmental stressors for each service.