Comparison of remote sensing and plant trait-based modelling to predict ecosystem services in subalpine grasslands

There is a growing demand for spatially explicit assessment of multiple ecosystem services (ES) and remote sensing (RS) can provide valuable data to meet this challenge. In this study, located in the Central French Alps, we used high spatial and spectral resolution RS images to assess multiple ES based on underpinning ecosystem properties (EP) of subalpine grasslands. We estimated five EP (green biomass, litter mass, crude protein content, species diversity and soil carbon content) from RS data using empirical RS methods and maps of ES were calculated as simple linear combinations of EP. Additionally, the RSbased results were compared with results of a plant trait-based statistical modelling approach that predicted EP and ES from land use, abiotic and plant trait data (modelling approach). The comparison between the RS and the modelling approaches showed that RS-based results provided better insight into the fine-grained spatial distribution of EP and thereby ES, whereas the modelling approach reflected the land use signal that underpinned trait-based models of EP. The spatial agreement between the two approaches at a 20-m resolution varied between 16 and 22% for individual EP, but for the total ecosystem service supply it was only 7%. Furthermore, the modelling approach identified the alpine grazed meadows land use class as areas with high values of multiple ES (hot spots) and mown-grazed permanent meadows as areas with low values and only few ES (cold spots). Whereas the RS-based hot spots were a small subset of those predicted by the modelling approach, cold spots were rather scattered, small patches with limited overlap with the modelling results. Despite limitations associated with timing of assessment campaigns and field data requirements, RS offers valuable data for spatially continuous mapping of EP and can thus supply RS-based proxies of ES. Although the RS approach was applied to a limited area and for one type of ecosystem, we believe that the broader availability of high fidelity airborne and satellite RS data will promote RS-based assessment of ES to larger areas and other ecosystems.


INTRODUCTION
Human societies benefit from a multitude of resources and processes that are supplied by ecosystems.They provide a vast range of services such as food, timber or clean water production, they regulate climate and diminish natural hazards, and they offer nonmaterial cultural assets (Costanza et al. 1997, de Groot et al. 2002, MEA 2005, Burkhard et al. 2009).The concept of ecosystem services (ES) presents a way to quantify, analyze and manage the benefits obtained from ecosystems (MEA 2005).Although the interest in ES has grown exponentially since the 1990s (Costanza et al. 1997, de Groot et al. 2002, MEA 2005, Fisher et al. 2009), classification and consistent quantification of ES still remain a challenge (Wallace 2007, Fisher et al. 2009, de Groot et al. 2010, Martı ´nez-Harms and Balvanera, 2012).Nevertheless, recent studies indicate wide international efforts towards a common, global system for monitoring of ES (Tallis et al. 2012, Haines-Young andPotschin 2013).Spatially explicit mapping of ES at different scales is required for sustainable land use planning and environmental decision making (Burkhard et al. 2012).Moreover, there is a growing interest in mapping of multiple ES and identification of areas with concentrated ES supply (Naidoo et al. 2008, O'Farrell et al. 2010, Lavorel et al. 2011).
The biggest challenge in spatially explicit mapping of ES is often a limited availability of primary data on ES (Eigenbrod et al. 2010).It is thus more common to use ES proxies-ecosystem properties, environmental variables or land cover / land use maps-which are naturally linked to ES and, at the same time, easier to obtain (Chan et al. 2006, Egoh et al. 2008, Eigenbrod et al. 2010, Lavorel et al. 2011).The mapping of ES is often based on land use or land cover data by transferring a single value of ES to each class (Costanza et al. 1997, Metzger et al. 2006, Burkhard et al. 2009, Eigenbrod et al. 2010).This approach is practical, but also less accurate as it neglects variability within a single class.Eigenbrod et al. (2010) showed that land cover based proxies of ES poorly correlated to primary data of biodiversity, recreation and carbon storage ES in England.At the other end of a gradient in model complexity with respect to biodiversity representation (Martı ´nez-Harms and Balvanera 2012), Lavorel et al. (2011) showed that ecosystem properties (EP), underpinning ES supply, can be better estimated when including spatial variations in environmental variables and plant traits than using a pure land use based model only.This spatially explicit approach required intensive field measurements to calibrate linear regression models of plant traits and EP.Considering large variability of vegetation properties within land use classes (Garnier et al. 2007) and across landscapes due to environmental gradients (Albert et al. 2010), field data can be a limiting factor in EP and ES mapping at larger spatial scales.Remote sensing (RS) offers a possibility to map spatial and, for some type of RS observations, temporal variation of ecosystem properties in a more extensive manner than individual field measurements (Kokaly et al. 2009, Ustin and Gamon 2010, Mulder et al. 2011, Alcaraz-Segura et al. 2013).
An advantage of the current RS technology is that it supplies a large variety of spectroradiometers, which acquire suitable data for ecosystem mapping at different spatial, spectral and temporal resolutions (Turner et al. 2003, Grace et al. 2007, Malenovský et al. 2009).Specifically, Ayanu et al. (2012) reviewed the suitability of existing satellite-based spectroradiometers and RS methods for ES assessment and concluded that more studies are needed to evaluate the validity of RSbased proxies for ES mapping.A disadvantage, mainly for airborne based spectroradiometers, can be higher operational costs as compared to field sampling of EP.Nevertheless, airborne RS data are prerequisite for development and testing of new RS-based methodologies, as demonstrated in this study.
In general, RS data can be used in two ways to support spatial assessment of ES.First, RS data serve to generate accurate and up-to-date land cover maps (Friedl et al. 2010), whereby different ES values are assigned to individual classes (Konarska et al. 2002, Hu et al. 2008, Lautenbach et al. 2011).Besides major limitations of land cover based assessment of ES, which were discussed earlier, Verburg et al. (2009) suggested that ES should be better linked to actual land use, which often differs from the actual land cover.Furthermore, the spatial resolution underlying land cover maps influences both the extent and the valuation of ES (Konarska et al. 2002, Di Sabatino et al. 2013).
The second way of using RS to support ES assessment is that RS data are directly used to estimate EP (e.g., biomass), which in turn serve as proxies for the actual ES (e.g., forage production) (Malmstrom et al. 2009).Those vegetation properties that play a key role in radiation absorption and scattering can be well retrieved from the optical RS data.These are for example chlorophyll content (Dash and Curran 2004, Zarco-Tejada et al. 2004, Schlerf et al. 2010), water content (Serrano et al. 2000, Colombo et al. 2008, Clevers et al. 2010) and structural properties such as leaf area index (Fernandes et al. 2004) or fractional vegetation cover (Asner and Heidebrecht 2002).
The most common RS method to study vegetation properties is to use field data to interpret RS images, empirical RS methods that build a statistical regression model between a limited number of field observations and a subset of spectral bands or their combinations (Verstraete and Pinty 1996, Serrano et al. 2002, Huber et al. 2008, Psomas et al. 2011, Tian et al. 2011).By nature, they are computationally fast, easy to implement and therefore suitable for feasibility case studies of smaller spatial extent.The major limitation is that established relationships are less transferable to the same area at a different time, to another location or to another sensor (Grossman et al. 1996, Asner et al. 2003, Colombo et al. 2003).
The objective of this research was to test the potential of high resolution airborne RS data for mapping of EP and ES in subalpine grasslands.Using a common assessment of the supply of services by subalpine grasslands based on their EP, we estimated EP directly from RS data using empirical retrieval methods.The RS-based maps of EP and ES were compared with estimates from the plant trait-based modelling approach combining mapped land use, abiotic variables and plant trait field data (Lavorel et al. 2011).Furthermore, we discussed advantages and disadvantages of both approaches.

MATERIAL AND METHODS
For clarity reasons, we show the main methodological steps in Fig. 1.First, we introduce a study area, input field measurements and remote sensing data.Then we describe the RS approach (this study) and the plant trait-based modelling approach (adopted from Lavorel et al. 2011) to map multiple ecosystem properties and services provided by subalpine grasslands.The mapped v www.esajournals.orgecosystem properties are green biomass, litter mass, crude protein content, species diversity and soil carbon content.The evaluated services are agronomic, cultural and carbon sequestration services, as well as the total delivery of ES.

Study site and field measurements
The study area is located in the Central French Alps, on the south-facing slope of the Villar d'Are `ne municipality (45802 0 22 00 N, 06822 0 07 00 E).It spans an altitude gradient ranging from 1800-2100 m a.s.l.Floristic composition of the subalpine grasslands has been affected by a long history of agricultural and pastoral land use.A detailed description of the study area is given in Que ´tier et al. (2007).Ecosystem properties (green biomass, litter mass, crude protein content, plant species diversity and soil carbon content) and plant traits (vegetation height, leaf dry matter content, nitrogen content and phosphorus content) were measured in 30 3 30 m permanent plots stratified by land use and altitude during years 2004-2009.Detailed protocols for field data collection are provided in Lavorel et al. (2011).Twenty-five field plots out of 57, which were distributed across the entire study area, were covered by the RS images (Fig. 2a).Therefore it is important to mention, that the plant trait-based modelling could benefit from a larger pool of field data than the RS approach.

Remote sensing data
Airborne imaging spectroscopy (or hyperspectral RS) data were acquired using the AISA Dual system (Specim, Finland) on 23 July 2008.We acquired six flight lines covering an area of about 4 3 2.5 km.The resulting ground pixel size was 0.8 m.In total 359 spectral bands covered the range of 400-2450 nm and an average spectral sampling interval between bands varied from 4.3 to 6.3 nm.We discarded about one third of the available spectral bands due to poor radiometric response (i.e., 401-444, 876-1140 and 2452 nm) and atmospheric disturbances (i.e., 752-771, 1334-1485, 1786-1968 nm).In total 240 bands were retained for further analysis.The AISA images were pre-processed in three steps.The AISA flight lines were first calibrated to at-sensor radiance values using sensor-specific radiometric calibration coefficients using the AISA processing toolbox CaliGeo v.4.6.4 (Specim 2009).Second, they were geometrically corrected using onboard navigation data and a local digital elevation model (DEM) of 10-m spatial resolution and subsequently orthorectified to a Universal Transverse Mercator (UTM, Zone 32N, WGS84 datum) map projection.Third, at-sensor radiance data were converted to hemisphericaldirectional reflectance factors (see Schaepman-Strub et al. [2006] for terminology) by applying an atmospheric correction using the ATCOR-4 software (Richter and Schlapfer 2002).A true color composite of the fully corrected AISA Dual data is shown in Fig. 2a.
We applied maximum likelihood classification (ENVI 4.8, ITT Visual Information Solutions) to distinguish grasslands from other land cover classes (trees and bushes, roads and urban areas, and bare soils and rocks).Subsequently, we applied a fully constrained spectral unmixing procedure (Zurita-Milla et al. 2007) on the classified grassland pixels to decompose the mixed spectral signal into relative fractions of three components: fractions of photosynthetic vegetation, non-photosynthetic vegetation, and soil.The actual land cover classification combined with the results of the spectral unmixing is shown in Fig. 2b.This allowed us to mask out all pixels with a green vegetation fraction lower than 25% (according to Okin et al. [2001]) and pixels classified as trees/bushes, roads/urban areas, soils/rocks and heavily grazed or mown grasslands were masked out.In total, about 25% of all AISA pixels were excluded from further analysis.Moreover, field plots located on the masked image area were also excluded and only 18 plots remained for RS data interpretation.

Estimation of ecosystem properties
Five EP, namely green biomass, litter mass, crude protein content, species diversity and soil carbon, were addressed in this study.The plant trait-based modelling approach (further denoted as the modelling approach) was adopted from Lavorel et al. (2011) and here we briefly introduce its main methodical steps (Fig. 1).Spatial distribution of EP was modelled using generalized linear models combining the maps of land use and field plot measurements of abiotic variables (topography and soil properties) and plant traits (vegetation height, leaf dry matter content, nitrogen and phosphorus content).In these EP models, spatial variations of abiotic and plant trait inputs were estimated based on land use and possible modifications reflecting effects of altitude or soil heterogeneity.The modelled EP were normalized between 0 and 100 using 5th and 95th data percentiles as boundary values and the resulting EP maps at a spatial resolution of 20 m were clipped to match the spatial extent of the AISA images (ArcGIS 10.0, ESRI).
We evaluated all possible two-pair band combinations in the range from 440 up to 2450 nm (R i , R j are reflectance values at wavelengths i and j ) to build two types of narrow-band vegetation indices: SML regression was applied to the full reflectance spectra (240 bands) and we evaluated all possible band combinations up to a maximum of four spectral bands in order to avoid model overfitting (Psomas et al. 2011).Regression models with high multicollinearity among selected spectral bands (variance inflation factor VIF .10) were eliminated.SML regression was implemented in R (version 2.15.0;R Development Core Team 2012) using the packages leaps and regr0.
PLS regression, a multivariate technique that reduces the large number of collinear spectral variables to a few non-correlated latent variables, was applied to full, non-transformed reflectance spectra using Matlab 7.11.0,function plsregress.
The low number of field EP data (n 18) prevented using an independent subset for model calibration and validation.Therefore, all EP data were used to build empirical regression models.The overall capability of each empirical model to explain the variability in measured EP was evaluated by the coefficient of determination (R 2 ).The predictive power of a model has been assessed by estimating the root mean square error of prediction (RMSEP) using the leave-oneout cross-validation approach.
The best regression model estimating EP from RS data was selected according to the highest R 2 and lowest RMSEP and it was then applied to the entire RS image at the original spatial resolution of 0.8 m.Resulting maps of EP from RS images were normalized between 0 and 100 using 5th and 95th data percentiles as boundary values and only then resampled to 20-m pixel size using the nearest neighbor resampling method to match the spatial resolution of the modelling approach.

Estimation of ecosystem services
For both, the remote sensing and the modelling approaches, agronomic, cultural and carbon sequestration ecosystem services were calculated in the same manner.Agronomic and cultural services were recognized as the two most relevant services by local stakeholders and carbon sequestration (a proxy of climate regulation service) is a service of global relevance (Lamarque et al. 2011).These three ES were calculated as linear combinations of EP according to coefficients adopted from Lavorel et al. (2011) (Table 1).These coefficients were estimated based on social evaluation on how local stakeholders perceive services provided by subalpine grasslands (Que ´tier et al. 2010, Lamarque et al. 2011).The agronomic service was calculated as the sum of fodder quantity (green biomass) and fodder quality (crude proteins).The cultural ES was calculated as the sum of the positive effect of species diversity, which increases aesthetic and cultural values of the area, and the negative effect of litter mass, i.e., accumulation of litter mass, which on the contrary, according to stakeholders, decreases the aesthetic and cultural value by demonstrating poor management.The carbon sequestration service was approximated by soil carbon content.The total ecosystem service supply was calculated as the sum of individual The location of hot spots and cold spots of total ecosystem service supply was analyzed by applying the Getis-Ord-Gi* spatial statistic in ArcGIS 10.0 (ESRI).The algorithm identifies statistically significant spatial clusters of high (hot spots) and low (cold spots) values by testing a null hypothesis whether the observed clustering is more pronounced than one would expect in a random distribution of those same values.

Spatial comparison
Spatial assessment of differences between the RS and the modelling approaches was implemented using Map Comparison Kit 3.2 (Visser and de Nijs 2006), a toolbox allowing spatial comparison of categorical maps using the fuzzy set theory (Hagen 2003).The fuzzy approach has two great advantages as compared to a direct pixel-by-pixel comparison.It allows assessing the spatial neighborhood of a pixel (i.e., closer pixels are more similar than distant pixels) and proximity among classes (i.e., some classes are more similar to each other than other classes) (Hagen 2003).For the purpose of spatial comparison, all maps were classified by grouping values to equally spaced bins.The spatial comparison was repeated at two spatial resolutions: 20 m and 100 m pixel size in order to identify how the scale of heterogeneity affects the comparability of the RS and the modelling approaches.
Similarity between the RS and the modelling approaches was expressed spatially (similarity maps) and by calculating the percentage of average similarity (a fuzzy equivalent of overall classification accuracy obtained from a confusion matrix).Moreover, both approaches were compared by using descriptive statistics and frequency distributions calculated for the entire image area, as well as for the individual land use classes.We used the Wilcoxon-Mann-Whitney U test ( p 0.05) to identify statistically significant differences between the approaches.
The land use classification (shown in Fig. 2c), which served to compare the RS and the modelling approaches, was fully adopted from Lavorel et al. (2011).Seven types of land use (LU) were found within the area covered by the AISA images.Three LU classes were located on terraces that have been previously cultivated: fertilized and mown terraces (LU1), non-fertilized but mown terraces (LU2), and unmown terraces but grazed in spring and autumn (LU3).Two LU classes were located on never cultivated permanent meadows: mown grasslands (LU4), and unmown grasslands but grazed in summer (LU5).Finally, one class covered alpine-grazed meadows at altitudes above 2000 m a.s.l (LU6) and the last class represented grasslands growing on steep and grazed slopes (LU7).

Ecosystem properties: comparison of remote sensing and modelling approaches
The best empirical models for predicting ecosystem properties from RS data are summarized in Table 2.The best model for green biomass estimation was based on a narrow-band ND index, whereas SML provided the most accurate solution for other EP.Interestingly, PLS v www.esajournals.orgdid not outperform SML despite the fact that some studies found PLS more powerful than SML or vegetation indices (Cho et al. 2007, Darvishzadeh et al. 2008).The prediction accuracy of the RS-based empirical models was moderate and of a similar order as those obtained from the modelling approach.Green biomass was estimated with the lowest accuracy among all EP (R 2 ¼ 0.54), which was lower than the modelling approach (R 2 ¼ 0.70).The prediction accuracy of the litter mass model (R 2 ¼ 0.60) was comparable to the modelling approach (R 2 ¼ 0.61).The prediction accuracy of the crude protein content model (R 2 ¼ 0.97) was higher than from the modelling approach (R 2 ¼ 0.62), but the reliability of the RS-based model is strongly reduced due to low degrees of freedom, which resulted in an over-adjusted model.The prediction accuracies of the species diversity model (R 2 ¼ 0.60) as well as of soil carbon content (R 2 ¼ 0.73) were higher than for the modelling approach (both achieving R 2 ¼ 0.31).Additionally, Fig. 3 shows scatter plots between measured and RS-predicted EP, indicating that estimated EP from RS had a comparable range as field measured EP.
Resulting maps of EP from the RS (left panels) and the modelling (middle panels) approaches at 20-m spatial resolution are presented in Fig. 4. The visual comparison between approaches emphasizes spatial consistency and continuity of the RS approach, contrasting with the land-use based pattern underpinning the modelling approach.The RS approach was thus highly sensitive to fine scale heterogeneity, which could be captured thanks to the high spatial resolution of the RS images.Therefore, the variability of RSbased EP within individual LU classes was up to ten times higher than the variability of the modelling approach (Appendix A).In contrast to the modelling approach, none the RS-based maps of EP clearly captured a distinct difference between grasslands growing on previously cultivated terraces (LU1-3 located in the southwest part of the image area) and grasslands at permanent meadows (LU4-7 located in the northeast part of the image area).Statistical analyses revealed that for all LU classes the mean EP value from the RS approach was always significantly different from the modelling approach (Wilcoxon-Mann-Whitney U test, p 0.05).For instance the modelling approach substantially provided greater values (resp.lower values) of green biomass for LU1 and LU5 (resp.LU3, LU6, LU7), and greater values of species diversity for all land use types except LU3 and LU5 as compared to the RS approach (Appendix A).
The spatial comparison between RS and modelling approaches based on the fuzzy set theory (Fig. 4) indicated that the average similarity was around 20% at the spatial resolu-Fig.4. Ecosystem properties estimated from the remote sensing approach and the modelling approach at a spatial resolution of 20 m.The most right maps show similarity between remote sensing and modelling approach.Frequency histograms show the distribution of values within the image.v www.esajournals.orgtion of 20 m.Crude protein content, for which there was the lowest number of field samples to calibrate RS images, was the EP with the lowest spatial agreement between the approaches (average fuzzy similarity calculated for entire image area was only 16.6%).Species diversity was the EP with the highest spatial agreement between the approaches (average fuzzy similarity reached up to 22.4%).In classes LU2, LU3 and LU7 the average similarity reached about 25%, whereas in LU1 (mown-fertilized terraces) it was only about 8%.We also investigated similarity of RS and modelling approaches at the spatial resolution of 100 m and the average fuzzy similarity increased from 20% up to 40% at the 100 m resolution (Appendix B: Fig. B1).

Ecosystem services: comparison of remote sensing and modelling approaches
Ecosystem services were calculated as simple linear combinations of ecosystem properties according to the weights presented in Table 1.Resulting maps of ES as derived from the RS (left panels) and the modelling (middle panels) approaches are presented in Fig. 5 (note that maps of carbon sequestration service are presented in Fig. 4m, n, o as soil carbon content).
The statistical analysis carried out at LU class level indicated that for almost all classes the mean ES value from the RS approach was significantly different from the modelling approach (Wilcoxon-Mann-Whitney U test, p 0.05).In particular, the agronomic value from the modelling approach was greater (resp.lower) for LU1, LU2 and LU6 (resp.LU3 and LU7) as compared to the RS approach.The cultural value from the modelling approach was greater for all land use types except LU3 and especially LU5 (Appendix A) than from the RS approach.
The spatial comparison between approaches based on the fuzzy set theory (right panels of Fig. 5 and Fig. 4o) indicated that the average Fig. 5. Ecosystem services estimated from the remote sensing approach and the modelling approach at a spatial resolution of 20 m.The most right maps show similarity between remote sensing and modelling approach.Frequency histograms show distribution of values within the image.
v www.esajournals.orgsimilarity for the agronomic and the carbon sequestration services was around 20%, but for the cultural ES and the total ES it was less than 10%.Analogous to EP, LU classes with highest spatial agreement across ES were LU3 and LU7, where the average similarity reached about 20%.In LU1 (mown terraces) the spatial agreement was the lowest (about 2%).At the spatial resolution of 100 m (right panels of Figs.B1o and B2 in Appendix B), the average similarity increased to a maximum of 50% for the agronomic service, to 40% for the carbon sequestration service, and to 30% for cultural and total ES compared to 20-m resolution.
Descriptive statistics for the total ecosystem service supply indicated that the modelling approach over-estimated (resp.under-estimated) total ES for LU1, LU2, LU 4 and LU6 (resp.LU3 and LU5) (Fig. 6).Further spatial analysis of hot and cold spots of total ES supply (Fig. 7) indicated generally low agreement between the approaches.The area of hot/cold spots from the modelling approach was about double compared to the RS approach.Moreover, the RS-based hot/ cold spots were small scattered patches across the study area, rather than compact areas associated with specific land use types.This contradiction reflected the decrease in contrast among land uses in the RS as compared to the modelling approach (Appendix A).
The modelling approach clearly identified alpine grazed meadows (LU6) as ES hot spots.These grasslands provided a synergy between high species diversity, soil carbon and crude protein content and low litter mass and the same synergy pattern was underlying the RS hot spots.The hot spots predicted by the RS approach appeared as a subset of those predicted by the modelling approach (40% of hot spots predicted by the RS approach).
The modelling approach clearly identified unmown-grazed permanent meadows (LU5) as cold spots.These grasslands exhibited high values for litter mass and green biomass and lower values for the three other EP.Similarly, many RS-based cold spots could be associated with the combination of high litter mass and low values for all other EP.The RS-based cold spots exhibited a small spatial overlap with the modelling approach (17% of cold spots predicted by the RS approach).They corresponded well to pixels with lower vegetation fraction that are often scattered around roads and terrace banks.

Quality of empirical RS models of ecosystem properties
In general, three major sources of uncertainties linked to the empirical RS models used in this study can be identified.They can be related to the quality of input RS data, to statistical methods and selected spectral bands, and to the quality and representativeness of field measurements used to interpret RS images.We minimized the uncertainties for RS data by removing noisy bands and pixels with low vegetation signal, but uncertainties depending on the quality of atmospheric and topographic corrections remained.Recent work of Laurent et al. (2011b) indicates that one way to effectively suppress the uncertainties from the atmospheric correction is to retrieve vegetation properties from the at-sensor radiance data.As it is beyond the scope of this work to investigate the atmospheric and topographic effects on the EP estimates, we focus our discussion on the other two types of uncertainties.
Empirical RS models predicting EP were solely based on the statistical selection of spectral bands.This approach is sometimes criticized because the selected bands are often not related to wavelengths with known absorption of leaf biochemical constituents (Grossman et al. 1996).The green biomass model was constructed from the ND index using one spectral band in the red chlorophyll absorption (682 nm) and the other one in the short-wave infra-red (SWIR) region (2441 nm).Psomas et al. (2011) also found that ND indices utilizing the SWIR bands were more accurate predictors of grassland biomass than ND indices based on red and near infra-red (NIR) wavelengths.In the litter mass model, a band in the NIR region (871 nm), which is mainly influenced by the canopy structure, and two bands in the SWIR region (1158 and 2308 nm) were combined.The first SWIR band is influenced by the water content and the other one can be associated with nitrogen and protein absorption (Kumar et al. 2001).In the model predicting crude protein content, two bands in the SWIR region (2038 and 2346 nm) were selected.Both Fig. 6.Comparison of the remote sensing approach and the plant trait-based modelling approach to estimate total ecosystem service supply per land use class (LU1 fertilized and mown terraces, LU2 non-fertilized but mown terraces, LU3 unmown terraces grazed in spring and autumn, LU4 mown grasslands, LU5 unmown grasslands grazed in summer, LU6 alpine-grazed meadows above 2000 m a.s.l., LU7 grasslands on steep and grazed slopes).The left panels show frequency distribution of the total ecosystem service for the remote sensing approach (black bars) and for the modelling approach (grey bars).The right panels show boxplots, where central line in a box is median, box height is interquartile range representing 50% of the data, whiskers are minimum and maximum unless the observed values exceeded 1.5 times the interquartile range and in that case they are marked with crosses as outliers.The star symbol indicates that differences in median values between two approaches are significant ( p 0.05, Wilcoxon-Mann-Whitney U test).
bands were located near a known absorption feature of protein and nitrogen (Kumar et al. 2001).Nevertheless, the reliability of this model is strongly reduced due to low degrees of freedom.The species diversity model combined a spectral band from the red chlorophyll absorption (706 nm) and from the SWIR region (2447 nm).There are only a few studies relating species diversity to spectral reflectance.For example, Carter et al. (2005) found that the best wavelengths for predicting species diversity prairie grass from AVIRIS data were located in the NIR region.In the soil carbon model the selected bands cannot be clearly associated with the soil property because of the vegetation cover.One spectral band was located in the red-edge region (734 nm), which can be linked to plant stress (Kumar et al. 2001), and two in the SWIR region (1328 nm associated with water absorption [Kumar et al. 2001] and 1686 associated with lignin absorption [Thenkabail et al. 2004]).
The last source of uncertainties in the RS models of EP is related to the quality and representativeness of field EP data.Broad variability in EP was captured by stratifying field sampling plots according to classes of past and present land use (Que ´tier et al. 2007) and altitude (Lavorel et al. 2011).The field sampling design, however, did not capture finer variations within the LU classes resulting from soil and microclimatic variability, specific LU history, and especially time since last mowing.Indeed, complementary sampling revealed that in particular within LU classes with grazing on former hay meadows (LU3, LU5), there was a significant effect of the date of mowing cessation, which influenced plant functional composition and associated EP such as green biomass, litter mass and species diversity (Caubet 2009).The variability in some EP was further reduced because of removing some plots that were located at areas with low green vegetation signal as a result of mowing or grazing at the time of the RS campaign.This negatively influenced mainly the litter mass model, because the plots that were removed had also higher values of litter mass (especially in LU5).Furthermore, it is also important that field sampling design reflects a spectral variability within RS images as well.Analysis of the spectral representativeness of field plots indicated that about 25% of the pixels were outside the spectral range covered by the calibration dataset.
Probably the strongest factor negatively influencing the quality of the RS empirical models was an asynchronous acquisition of field EP data and RS images.Field EP data were collected between years 2004 and 2009, which smoothened variability among the vegetation seasons.In contrast, RS images represent one snapshot in time capturing the actual ecophysiological and development stage of grasslands during the second half of July 2008, which cannot exactly be compared with inter-annual mean field EP data (although 2008 appeared to be representative of an average year [Lavorel and Grigulis, unpublished data]).Further, the date of RS acquisition was not optimal as compared to v www.esajournals.orgvegetation phenology and agricultural use.A good example that illustrates the temporal disparity is green biomass estimation in unmown-grazed summer meadows (LU5).Average green biomass value in this LU class was the highest for the modelling approach, but the lowest for the RS approach.This can be attributed to a large proportion of grassland pixels that were grazed at the date of RS acquisition, whereas the modelling approach was based on peak biomass measurements prior to mowing or grazing (first half of July).Likewise, field crude protein content was estimated at peak vegetative growth in the middle of June, where differentiation among vegetation types is known to be greatest, whereas the later date for RS measurement resulted both in an average lower value (more advanced phenology) and a convergence among vegetation types (Duru 1997).We can assume that both green biomass and crude protein content are the most dynamic EP within and between years and can be highly sensitive to disturbances such as drought in 2004 or vole outbreaks in 2009-2010.

Alternative RS solutions for ecosystem properties estimation
There are clear limitations and uncertainties behind the RS-based models of EP as discussed in the previous section, but there is large variety of other types of RS data and methods, which might provide more suitable solutions to predict currently investigated or other EP, as well as some plant traits (Homolova ´et al. 2013).
Green biomass and litter mass can be possibly related to the results of spectral unmixing, which allows subdividing the mixed nature of the spectral information into relative fractions of soil, photosynthetic vegetation (PV) and non-photosynthetic vegetation (NPV) (Roberts et al. 1993, Asner et al. 2005).We found a positive correlation between field data on green biomass and the PV fraction (r ¼ 0.57, results not shown), and between litter mass and the NPV fraction (r ¼ 0.81, results not shown).This finding and similar patterns found by Numata et al. (2007) encourage a hypothesis that PV and NPV fractions can serve as proxies for green and dry vegetation biomass, respectively.
Crude protein content is essentially proportional to nitrogen content and a moderately strong relationship exists between nitrogen and chlorophyll content (Homolova ´et al. 2013).Chlorophyll content is one of the vegetation properties that nowadays can be estimated from RS with high fidelity using a variety of chlorophyll-sensitive vegetation indices (Sims andGamon 2003, Gitelson et al. 2006) and physical retrieval methods based on radiative transfer modelling (Zhang et al. 2008, Jacquemoud et al. 2009, Malenovský et al. 2013).Therefore it can be used as an operational proxy for crude protein content estimation.
In the case of species diversity, some authors suggested to approximate it by indicators of spectral heterogeneity derived from RS images (Rocchini 2007, Rocchini et al. 2010).This type of analysis usually requires RS data of lower spatial resolution (tens of meters) and more spectrally contrasting land cover types (i.e., grasslands vs. forests).Therefore, we expect that the applicability of this alternative approach for our data will be limited, although it can be a suitable on large scale mapping of species diversity.
Mapping of soil properties is generally difficult in areas covered by vegetation (Bartholomeus et al. 2011).Some studies suggest that soil carbon content correlates with other vegetation properties such as leaf area index or canopy height (Li et al. 2010), but screening of our field data did not indicate any significant relationship with measured plant traits or EP (results not shown).Here, soil carbon content was used as a proxy for the carbon sequestration service, given that in mountain grasslands aboveground biomass senesces and largely decomposes every year.RS based solutions that directly estimate carbon stocks sequestered by vegetation, i.e., primary productivity (Running et al. 2004, Zhao et al. 2005), may be considered for other types of ecosystems such as forests.Preferably still, RS data can also serve as an input into ecosystem models predicting carbon stocks (Porfirio et al. 2010).

Comparison of remote sensing and modelling approaches
The visual comparison of EP as well as ES at the spatial resolution of 20 m (cf.left and right panels in Fig. 4 for EP and in Fig. 5 for ES) contrasted the spatial consistency and continuity of the RS approach with the land-used based v www.esajournals.orgpattern of the modelling approach.Advantages and disadvantages of both approaches to quantify EP and consequently ES are closely related to attributes of field data, as well as to underlying spatial data and their spatial and temporal resolutions.The advantages/disadvantages are summarized in Table 3 and further discussed in following sections.
The spatial resolution of underlying spatial data plays an important role in assessment of ES.A coarser spatial resolution of the modelling approach (20 m) did not fully reflect spatial heterogeneity within the study area as it was discussed in section Quality of empirical RS models of ecosystem properties.Moreover, it might have easily included non-vegetated areas assigning them an unrealistic value or left out small vegetation patches with high ES supply.Di Sabatino et al. (2013) demonstrated that land cover data at a spatial resolution of 20-30 m can significantly underestimate total area of highly fragmented ecosystems such as water bodies.The high spatial resolution of RS data helped to eliminate non-vegetated areas, which usually have low EP and thereby ES values.EP models based on RS data are much more sensitive to spectral variations caused by micro-topography and soil background.Small-scale case studies like this one can certainly benefit from the high spatial resolution, airborne RS data, but they are less suitable for larger scale studies.Therefore satellite RS instruments such as Landsat or Sentinel-2, which provide a good compromise between spectral (;10 bands) and spatial (10-60 m) resolutions, are more appropriate for larger scale assessment of ES.
Both approaches addressed a temporal aspect of the analysis differently.As discussed in section Quality of empirical RS models of ecosystem properties, the modelling approach is less dependent on inter-annual variability, because field data on plant traits and EP were averaged across several vegetation seasons.Moreover, the underlying LU classification is almost time independent as it reflects current and past land use management.In contrast, the RS approach reflects an actual development stage and management activities and therefore it is less generic than the modelling The dimensionality of the underlying spatial data to model EP differed between the approaches.The modelling approach builds only on two spatially explicit, independent landscape input variables (LU and DEM maps), whereas RS images offered initially 240 spectral bands for EP estimation.However, spectral bands can be highly correlated and therefore the spectral dimensionality can be significantly reduced, which allows independent estimation of a maximum of 3-12 variables from a single-view image (Verhoef 2007, Laurent et al. 2011a).
The RS approach allowed skipping one important step of the modelling approach, in which plant traits and abiotic properties were estimated from LU and DEM.By skipping this step, the number of field data on plant traits and abiotic properties can be significantly reduced.Nevertheless, field data collection of EP will still remain a limiting factor for large scale and long term applications.
A big disadvantage of the RS empirical methods is that transferability of regression equations to other test sites or other type of RS data is limited (Grossman et al. 1996, Asner et al. 2003, Colombo et al. 2003).Similar constraints apply to some extent to the modelling approach, although trait based models can provide robust inter-site predictions of EP (Fortunel et al. 2009).Luckily, RS offers more universal solutions based on physical radiative transfer models, though their retrieval capacity is limited to a few vegetation properties, such as leaf area index, chlorophyll content and water content (Zarco-Tejada et al. 2004, Schlerf et al. 2005, Baret and Buis 2008, Colombo et al. 2008, Jacquemoud et al. 2009).
Finally, the scheme combining the underlying EP of subalpine grasslands into ES was the same for both approaches.Although it can raise concerns, it was based on a social evaluation as perceived by local stakeholders (Que ´tier et al. 2010, Lamarque et al. 2011).This makes it very specific for the actual study area, but we can assume that the scheme can be transferred to other subalpine regions with similar landscape management and ES demand (Lamarque et al. 2011).However, its applicability to more diverse landscapes needs to be further tested.

CONCLUSIONS
In this case study, we tested high spatial and spectral resolution RS data for the spatial assessment of ecosystem properties (EP) and multiple ecosystem services (ES) in subalpine grasslands.Furthermore, the RS-derived maps of EP and ES were compared with results of the plant trait-based modelling approach proposed by Lavorel et al. (2011), which predicted EP using general linear models combining maps of land use, abiotic variables and plant traits.The prediction accuracy of the RS-based empirical EP models was good (R 2 ¼ 0.54-0.73)and it was similar or better than those obtained from the modelling approach.
The average similarity between the RS and the modelling approaches at the 20 m spatial resolution varied between 16 and 22% for individual EP, but dropped to about 7% for the total ecosystem service supply.The spatial agreement increased by about 20 percentage points when increasing pixel size from 20 to 100 m, showing that some of the discrepancies among approaches resulted from their differences in intrinsic spatial resolution.Despite the fact that we found similarities in the bundles of EP underpinning the total ES in both approaches, the exact location differed as the RS-based hot/ cold spots were smaller scattered patches within the study area, owing to decreased variability in EP within land use types in the RS as compared to the modelling approach.
RS-based maps provided more detailed insight into the fine-scale spatial distribution of EP and thereby ES, because RS data are much more sensitive to small scale heterogeneity than the LU map.RS-based results of EP exhibited about ten times higher variability within individual LU classes than the modelling approach, which reflected by construction the underlying LU classification.
We are aware of the limitations of the RS approach presented in this study, which mainly originated from the asynchronous acquisition of RS and field data, both in terms of specific year, and especially of timing within the vegetation season.But without any objective reference, it is hard to decide which of the two approaches is better, as both are imperfect assessments of a reality and both having their advantages and disadvantages.The underlying spatial data, RS images for the RS approach and LU map for the modelling approach, are the major source of the discrepancies between the two approaches.Therefore an approach that combines RS images, LU and DEM maps with field EP observations would probably be the best compromise for accurate spatial assessment of multiple ES supplied by subalpine grasslands.
In general, future research should look for reliable links between ES and RS proxies.This will require more case studies, where RS and ES data (and EP determining ES) are acquired and analyzed together.The assessment of ES can already benefit from those vegetation properties, for which RS methods are well established, e.g., chlorophyll content, leaf area index, fraction of absorbed photosynthetically active radiation, fractional vegetation cover.Besides the spatial and spectral domains of RS data that were explored in this study, future assessment of ES should also consider multi-temporal RS data.Multi-temporal data can help to address the phenology of subalpine grasslands, which is a relevant property for the assessment of the agronomic service.

Fig. 2 .
Fig. 2. (a) True color display (R ¼ 640 nm, G ¼ 552 nm, B ¼ 462 nm) of the airborne imaging spectroscopy data acquired in the spatial resolution of 0.8 m by the AISA Dual sensor on July 23, 2008.The AISA image captures the study area in the Central French Alps (Villar d'Are `ne, 45802 0 22 00 N, 06822 0 07 00 E) and yellow dots indicate the location of field plots.(b) Current land cover classification derived from supervised image classification and spectral unmixing applied on AISA images (pixel size 0.8 m).(c) Land use classification as adopted from Lavorel et al. (2011) is shown at the spatial resolution of 20 m.Unclassified areas are in black.

Fig. 3 .
Fig. 3. Comparison of field measured and RS-predicted ecosystem properties: (a) green biomass (g/m 2 ), (b) litter mass (g/m 2 ), (c) crude protein content (g/kg), (d) species diversity expressed as Simpson's diversity index (unitless), and (e) soil carbon content (%).Predicted EP values represent here a mean value per plot calculated from pixels within a plot, which have green vegetation fraction .25%.The solid line represents the one-to-one line.

Fig. 7 .
Fig. 7. Location of hot spots (a) and cold spots (b) of the total ecosystem service supply in subalpine grasslands as derived by the remote sensing and the plant trait-based modelling approaches at a spatial resolution of 20 m obtained by using the Getis-Ord-Gi* spatial statistics.Grey areas are grasslands with an average supply of total ecosystem services.

Fig. A1 .
Fig. A1.Descriptive statistics as derived from the spatial maps of ecosystem properties and services at 20 m spatial resolutions show the comparison of the remote sensing approach and the plant trait-based modelling approach per land use class.The evaluated ecosystem properties are (a) green biomass, (b) litter mass, (c) crude protein content, (d) species diversity and (e) soil carbon content, and the ecosystem are (f ) agronomic, (g) cultural, and (h) total supply.Symbols are as in Fig. 6.

Table 2 .
Summary of the best regression models estimating ecosystem properties from the AISA remote sensing data.For comparison purposes, the last column shows R 2 as obtained from the plant trait-based modelling approach.ND ¼ narrow-band normalized difference ratio, SML ¼ stepwise multiple linear regression, R i ¼ reflectance at wavelength i, df ¼ degrees of freedom, R 2 ¼ coefficient of determination, RMSEP ¼ root mean square error of prediction relative to mean observed value.

Table 3 .
Overview of advantages (þ) and disadvantages (À) of the remote sensing approach and the plant traitbased modelling approach for assessment of ecosystem properties and services in subalpine grasslands.