Analysing carbon sequestration and storage dynamics in a changing mountain landscape in Portugal: insights for management and planning

ABSTRACT We assessed the effects of landscape change on the climate regulation ecosystem service in a mountain river basin of Portugal, through the quantification, valuation and mapping of carbon sequestration and storage. The analyses were based on land use and land cover (LULC) changes that took place between 1990 and 2006 and on expected changes defined by three LULC change scenarios for 2020. We used the Integrated Valuation of Ecosystem Services and Tradeoffs model for scenario building and carbon assessment and valuation, and several modelling tools to assess past, current and future carbon in four different pools. Soil organic carbon data was obtained through an extensive sampling scheme across the entire study area. Recent (1990–2006) and expected landscape changes (2006–2020) affected considerably carbon sequestration and storage. Observed landscape changes generally promoted carbon sequestration and storage, and had a positive effect on the climate regulation ecosystem service, both biophysically and economically. Expected LULC changes further extend the capability of the landscape to increase carbon sequestration and storage in the near future. The carbon sequestered and stored in vegetation and soil contributes to avoid socio-economic damages from climate change, while increasing the economic value of particular LULC classes and the whole landscape. These results are essential to inform land planning, especially on how, where and when changes in landscapes may affect the provision of the climate regulation ecosystem service. EDITED BY Sandra Luque


Introduction
Carbon flows naturally in the Earth System through the atmosphere, biosphere and lithosphere in an ensemble of processes known as the carbon cycle (Ciais et al. 2013). However, emissions of carbon in the form of carbon dioxide, one of the major greenhouse gases, have increased over time both due to the use of fossil fuels for energy and due to historic anthropogenic land use and land cover (LULC) changes. These processes have largely increased its atmospheric concentration, contributing to climate change and increasing the likelihood of environmental and economic losses in the future (Noble et al. 2005;EEA 2012;Ciais et al. 2013). On the other hand, increasing carbon concentration in the atmosphere and global warming may act as drivers of LULC change at longer timescales due to vegetation shifts poleward and changes in the distribution of patterns of agricultural areas (Davies-Barnard et al. 2015), as well as modifiers of the carbon dynamics in the Earth system by CO 2 -fertilization, which may contribute to higher carbon uptake in terrestrial ecosystems (Devaraju et al. 2016). However, uncertainties still persist regarding the effect of CO 2 -fertilization and the increase of temperatures on the balance between ecosystem productivity and respiration losses (Chen et al. 2014), as well as other possible effects of climate change on ecosystems (e.g. the increase of pests and diseases) (Devaraju et al. 2016), thus not ensuring long-term compensation of recent carbon emissions. Hence, the reduction of anthropogenic emissions of greenhouse gases is of utmost importance to balance the composition of the atmosphere and mitigate future damages, as underlined recently in the Paris climate conference (COP21) agreement (Paris Agreement).
The capacity of ecosystems to influence the regulation of the concentration of greenhouse gases in the atmosphere, and therefore climate, is an essential ecosystem service with many benefits for human societies (de Groot et al. 2002(de Groot et al. , 2010MEA 2005), from the mitigation of socioeconomic damages associated with climate change (Conte et al. 2011) to the maintenance or improvement of local economies (De Koning et al. 2011;IIED 2012;Bottazzi et al. 2013;Clark et al. 2014). Terrestrial ecosystems, including forest, seminatural and agricultural ecosystems, play an important role in carbon cycling. They act both as carbon sinks, through carbon sequestration and storage in plant biomass, litter and soil organic matter, and as carbon sources, through emission of carbon from biological processes such as respiration (Chen et al. 2014) and from anthropogenic activities such as land use and its dynamics (Ciais et al. 2013) or the occurrence of disturbances (e.g. wildfires).
The supply of the climate regulation service is directly affected by LULC change (Foley et al. 2005;Nelson et al. 2005) since it acts directly on the exchanges of greenhouse gases between terrestrial ecosystems and the atmosphere, affecting carbon stocks positively or negatively, according to the type of LULC conversion and on the management measures carried out (Zhu et al. 2010;Houghton & Goodale 2013). LULC change is a major driver of carbon emissions contributing to nearly 12.5% of the total global anthropogenic carbon emitted (Houghton et al. 2012). LULC change can affect the economic value of the climate regulation ecosystem service supplied by terrestrial ecosystems and biomes. Recent estimates indicated an approximate value of 2250 int$ ha −1 yr −1 (2007 prices) among the main types of terrestrial biomes (forest, woodland and grassland), almost 50% of the total monetary value estimated for the climate regulation service among all biomes (de Groot et al. 2012). However, this value decreased nearly 29% as a result of LULC changes between 1997 and 2011, mainly due to the loss of tropical forests (nearly −34%), despite the slight increase of temperate and boreal forests (nearly 1.7%) and grasslands and rangelands (nearly 13%) over this period (Costanza et al. 2014).
The assessment of the climate regulation service in a spatially explicit manner is widespread in the scientific literature (e.g. Egoh et al. 2012;Martínez-Harms & Balvanera 2012;Crossman et al. 2013). Also tools combining biophysical quantification with economical valuation into maps, facilitating spatially explicit assessment and modelling, are available for general use (e.g. Nelson et al. 2009;Polasky et al. 2011;Bagstad et al. 2013;Leh et al. 2013). Often, carbon storage and sequestration are used as indicators in the spatial assessment of ecosystem services, namely climate regulation (de Groot et al. 2010;Reyers et al. 2010;Maes et al. 2016). The integration of spatially explicit indicators in the assessment of ecosystem services, as well as the analysis of potential impacts of future land management scenarios (Nelson et al., 2010), is of high importance to increase our understanding of the ecological processes that support the capacity of ecosystems to provide services, both currently and in the future, and to identify and evaluate the benefits that human societies can obtain from them (de Groot et al. 2010;Seppelt et al. 2012). Such indicators and scenarios are also useful to assess key policy questions on the maintenance and restoration of ecosystems and their services, such as the spatially explicit identification of synergies and trade-offs among different ecosystem services, and between ecosystems services and biodiversity (Reyers et al. 2010;Maes et al. 2013) and to better integrate the concept and the value of ecosystem services in landscape planning and management (de Groot et al. 2010).
Various approaches have been used to model and map the supply of climate regulation ecosystem services, the most common and simple being the quantification of carbon stocks associated with soil and vegetation pools, and their variation over space and time due to changes in LULC (Crossman et al. 2013). More complex approaches usually involve processbased models which simulate carbon dynamics in the vegetation (Kim et al. 2015) and in the soil (Smith et al. 1997;Li-Xia & Jian-Jun 2003), considering interactions between soil, climate, vegetation growth, disturbances and/or land-use change. The economic valuation of the service is usually conducted under cost-based methods, often using the avoided damage costs approach (Brouwer et al. 2013, de Groot et al. 2012 to estimate the economic benefits from carbon sequestration and storage. This approach relies on a societal perspective of carbon valuation, instead of a market-based approach, using the social cost of carbon (SCC) as a reference for carbon price (Valatin 2011). SCC estimates the economic cost of an additional ton of carbon dioxide emitted to the atmosphere (Nordhaus 2011) or the monetized benefit and its present value obtained by the sequestration of an additional metric ton of carbon (Nelson et al. 2009).
The assessment of the supply and value of ecosystems services, such as climate regulation and others, in mountain ranges and their ecosystems is of extreme importance since the benefits produced are many and extend beyond their borders (Körner et al. 2005;Fischlin et al. 2007). At the same time, mountain areas are highly sensitive to processes of change, namely climate and LULC change, from both ecological and socioeconomic perspectives (EEA 2010), which in turn create uncertainty on the level and value of ecosystem services these systems will be able to supply in the future. Studies on carbon dynamics at the landscape level are still lacking (Chen et al. 2014), but this is even more evident in mountain areas where detailed evaluations of carbon vertical (pools) and horizontal (landscape) dynamics are currently insufficient. Assessing climate regulation through carbon sequestration and storage in space and time is also important from a planning point of view. The definition of aims and uses of different land types in mountain areas are today possibly changing towards the provisioning of ecosystem services affecting communities at larger scales, namely those related to biodiversity conservation, aesthetics and carbon storage. A thorough evaluation of roles and values of particular land uses is today hard to conduct considering the insufficient current knowledge on the effects of landscape change on ecosystem services in supply and value. In most of the cases, detailed data and/or detailed analysis of changes limit further planning efforts in mountain areas. Reliable indicators for the effects of LULC change on climate regulation services require testing in particular conditions of mountain ranges.
In this study, we focused on the assessment of the impacts of recent past and potential near-future landscape changes on the supply and value of the climate regulation ecosystem service in a mountain area of Portugal. We used as indicators the amount of carbon sequestered and stored in the landscape and its social value. Based on detailed data on carbon ecosystem partitioning and vegetation growth, we aimed to expand our understanding of vertical (among carbon pools) and spatial (among landscape classes) carbon dynamics as a result of LULC changes. We aimed also to understand how different landscape change trajectories may affect the future supply and value of the climate regulation ecosystem service in mountain rural areas, thereby providing background information and knowledge for land planning in changing mountain landscapes.

Study area
The research was conducted in the Sabor River's upper basin (Figure 1), a mountain area approximately 30,650 ha in size located in the district of Bragança, northeastern Portugal (lat 41.9893°to 41.7691°, long −6.5747°to −6.82292°). Current population is 26,000 mostly in the city of Bragança. Despite depopulation and agricultural abandonment observed over the last decades in rural parishes of the study area Pinheiro et al. 2014), the main socio-economical activities in the area are still related to agriculture and forestry even if major sources of income of land owners are in the services sector in the city of Bragança.
The Sabor River's upper basin is part of the Montesinho Natural Park and of the Natura 2000 network under the Birds and the Habitats EU Directives. The area covers a wide range of topographic and climatic conditions: a mountainous area (Serra de Montesinho) with the highest elevation at 1486 m (average annual temperature 8.5°C; average annual total precipitation above 1200 mm) in the west; a plateau on average at 900 m (average annual temperatures from 10°C to 12.5°C; average annual total precipitation from 800 to 1000 mm) in the east; and the valleys of the Sabor River and its tributaries at elevations lower than 650 m (average annual temperature from 12.5°C to 14°C; average annual total precipitation from 800 to 1000 mm) in the centre (Aguiar 2000;Pereira 2006;IPB/ICN 2007). The most representative soil types are umbric leptosols associated to granite and schist lithology, and dystric leptosols. Other soil units are eutric leptosols associated with igneous basic and ultrabasic lithology, umbric and dystric cambisols, chromic luvisols and fluvisols associated with igneous basic lithology and haplic alisols associated with sedimentary lithology (Agroconsultores&Coba 1991;Aguiar 2000;Pereira 2006;IPB/ICN 2007).
The diversity of topographic, climatic and geological conditions and the diversity of human activities in the area are reflected in the diversity of vegetation and land use/cover types across the landscape. Seminatural habitats represent almost two-thirds of the area, comprising a high diversity of scrub communities (dominated by Cistus spp., Cytisus spp. and Erica spp.) and open forests. Forest areas comprise natural forests of Quercus pyrenaica (deciduous) and Quercus rotundifolia (evergreen), riparian forests of Alnus glutinosa and Fraxinus angustifolia, and planted Pinus pinaster stands. Despite a recent strong decrease in agriculture, agroforestry systems of Castanea sativa, temporary and permanent crops, meadows and pastures are still well represented in the area.

Research framework
We assessed the effects of landscape change on the climate regulation ecosystem service in the study area through the quantification, valuation and mapping of carbon sequestration and storage, considering (1) the LULC changes that took place between 1990 and 2006, and (2) the changes likely to take place in the near future as a result of socio-economic drivers based on three alternative landscape scenarios projected for 2020. We used the InVEST (Integrated Valuation of Ecosystem Services and Tradeoffs) model (Tallis et al. 2013;Sharp et al. 2015) for scenario building and for carbon sequestration and storage assessment and valuation. The model was fed with LULC spatial data for the three dates under analysis and with carbon data on live above-and belowground biomass, litter and soil organic carbon (SOC) for each of the major LULC classes. Carbon data for each of these pools were obtained from data collection and/or modelling.
The InVEST carbon module compares levels of carbon in those four different pools over time based on LULC spatial data. The application of the InVEST carbon module results therefore in a three-dimensional comparison of stored carbon in the landscape, and the dynamics of which is driven by LULC change. From that comparison, sequestration rates (SRs) and the value of the climate regulation ecosystem service are derived. Although processes responsible for carbon sequestration and release fluxes are not explicitly incorporated in the module, it can provide reliable estimates of carbon stored in the landscape and of the role of landscapes and particular LULC classes to supply the ecosystem service under consideration. This approach has been widely used in research (see the review by Crossman et al. 2013). Uncertainty could not be addressed in the physical assessment of carbon dynamics in the landscape. The results of the analyses refer, therefore, to average input conditions for the pools considered.

LULC cartographic basis and scenarios
We used LULC spatial databases of the Sabor River's upper basin for years 1990 and 2006, created based on the methodological approach described in Carrão et al. (2008) and based on interpretation of orthophoto maps (Amorim 2015) ( Figure 1). Additionally, three alternative landscape scenarios projected for 2020 were developed using the InVEST Scenario Generator module (Table 1). These scenarios aimed to understand how different potential changes in the landscape may impact the future supply and monetary value of the climate regulation ecosystem service. The three scenarios were based on observed and potential trends of LULC changes for the study area described in Azevedo et al. (2011) and Pinheiro et al. (2014) (Table 2).
• Agriculture abandonment (Abandonment): simulates a moderate management scenario, where the LULC change trend observed in the recent past (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006) will persist in the future; reflects a decrease of agricultural land use due to depopulation and abandonment of rural areas; agricultural land tends to be replaced by shrubland, natural regeneration forest and/or forest plantations. • Forest expansion (Forest): simulates a more intensive management scenario, including investments in nature conservation and carbon emissions mitigation measures such as reforestation, afforestation and disturbances management (e.g. wildfires); forests tend to increase at the expense of abandoned agriculture fields and shrublands. • Shrubland expansion (Shrubland): simulates a less intensive management scenario, where the effects of agricultural abandonment and the increase of wildfires will result in the decrease of agriculture and forest areas that will tend to be replaced by shrublands.
Furthermore, we added information on physical and environmental factors (e.g. elevation, aspect, distance to rivers, distance to urban areas and fire risk cartography), which was combined with the transition likelihood matrix in the InVEST Scenario Generator module to determine the spatial suitability for the allocation of LULC changes in the landscape (see Appendix A for details).

Vegetation and litter carbon
To estimate the amount of carbon stored in aboveand belowground biomass and litter, we used previously published data describing quantitatively biomass and carbon storage at local, regional and national scales, but for land uses with higher levels of carbon content, such as forests, we used biomass equations and growth and yield models (Appendix B).
In the case of forests, we used data for stand variables (e.g. age, dbh, dominant height and tree density) available for the study area (Pires 1998;Xavier 1999). The estimation of above-ground biomass and carbon was done using forest growth and yield functions and models (e.g. Pérez-Rodríguez et al. 2016). Since age was known for most of the stands and since tree growth and stand biomass strongly affects the rate and amount of carbon sequestered and stored in the landscape, particularly at young ages (such as those observed in the study area), we included growth throughout the period of analysis at both above-and belowground levels. This procedure, however, has only been carried out in forests, while for the remaining LULC classes we assumed that biomass was constant over time. Belowground biomass and carbon were estimated using mostly the allometric functions of Montero et al. (2005). Biomass was converted into carbon content using conversion factors as described in Appendix B. Average carbon stocks for above-and belowground biomass and for litter were calculated for each major LULC class. These data were later used to build the carbon transition matrices in InVEST.

Soil organic carbon
Spatial sampling design and field surveys We used a two-stage sampling design (Gruijter et al. 2006) to select locations for SOC data collection. In the first stage, we followed a stratified random sampling approach to selected 25 Primary Units (PU) registered in a 1 × 1 km regular grid in the WGS 1984 UTM 29 N coordinate system. The stratification layer included spatial data related to climate, topography, soil types, fire regime and the distribution of areas with different nature protection status. Stratification was obtained by the Hard Competitive Learning clustering algorithm (Rolls & Deco 2002) with a total of five strata implemented in the R statistical software (R Development Core Team 2014). To reduce surveying costs in the entire PU area, in a second stage we used a systematic sampling approach to select five 200 × 200 m Secondary Units (SU) per PU, located at the corners and the centre of each PU.
From the 125 SUs selected for field surveying, soil samples were collected in 120 of these SUs since one PU was of very difficult access. At the centre of each SUs, one soil core was collected in the 0-5, 5-10, 10-20 and 20-30 cm layers. Soil sampling was only possible in the 0-30 cm depth interval because soils in the surveyed sites were usually very shallow and stony. However, the mineral soil surface layer is the main volume relevant for spatial survey of carbon storage, according to the Kyoto Protocol requirements (Smith et al. 2002;Schulp et al. 2008;Vesterdal et al. 2008). To determine bulk density, undisturbed soil samples were taken using a core sampler with a volume of 100 cm 3 .

Lab determination
Total SOC (down to 30 cm depth) in each plot was calculated from carbon concentration and soil bulk density at the same depth. Samples for soil carbon evaluation were air-dried and sieved to determine the coarse fraction (>2 mm) and analysed for carbon concentration by dry combustion (ISO 1995). SOC content (kg C m −2 ) was calculated for the mineral soil layer using Equation (1) (Tate et al. 1997;Percival et al. 2000): where CC is the carbon concentration of the mineral soil layer: CC (g kg −1 ) equivalent to CC (kg t −1 ), BD is the bulk density of the mineral soil layer: z is the thickness of the mineral soil layer: z (cm)/100 equivalent to z (m) and CE is a correction factor for coarse elements content Total soil carbon storage per unit area (kg C m −2 ) was finally estimated by summing the mean SOC values obtained for each of the three soil layers.

Carbon modelling
Estimates of SOC stocks for particular combinations of major LULC class and other land variables were calculated based on an exploratory cluster analysis performed with data from the 120 samples collected in the field and data on 'altitude'. The selection of the variable 'altitude' was based on previous work on SOC modeling at the landscape level carried out in the study area, where it was found to be a significant predictor of SOC in the landscape  Before the cluster analysis was carried out, the 120 samples were aggregated into four major LULC classes, namely agricultural areas (n = 9), seminatural areas (n = 81), coniferous forest (n = 11) and broadleaved forest (n = 19). Then, the Ward's hierarchical clustering method was applied with the squared Euclidean distances for the variable 'altitude', followed by a k-means cluster analysis based on the number of clusters identified as meaningful in the dendrogram of the hierarchical analysis conducted in each major LULC class. This process allowed the discretization of seminatural systems and broadleaved forests in altitudinal carbon clusters. Samples within each of these two major LULC class were then classified according to the k-means clustering results. SOC values were later averaged for each group within major LULC classes as shown in Table 3. In coniferous forests, however, clustering was not observed and SOC was grouped according to the forest-type information, namely into the categories 'maritime pine stands', 'maritime pine plantations' and 'other coniferous' (Table 3). In the case of agricultural area, due to low variability of agricultural cover types, we assumed an average value of 24.52 Mg ha −1 throughout the study area.
With the purpose of avoiding abrupt and unlikely changes of SOC stocks in LULC transitions, we considered a gradual change of SOC stocks between LULC classes according to the soil carbon SRs and temporal patterns for LULC conversions described in Deng et al. (2016) (Table 4) applying the following equation: where SOC f is the SOC stock of the new LULC (Mg ha −1 ), SOC i is the carbon stock of the LULC converted (Mg ha −1 ), SR is the soil carbon SR corresponding to the major LULC transitions (Mg ha −1 yr −1 ) listed in Table 4 (adapted from Deng et al. 2016), and ΔT is the time period of the LULC transition (yrs).
In the cases where forests and seminatural areas did not change over time, we assumed that carbon in the soil varied depending on a factor calculated based on a 50-yr period, which represents the annual increment of soil carbon. This rate was used to add or to discount carbon to SOC values previously estimated. However, for agricultural and agroforestry systems, we assumed soil carbon to be constant over time in the absence of LULC transition considering that the effects of tillage balanced the biomass incorporation into the soil.

Invest input data preparation
Simulations of carbon dynamics with the InVEST carbon module require for the same area pairs of LULC maps of consecutive dates and transition matrices for the four carbon pools between these dates. The establishment of the LULC classification system used to guide the construction of input maps and the corresponding carbon pools matrices was a complex process involving several steps. First, we created maps of the area classified according to the six major LULC classes considered in the preparation section LULC cartographic basis and scenarios Table 1. The SOC pattern observed and assessed through cluster analysis forced the LULC classification to include extra classes to take into account observed differences in soil carbon due to elevation and age of forest stands, resulting in 11 different classes (Table 3). In a subsequent step, with the information of Pires (1998) . It was, however, the best way of dealing with the specific land-use transitions in the areas in a detailed manner taking into account the variables that were found relevant to explain the distribution of SOC in the study area and the variables directly related to the above-ground standing biomass. Based on these landscape components, we recoded the InVEST carbon pools matrices in order to incorporate not only possible LULC transitions but also the actual condition in any evaluation period of all pools. All maps were prepared with the 25-m spatial resolution allowed by the underlying land cover maps.

Economical valuation
We used the avoidance damages costs approach (Pascual et al. 2010) taking the SCC prices as a proxy of the economic damage avoided, or the monetary benefits that society can obtain, due to carbon emission reduction, conventionally related to an additional metric ton of carbon sequestered from the atmosphere (Nelson et al. 2009).
Given the uncertainties about the appropriate SCC value and discount rate to value the climate regulation ecosystem service (Conte et al. 2011), we carried out a sensitivity analysis (Boardman et al. 2001) applying three SCC values, $23/Mg C, corresponding to the estimated SCC mean value from a meta-analysis of 211 estimates assuming a conservative assumption on greenhouse gas emission reduction policy (Tol 2008), $44/Mg C, corresponding to the estimated SCC mean value assuming an optimal economical and climate trajectory according to a limited global temperature rise of 2°C above the 1900 average (Nordhaus 2011), and $312/Mg C, corresponding to the estimated SCC mean value assuming a 'business as usual' trajectory (Stern 2007), and three discount rates, 1.4% (Stern 2007), 3.5% (Valatin 2011) and 7% (Sharp et al. 2015). Due to constraints of the InVEST carbon module, we assumed in the simulations the nearest integer value of the discount rates values mentioned earlier (1%, 3% and 7%). The annual rate of change in the price of carbon was established as 5% (Nelson et al. 2009) for all the simulations. A total of 36 (3 discount rates × 3 SCC prices × 4 scenarios) simulations were performed in the InVEST carbon model in order to cover all possible combinations of LULC scenarios and economic parameters.

Results
Recent carbon storage and sequestration (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006) Carbon stored in the study area increased in 49.20% over a 16-yr period (Figure 2, Figure 3), corresponding to an SR of 1.45 Mg C ha −1 yr −1 (Figure 4). In 1990 and 2006, soil was the most important carbon pool in the landscape, increasing its storage capacity during this period and was responsible for maintaining near 63% of all carbon in the landscape in both dates ( Figure 2). All carbon pools showed a positive carbon SR during this period. Soils presented the highest rate, followed by above-and belowground biomass, and finally by litter (Figure 4). Carbon density was also higher in soils (approximately three times higher than in above-ground vegetation) in both 1990 and 2006 (Table 5).
Carbon was stored mainly in seminatural systems, both in 1990 and 2006 (Table 6). Although these areas increased their storage capacity in the landscape, the relative contribution of seminatural areas for total carbon in the landscape in 2006 decreased mainly due to an increase of carbon storage in forests areas (Tables 6 and 7). Carbon stored in agriculture areas was reduced during this period, which resulted in carbon emissions over time (Tables 6 and 7).
In 1990 and 2006, agriculture and seminatural systems stored carbon mainly in the soil (Table 6). While in agriculture areas carbon stored in the soil decreased slightly, it increased in seminatural areas. Moreover, soils in seminatural areas were the only carbon pool sequestering carbon in this period, presenting the highest SR of all carbon pools in all major LULC classes, which balanced the carbon losses in these systems derived from carbon emission in vegetation pools (Table 7).
Contrarily, in forest areas, carbon was stored mainly in the vegetation pools, especially above ground. Although all pools in forest areas gained carbon, above-ground biomass and soil pools showed the highest rates (Tables 6 and 7). The contribution of the vegetation pools to the overall carbon in forest areas, however, decreased slightly due to a decrease in belowground biomass, while the contribution of the soil carbon pool increased (Table 6).
In the forest and seminatural LULC classes, carbon density increased considerably from 1990 to 2006, but decreased slightly for agriculture (Table 8). Higher carbon densities in agriculture and seminatural classes were found in soil pool, whereas in forests the highest density was found in the above-ground biomass pool.

Carbon sequestration and storage under future LULC change scenarios (2006-2020)
All three LULC change scenarios would increase the total carbon stored in the study area comparatively to the 2006 baseline conditions (Figures 2 and 3). Despite the similarities observed among scenarios, the forest expansion scenario showed an SR above the rate observed for the 1990-2006 period ( Figure 3).
The overall distribution of stored carbon across pools in the 2020 scenarios followed the patterns observed for 1990 and 2006: carbon is stored manly in the soil, followed by the above-and belowground biomass, and litter pools (Figure 2).
In the abandonment and shrubland scenarios, carbon in the landscape was stored mainly in seminatural systems, followed by forests and agricultural areas, in a pattern similar to that observed in 1990 and 2006 (Table 6). Contrarily, in the forest expansion scenario, forests stored most of the carbon in the landscape, followed by seminatural areas and agricultural areas (Table 6).
Forest areas gained carbon in all scenarios and in the majority of carbon pools, especially the soil and the above-ground biomass carbon pools in the forest scenario, losing carbon in soil and litter pools in the shrubland scenario (Table 7). Carbon in seminatural   areas increased in the abandonment and shrubland scenarios, but decreased in the forest scenario (Table 7). Moreover, all carbon pools of seminatural areas gained carbon in the abandonment and shrubland scenarios, while in the forest scenario only the soil carbon pool gained carbon (Table 7). Finally, agriculture areas emitted carbon in all scenarios, especially from the soil pool in the abandonment scenario (Table 7). Forests and seminatural areas increased their carbon content in all scenarios (Table 8), whereas agricultural areas maintained the 2006 levels in all scenarios.

Spatial distribution of carbon storage and sequestration
The maps generated from the InVEST carbon module revealed the impacts of LULC changes on carbon stocks across the landscape, allowing the identification of areas that sequestered or emitted carbon over time (Figures 5 and 6). Moreover, these maps allowed to identify areas that can potentially supply this ecosystem service and areas of net carbon emissions where the supply of the service is uncertain according to future LULC scenarios.
In general, between 1990 and 2006, higher levels of carbon sequestration were observed mostly in the highest elevation areas (above 900 m) located in the northwest of the basin, where seminatural areas dominate, and in the northeast part of the basin (elevation from 650 to 900 m), an area dominated mostly by forest areas and seminatural areas. Lower levels of sequestered carbon, or even carbon emissions, were observed in lower elevation areas (below 650 m) corresponding to the central part of the basin, dominated mainly by agriculture and urban areas.
The pattern observed for the 1990-2006 period was also, in general, observed for 2020 in the three scenarios ( Figures 5 and 6). In addition, the forest expansion and the agriculture abandonment scenarios showed an increase in carbon sequestration in riparian areas. Contrarily, in the shrubland expansion scenario, although carbon sequestration areas were spatially distributed similarly to 1990-2006 patterns, there were also carbon emission in areas in the northern part of the basin where forests were predicted to be replaced according to this scenario.
Over the simulated period (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020), the monetary value of sequestrated carbon also increased significantly (Table 9) although there were differences among scenarios. The forest expansion scenario presented the highest increase in value, followed by abandonment and shrubland expansion. As observed for the 1990-2006 period, a low societal preference for immediate benefits of carbon sequestration lead to higher monetary values of the service, decreasing as the discount rate increased (Table 9). Considering the entire 30-yr period, the total economic costs avoided due to carbon sequestration were higher in the landscape moving towards the forest scenario. For the worst-case scenario (1% discount rate, SCC of $23/Mg C) this value was $21.22 million ($23.08 ha −1 yr −1 ). For the shrubland scenario, the minimum total economic costs avoided was $18.72 million ($20.37 ha −1 yr −1 ) ( Table 9).

Landscape change and carbon dynamics
The total carbon SR in our study (see Figure 4) is comparable with others reported for similar conditions. Padilla et al. (2010), for example, estimated a rate of 1.27 Mg C ha −1 yr −1 based on LULC changes over a 74-yr time period in a Mediterranean mountain area in South-East Spain. The similarities with this case can be explained by an identical LULC change trend (expansion of forest areas and contraction of agriculture and shrubland areas).
Considering just vegetation, however, our results were higher than those estimated by Pinheiro et al. (2014) in an area adjacent to our study area based in changes over a period of 48 yrs (0.272 Mg C ha −1 yr −1 ). Our results are much higher than those obtained by Muñoz-Rojas et al. (2011) in an analysis of LULC changes over 51 yrs in southern Spain (0.039 Mg C ha −1 yr −1 ) and those of Tappeiner et al. (2008) covering 140 yrs of LULC changes in an Alpine valley area in Switzerland (0.04 Mg C ha −1 yr −1 ), which may reflect differences in cover types present in each area, as well as in types of changes observed over time, and also the different  environmental and geographical context that can influence the carbon SRs in each vegetation cover type.
The estimates we obtained for total landscape carbon density in vegetation (see Table 5   scenarios and the results of Muñoz-Rojas et al. (2015) can be explained, at least partially, by the expansion in forest areas (forest scenario), and forest and shrubland areas (abandonment scenario) that was higher in our case than in southern Spain where, despite the increase of forest areas, agriculture intensification was also observed, which may have decreased the total landscape soil carbon density over time.
Rates of carbon sequestration in forest biomass in our study (see Table 7) revealed to be low when compared with other estimates for Portugal. Pereira et al. (2009), for example, reported rates from 4.09 to 7.09 Mg C ha −1 yr −1 for maritime pine, 4.09 to 8.73 Mg C ha −1 yr −1 for eucalyptus and 0.54 to 2.18 Mg C ha −1 yr −1 for Pyrenean oak forests. These rates reflect the productivity of the regions where data were collected and the geographic scale at which they are reported (the cases of maritime pine and eucalyptus) or the carbon pools under consideration (the case of Pyrenean oak that includes SOC). Ribeiro et al. (2011) estimated carbon sequestration to occur at an average rate of 1.13 Mg C ha −1 yr −1 in the forests of all northern Portugal, including the coastal region. In a mountain area in southeastern Spain, Padilla et al. (2010) reported an SR of 2.7 Mg C ha −1 yr −1 . Our results, however, reveal that forest growth rates in this study region are much lower than regions under other ecological conditions. Using National Forest Inventory data for the entire northeastern region of Portugal, where our study area is located, Azevedo (2012)  The improvement of the climate regulation ecosystem service in the Sabor River's upper basin from 1990 to 2006 was due to changes in vegetation (see Table 1) and in soil carbon storage, the largest carbon pool across the landscape, and to the spatial pattern of soil carbon and LULC changes (see Figure 2). Changes in LULC involved the reduction of low carbon content agriculture crops and the expansion of forests and other tree-dominated systems of higher carbon density (see Table 1). On the other hand, there is a spatial pattern of SOC strongly related to elevation (see section on "Soil organic carbon"). Indeed, soils at higher elevations have much higher carbon contents than soils at lower elevations. In addition, shrubland and forest LULC classes usually add carbon to the soil (see section on "Carbon modelling"). The expansion of LULC classes that combine higher carbon content in vegetation pools and highest SOC levels at higher elevation increases carbon storage in most of mountainous areas of the study area.
In general, negative balances of carbon sequestration and storage (net carbon emissions) occurred both in soil and vegetation pools when LULC changed, e.g. when seminatural or forest areas that were converted to agriculture areas ( Figure 6). Although losses in stored carbon derived from these LULC conversions were balanced in the overall landscape carbon stocks, it is of utmost importance to take these conversions into account in landscape change scenarios, since they can affect negatively the provision of the climate regulation service (Zhu et al. 2010;Scharlemann et al. 2014).

Economical valuation
Besides contributing to avoid economic damages for global society due to climate change, as a climate regulation service, carbon sequestration and storage Table 9. Total economic costs avoided (Million $) and rate ($ ha −1 yr −1 ) due to carbon sequestration among LULC scenarios per SCC price ($ (Mg C) −1 ) and discount rate. in the Upper Sabor landscape adds value to natural resources in this mountain area. The results obtained in our study (see Table 9) are similar to those obtained from Padilla et al. (2010) despite the different carbon price used ($67/Mg C). These authors reported values that ranged from a minimum of $23.09 ha −1 yr −1 , in a landscape scenario where forest areas were converted to shrubland/grassland, to a maximum of $118 ha −1 yr −1 , when shrubland/grassland areas were converted back to forest, and an intermediate value of $85.16 ha −1 yr −1 corresponding to the abandonment of agriculture and shrubland areas and expansion of forest areas.
In their meta-analysis, Brenner et al. (2010) estimated that the gas/climate regulation service of temperate forests and grasslands in the coastal zone of Catalonia (Spain) could generate $140 ha −1 yr −1 , which is a higher value than our estimates for an SCC = $23 and $44/Mg C, but lower when we use an SCC = $312/Mg C. Mendes (2005) estimated that Portuguese forests (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) could generate a monetary value of $0.80 ha −1 yr −1 (for a carbon price of $18.51/Mg C), which is very low compared to our estimates, even using a carbon price similar to our lower SCC value. Finally, Moore et al. (2011) based on a value transfer analysis, estimated that private non-urban forests in Georgia (USA) can generate $11.33 ha −1 yr −1 , using an average carbon price of $21/Mg C and a 7% discount rate, a value slightly lower than ours for a similar carbon price (SCC = $23/Mg C) and discount rate (7%).
The differences shown between our results and those from the studies referred earlier reflect the uncertainty about the economical value of this ecosystem service, which mainly derives from different carbon prices and emission trajectories used in the valuation, as well as the uncertainties about future society preferences regarding the climate mitigation actions. Nevertheless, the assessment of the economical value of this and other ecosystem services will be an important tool to strengthen the value of natural resources and to better support planning and management actions.

Insights for management and planning
Although based on a deterministic and non-mechanistic methodological approach, this study has provided evidence of the importance of landscape change as a major process driving spatiotemporal carbon dynamics as well as the provision of the climate regulation ecosystem service in a mountain rural landscape. The study has also confirmed that there is a considerable potential monetary value associated to this ecosystem service in the study area, and that future landscape changes will affect that value. This ensemble of results can be valuable to inform decision-making regarding planning and management in this and other similar mountain landscapes.
One of the major implications of our results is in terms of the role that mountain areas should play at regional and national levels within a climate change mitigation policy framework (Kohler & Maselli 2009). The replacement of agriculture by forest and shrubland LULC classes in higher elevation areas increase soil carbon sequestration at rates higher than in lowland areas, besides adding carbon in above-and belowground biomass. Due to high precipitation, but mostly to low temperatures during most of the year, sites at higher elevations retain higher levels of undecomposed organic matter. For this reason, soils in mountain areas are important carbon sinks (Aguiar et al. 2009), and land management in these areas should take that into consideration to decrease carbon emissions. Regional and national planning frameworks should therefore promote the establishment and persistence of LULC types that not just maintain existing carbon stocks but also promote higher sequestration at the vegetation and soil levels. Policy and planning in forestry, biodiversity conservation (Protected Areas and Natura 2000 network sites, often located in mountains), game management, fire risk management, rural development and other sectors in mountain areas thus require the incorporation of carbon sequestration and storage objectives and targets, for which our results can be of great importance.
Seminatural areas (shrubland) have a particular role in higher elevations, where they are dominant across the landscape and in practice the only alternative since forest plantations are difficult and expensive to establish. Shrublands have, however, a negative connotation among local communities, including land owners and managers, due to their low productivity and high fire-proneness, and due to the low diversity of income generation activities that can be based on these ecosystems. The evidence provided in this research that shrublands are able to maintain high carbon densities, especially in the soil pool, contributes significantly to the existing knowledge on shrub formations in terms of sequestration and storage (Fonseca et al. 2011(Fonseca et al. , 2012Nogueira et al. 2014), and supports the integration of these areas in carbon-related policy and planning. In addition, the economic valuation of the climate regulation ecosystem services supplied at higher elevations, mostly provided by these seminatural systems in the study area, offers relevant quantitative information that can be used in land planning in the region and in the country. It can also be used as a way of valorization of these systems that tend to be overlooked locally.
In addition to the changes observed in the recent past decades, there is still potential to further rise the supply of the climate regulation ecosystem service through carbon sequestration and storage. Any of the scenarios analysed in this study predicts increasing carbon storage in the landscape. Also the monetary value of the landscape from this point of view increases considerably, at rates similar to those observed for the 1990-2006 period. Although the results of the projections for 2020 suggest that naturally carbon in the landscape and in particular systems will increase in the future, the implementation of soil conservation measures should be considered a priority given the role that soils play in sequestration and storage of carbon in mountainous areas. These measures are particularly relevant in fire-prone landscapes, such as in this study area, where wildfires can lead to direct carbon emissions by wildfires and subsequent erosion processes (Hurteau & Brooks 2011;Vågen & Winowiecki 2013;Scharlemann et al. 2014).
Our results on carbon storage and sequestration in forest and seminatural areas were in part due to the maintenance and growth of these units over time, and due to the disturbance regime in the area that tended to be of small scale even if sometimes locally frequent. This is related to decreasing pressure of grazing in most of the area and to the infrequent occurrence of wildfires comparatively to other areas in Portugal. Maintaining these low disturbance regimes will have also a strong contribution to the maintenance of the carbon storage capacity of the landscape. Although fuel buildup in the area increases fire hazard , management practices that can reduce fuel loads and break the spatial continuity of highly flammable areas will have a strong positive contribution to resilience to fire in the area (Fernandes 2013). In fact, prescribed burning or mechanical removal of vegetation, although affect locally above-ground biomass, will preserve belowground biomass and SOC, decreasing the susceptibility of soils to erosion as well.
Our results incorporate fire disturbance for the 1990-2006 period and for one of the future scenarios (shrubland). Future shifts in fire frequency/intensity from the levels assumed in this work, particularly at higher elevations, may affect the performance of particular LULC classes and pools in carbon sequestration and storage. This outcome is however still uncertain. Similarly, changes in climate, not addressed in this research, are likely to affect the net carbon balance of many LULC classes and of the landscape as a whole in an uncertain way. Although an increase in above-and belowground biomass is likely to occur due to longer growth seasons and higher temperatures, it is also expected that carbon emission from the soil pool will increase, thus decreasing net carbon sequestration. The potential magnitude of these effects is still unclear. These effects should be higher in elevated areas where currently dominant ecosystems such as shrublands tend to play an important role as carbon sinks (Table 8, Figures 5 and 6). Interactions between climate change, vegetation and fire should also be considered from a practical perspective. Warmer, dryer and longer summers combined with higher vegetation growth rates may increase fire hazard in the region, therefore potentially increasing overall emissions. The uncertainty around the types and magnitude of these and other interactions is, however, very high.
The spatially explicit assessment of the past and future trends of this ecosystem service can contribute to a better identification of the most suitable areas for the provision of this regulation service as well as to the quantification of their value. It can also allow the identification of sensitive areas where carbon emission takes place or where it is likely to occur in the future due to LULC change or other driver. Considering that these areas correspond to a certain pattern of previous LULC and soil carbon, this can inform planning and management of this and other areas.
Although the benefits and the beneficiaries of the climate regulation service are global, carbon sequestration and storage in ecosystems and in the landscape must be considered in planning for specific mountain areas. This would allow supporting and exploring the development of schemes of payment of ecosystem services in the area as a way of valorization of mountains areas by the services they provide to the society as a whole.
Finally, with the efforts to increase knowledge on carbon dynamics and on the provision of the climate regulation ecosystem service at the landscape scale, we consider that continuous monitoring actions at this scale level are necessary, both on carbon pools dynamics and on land-use patterns, in order to validate model predictions, reduce uncertainty and increase the ability of an ecosystem services-based approach to provide better support in planning and management actions in this and similar areas.

Conclusion
This study assessed the provision of climate regulation ecosystem service, both biophysically and economically, through the analysis of the carbon storage and sequestration dynamics at the landscape level as a result of LULC changes occurred between 1990 and 2006 and in three scenarios for 2020 in a mountain area in northeastern of Portugal.
In general, recent past LULC changes in the study area had a positive effect on carbon storage and sequestration, increasing the total carbon storage in 49.20%, at a rate of 1.45 Mg C ha − 1 yr −1 , mostly due to the expansion of forest areas at the cost of areas of low carbon density such as agriculture areas, increasing the storage capacity among pools, and due to the high levels of carbon stored in soils, that maintained 63% of all carbon in the landscape, especially in high elevation shrubland communities. Hence, conservation actions on forest areas as well as in seminatural systems, especially in high elevation shrubland areas, should be simultaneously considered and balanced.
Besides the evidences of the positive impact that LULC changes had on carbon dynamics, those were also responsible for a considerable monetary value associated with this ecosystem service that ranged from a maximum of $148.26 million ($4837.95 ha −1 ), assuming the highest carbon price and a lower societal preference for immediate benefits over future benefits, to a minimum of $7.86 million ($256.63 ha −1 ), assuming the lowest carbon price and a higher societal preference for immediate benefits over future benefits. Despite the uncertainties associated with the proper carbon price and societal choices, these results can also contribute to reinforce the natural value of the mountain ecosystems.
Furthermore, although the uncertainty associated with the supply of the climate regulation ecosystem service in the future was evidenced by the three alternative landscape scenarios projected for 2020, the potential for carbon storage and sequestration in the study area seems to be high. However, the forest expansion scenario seems to be the most sustainable scenario for the maintenance of this ecosystem service in the landscape, with an SR of 1.49 Mg C ha −1 yr −1 . Additionally, the Sabor River's upper basin can continue to contribute to the avoidance of economical cost associated to carbon emissions with values ranging from a maximum of $287.82 million ($313.06 ha −1 yr −1 ) in the best case and a minimum of $21.22 million ($23.08 ha −1 yr −1 ) in the worst case.
The estimates on carbon stocks and SRs obtained in this research, associated with a spatially explicit analysis, provide baseline information for future studies on carbon dynamic at the landscape level and mitigation scenarios associated to LULC changes, as well as for better supporting future planning and management actions carried out in this and other mountain areas that experienced similar LULC changes. Moreover, the estimates obtained for the economic value of this ecosystem service can also be useful to support schemes for payment of ecosystem services improving local economies in these regions.

Disclosure statement
No potential conflict of interest was reported by the authors.