The Integration of Ecosystem Services in Planning: An Evaluation of the Nutrient Retention Model Using InVEST Software

Mapping ecosystem services (ES) increases the awareness of natural capital value, leading to building sustainability into decision-making processes. Recently, many techniques to assess the value of ES delivered by different scenarios of land use/land cover (LULC) are available, thus becoming important practices in mapping to support the land use planning process. The spatial analysis of the biophysical ES distribution allows a better comprehension of the environmental and social implications of planning, especially when ES concerns the management of risk (e.g., erosion, pollution). This paper investigates the nutrient retention model of InVEST software through its spatial distribution and its quantitative value. The model was analyzed by testing its response to changes in input parameters: (1) the digital terrain elevation model (DEM); and (2) different LULC attribute configurations. The paper increases the level of attention to specific ES models that use water runoff as a proxy of nutrient delivery. It shows that the spatial distribution of biophysical values is highly influenced by many factors, among which the characteristics of the DEM and its interaction with LULC are included. The results seem to confirm that the biophysical value of ES is still affected by a high degree of uncertainty and encourage an expert field campaign as the only solution to use ES mapping for a regulative land use framework.


Introduction
The interest in understanding the implication of land use changes on different ecosystem services (ES) is garnering attention, using both biophysical and economic values as a proxies of soil degradation [1][2][3][4]. Among others, methodologies for ES accounting at different scales have been developed, and new tools for mapping ES are now available (e.g., Integrated Valuation of Ecosystem Services and Tradeoffs-InVEST, Artificial Intelligence for Ecosystem Services-AIRES, Land Utilization and Capability Indicator-LUCI) [5][6][7][8]. The assessment of ES through their spatial distribution is a key element to set policy targets for sustainability and resiliency because their visualization highlights the trade-off between different alternatives [9][10][11]. Nonetheless, weakness arises when ES mapping is provided at the local scale for regulative planning purposes when the assessment should support the decision-making process regarding land use properties and building rights. In that phase, it seems that the experience of mapping is still in the early stage of development [3,7,12] and the debate on how different methods of mapping ES influence the stakeholder' s interaction aimed at taking sharp decisions on land use configurations [2,13,14] has not been adequately opened.
Such a gap between general ES analysis and its effective implementation at the local stage is underlined by many publications [2,7,13]. The gap is mainly due to a different kind of issue: (i) the ES principles, concepts, and mapping are not yet currently used to define the land use regulation at regional, sub-regional, and local stages [9,[21][22][23]; (ii) the technical knowledge of the tools for mapping ES and their utilization during Strategic Environmental Assessments for plans and programs is limited [22,24,25] and; finally, (iii) even when ES are accounted for and considered at the local stage, it is difficult to define how each ES should contribute to achieving an efficient land use plan in terms of the whole, and integrated ecosystem gains rather than of a single function [16,26,27]. Above all considerations, different policy applications require different accuracies, scales, and spatial modeling approaches to set targets and objectives [1,17,22]. Thus, the in-depth comprehension of input-output relations of each ES model became a key point to integrate ES in planning.
The aim of this study is to demonstrate how specific ES values are reliable enough to be used in the decision-making process for land use planning at the local stage. An in-depth assessment of the nutrient retention model defines how the output map is influenced by inputs and, in particular, how the digital terrain elevation model (DEM) affects the delivery of nutrients along the hill slopes in a selected study area. Inasmuch as planners are not experts in the field of ES mapping, the knowledge of model sensitivity to input parameters is essential to create awareness on how critical the use of the output should be when maps are used to set territorial policies rather than land use regulations.
The influence of DEMs with respect to the final output was tested looking at both the spatial distribution of the biophysical value and the quantity of nutrients delivered to the landscape; similarly, the influence of land use/land cover (LULC) data has been tested observing how shifting the load/retain parameters linked to LULC classes affect the output values. After a general analysis, a test of sensitivity has been applied to DEM investigating how different resolutions (10,20,25, and 75 m) affect the total export of nutrients in the streams. The behavior of six topographic DEM-related indices have been analyzed selecting a portion of the catchment area characterized by high slopes with a pre-alpine environment and with mountains and hills (the northwestern part of the study area), thus to reduce the time for model processing and obtain a less time-consuming output. A probability distribution has been assigned to DEM variables and the response of the outputs are analyzed.
The results demonstrated how a non-expert utilization of mapping tools should generate possibly distorted evaluations during planning activities. The awareness of how a mapping tool is influenced by inputs is a key element to control the ES mapping process, especially when outputs are used to set planning prescriptions, such as limitation, mitigation, or compensation measures for land use transformations [15].

The Research Context
The research LIFE SAM4CP held by the DIST (Interuniversity Department of Regional and Urban Studies and Planning, Politecnico di Torino. DIST is a partner of the LIFE research concerning ES mapping activities. The research group includes the scientific coordinator Prof. Carlo Alberto Barbieri and collaborators Prof. Angioletta Voghera, Prof. Giuseppe Cinà, Prof. Carolina Giaimo, and the technical staff, Francesco Fiermonte, Gabriella Negrini Costanzo Mercugliano, and Marcella Guy) provides a detailed assessment of the spatial distribution of the principal ES using InVEST (Integrated Evaluation of Ecosystem Services and Tradeoffs) [28]. The aim of the research is to carry out the spatial distribution of biophysical values in planning processes, using maps to define an ecological and social assessment of the land use transformations at the local scale and, thus, providing better solutions for communities involved in local planning [11].
The research selected seven ES (habitat quality, carbon sequestration, water yield, nutrient retention, sediment retention, crop pollination, and crop production) to evaluate the sustainability of urban transformations [29][30][31]. InVEST software (InVEST is a suite of free, open-source software models used to map and value the goods and services from nature that sustain and fulfill human life. The software is a product of the partnership between Stanford University and the University of Minnesota, The Nature  Conservancy,  and  the  World  Wildlife  Fund.  To  download  the  software,  go  to: http://www.naturalcapitalproject.org/invest/) [28] is one the primary open access products of the Natural Capital Project, and even if it appears as an easy-to-use software (with on-line training guides) its utilization implies expert knowledge and a collection of numerous amounts of input data (both statistical, with .csv files, and associated geometries with relative raster cell values, rather than vector polygons in .shp format).
The LIFE SAM4CP uses InVEST as the most common and diffuse software to map and assess the biophysical values of the different ES [32,33], nonetheless, the outputs of InVEST were subsequently analyzed in a GIS environment to evaluate the spatial distribution and the quantity of biophysical maps [34][35][36].
When the territorial dimension of the project was defined (municipalities of Bruino, Settimo Torinese, Chieri, and None were selected as the testing areas), then the collection of input data for modeling each ES was setup. The mapping context has been extended to a broad squared selection of Piedmont territory (50 km × 50 km) centered in the pivot city of Turin, thus to avoid the so-called "edge effect" [28] in the testing municipalities ( Figure 1). Particularly, the mapping context is a selection of the larger catchment area of the Po River. The mapping context spans 312,525.28 ha and it constitutes 33.8% of the territorial dimension of the metropolitan city of Turin (682,701.11 ha). The selection comprises 26% of plain areas (altitudes below 100 m above sea level (a.s.l.)), 21% of hilly areas (altitudes ranges between 100 and 600 m a.s.l.), and 52% of mountainous area (altitudes above 600 m a.s.l.). Particularly, 17.7% is composed of anthropic areas, 52.7% by agricultural areas, 28.5% by natural and semi-natural zones, and 1.0% by water bodies ( Table 1).
The metropolitan area of Turin is composed by a mixed and heterogeneous landscape: from the suburban hills with a semi-dense built-up zone, to the dense and continuous built-up area of the Turin center. The natural mountain zones, and the ones with specific crop management and terraces, to plain floor valleys with a mixed rural development with intensive agricultural productivity, but subjected to high urban expansion. The settlement system of Turin is axially distributed, with a ribbon development along the axes to the northwest (France), northeast (Milan-Venice), and south (Genoa).
After the economic boom in the flat parts of the infrastructural corridors, rather than in the Turin valleys, the pressures of anthropic activities generated a high concentration of land take with the concentration of built-up areas associated with the presence of road networks at the primary and secondary level [37].

The Planning Outcomes
The nutrient retention model track, with the maximum possible reliability both in terms of geometrical precision and accuracy of the associated biophysical values, is a pixel-based routing of the nutrients that come from diffuse sources in the territory (the composition of the different land uses in the area of investigation) to the delivery areas. When it rains, water flowing over the landscape routes down the slopes carrying pollutants from diffuse sources into streams, rivers, or lakes. Delivery areas are the parts where the water routing model generates a stream, so that the nutrient, which has not been previously retained by the landscape, routes downslope via the surface or subsurface into streams, causing pollution in the streams [28].
Land-use planners need information regarding the contribution of ecosystems to mitigate water pollution. Specifically, they require spatial information on nutrient export and areas where the highest filtration happens. The nutrient retention model provides this information for non-point source pollutants, but it generates a spatially-explicit prefiguration of the areas where pollutants flow into streams; thus, it is possible to reinforce vegetation filtering services to avoid water pollution.
Considering the geometrical resolution of DEMs, many studies found that better model performances are reported in the range of four to 10 m, while coarser pixels generate poor performances [38][39][40]; however, at the same time, it seems that no linear relationship is guaranteed between finer DEM resolution and better reliability, in terms of model output.
As an example, the sediment retention model, which allows calculating the soil erosion, largely depends on the intermediate water routing process that generates the runoff index and, subsequently, additional inputs on land use management interacting with the water flow model. Some cases studies argue that a field campaign after a landslide confirms better feedback with predictive models generated with coarser DEM resolutions.
Analogously, the nutrient retention model consists of a first intermediate process which generates the runoff index which is DEM-dependent, then the results of the intermediate process are used to interact with nutrient loading and absorbing values dependent on LULC. Obviously, if the preliminary interaction with the DEM generates poor performances, then the model increases the errors influencing the final output, which will probably fail to represent a reliable and detailed output usable for decision-making processes.
Mapping nutrient retention is pivotal for the SAM4CP project because the effects of anthropic activities (both for rural and urban uses) on water quality are scarcely used to steer the decisionmaking process for planning purposes. The standard approach to water pollution control is the reduction of anthropogenic sources in the territory, while it is less considered that the soil has a high power to retain and abate the concentration of pollutants before they reach stream areas.
For instance, if the knowledge of load and retention areas are spatially defined by the nutrient path to streams, the new ecological compensations for urban transformations should consider the plantation of new vegetated areas to increase nutrient removal along the nutrient path, achieving a great result in terms of ecological connectivity, but also in terms of filtering pollutants.
The "filtering" service should be empowered by planning decisions because the vegetation can remove, store, or release the pollutants in other forms. A typical example of the service is provided by the creation of riparian zones. Riparian zones can slow the flows, enabling the removal of pollutants taken up by vegetation and serves as a barrier before pollutants reach the stream. If the knowledge of the nutrient export is spatially represented, then mitigation or compensation measures can be designed by the land use plan as specific contribution to improve the environmental condition of the territory.
Moreover, if biophysical values of the pollutant export are accounted for by the current state and the planned ones, then an in-depth observation of the contaminant trend for alternative land use options should be evaluated. Particularly, the conversion of agricultural land uses, the reduction of natural or semi-natural zones associated with the expansion of urban land uses, represent a threat to water quality. If one of the major targets of contemporary planning is to find sustainable paths for planning transformations, then the trade-off between the land use configuration and the particular variation of different ES should be considered.
Considering the above-mentioned purposes, the research required a technical validation of the nutrient retention model as a key action to understand its direct utilization.

The Methodology of Validation
In the mapping stage, the spatial distribution of biophysical values for each ES was used as a proxy of soil quality for planning purposes. Such an assessment was necessary since the project was aimed at defining spatial planning rules to limit, mitigate, or compensate for soil sealing [17,41], acting with prescriptive land use regulation [42]. For example, one of the targets was to decide whether or not confirming a potential land use transformation considered the soil ES: if a potential alteration was affecting soils with a high ES provision, then the option of limiting the transformation should be considered, whereas when transformation affects areas of medium ES quality, then different measures of restoration (mitigation or compensation) were promoted [43][44][45].
At that stage, the cause-effect relation between input and output of each single ES model has been verified to avoid mistakes and the potential to distort the evaluations of maps [46,47]. While the spatial distribution of the majority of ES was directly understandable (e.g., habitat quality, carbon sequestration, crop production, crop pollination, and water yield), the nutrient retention and sediment retention models presented a distribution dependent from the DEM, because the final output depends on gravitational models, vegetation, and soil characteristics [34,36,38,40,48,49].
The analytical interpretation of such ES with GIS tools was sometimes counterintuitive and the DIST research team was critical and skeptical about the utilization of models for decision-making processes. Indeed, the process concerning the interpretation of output that determines the design of land use rules has been considered critical.
Considering the above-mentioned perplexities, a period of validation for the nutrient retention model has been planned with the aim of finding accurate relations between input and output; particularly considering the influence of DEM and LULC [50]. Validation consisted of (i) the model being tested by changing the DEM and verifying its effect on the output; and (ii) the output was then used as a benchmark to verify the interaction with LULC.
The methodology of validation considered the analysis of intermediate files as a product of each single model operation. Such an approach was useful to better understand how the model relates the input to the output at each step, increasing the feedback with the final results, and simplified the judgment. Each process was then discussed by the group with technical brainstorming increasing the confidence until a final validation was reached, agreed, and approved.
According to the InVEST users guide [28], the selected inputs for nutrient retention are (see Table A1): (1) LULC. A GIS raster dataset, with an integer LULC code for each cell.
(2) The analysis has been carried using the Regional Land Cover raster dataset of 2010 (Land Cover Piemonte (LCP)) with a scale of representation of 1:10,000. This dataset is of high utility and precision. Its contents directly fit with the land use classes of the Corine Land Cover Pan-European datasets. Thus, its utilization guarantees a cross-checking analysis of input data between different scales of representation. The LCP has been used as LULC input data for calculation using the selection of the squared territorial borders.  Table. A .csv table of LULC classes containing data on the water quality coefficients used.
Inputs 3 (root restricting layer depth), 4 (precipitation), and 5 (plant-available water content) are dependent to the land capability classification map of Piedmont ( Figure 2) which covers the entire catchment area and defines the soil characteristics. The nutrient retention model assumes that nutrients moves towards streams, that the soil is saturated by water (rainfall), and run-off is then generated. The distribution of land capability classes is quite heterogeneous in the mapping area. Classes range from 1 (soil suitable for agricultural purposes) to 4 (soil with several limitations to agricultural uses). The soil depth and its composition (the mixture of sand, silt and clay) was a proxy to set the plantavailable water content using SPAW software (freely downloadable from the United States Department of Agriculture website) which was needed to set the values of input 6 ( Table 2).

Spatial Distribution of Flows and Infiltrations
The model output tracks the path load of a diffuse contaminant (nitrogen, by default), which reaches the stream level from the loading coefficient associated with land use. Accordingly, a sum of per-pixel values in a watershed basin gives, as output, the total amount of pollutants that flow into water streams, affecting the overall water quality. Nitrogen is assumed by the model as the nutrient due to agricultural practices which normally largely affect the water quality at the catchment level. Loading values are concentrated on class 2 (agricultural areas, see Table A1), while retention parameters are distributed in class 3 (natural and semi-natural areas).
The more planned land use changes increase the values of filtering instead of loading and the greater their interaction among the catchment area, the more the absolute nutrient retention value decreases, as a consequence of a better land use composition in the scenario representation.
Interpretation of results comes from an analysis of the "n_export.tif" pixel distribution which represents the biophysical output of the model. The reliability of output, especially for comparative studies between alternative land uses, depends on the interaction between the input data. Even if the DEM sensitivity remains obscure most of the InVEST users, the model generates intermediate folders where each process is placed: • The model elaborates a water yield calculation based on an evapotranspiration model. This procedure allows the calculation of an intermediate per-pixel representation where the fraction of the rainfall water that is not stored in the soil and removed by evapotranspiration processes moves to another pixel; • When water is not retained by pixels, the run-off process starts to track the flow-paths that interact with DEM defining the up-flows and stream areas by gravitation. The run-off does not consider if the water routes move via the surface, subsurface, or base flows; it just interacts with the land use model and related nutrient load/retention values. The balance among the load and retention indicates where the landscape contributes to filtering and retaining the nutrients and where nutrients reach the streams.
A preliminary evaluation of the two intermediate results is to check if the runoff model generates the streams in the area where the land use is composed of water bodies. The more streams produced by the runoff model that overlap the water bodies in the land use, the more the results will be reliable in the way the nutrient is routed to the contamination points.
This basilar condition is crucial to set the command "threshold flow accumulation" which defines the number of upstream cells that must flow into a cell before it is considered part of a stream. Indeed, if the model is used at the local scale (1:5,000/1:10,000), then the command has to be set with values ranging from 50 to 200, otherwise it will be sufficient to use values ranging from 500 to 1000 and the correspondence between streams and the water bodies of the large-scale land-use dataset will be achieved.
As explained, DEM is essential to route the flows of nutrients. Its reliability depends on (i) the geometrical definition; and (ii) the continuity and the distribution of the values (holes or gaps to fill in the DEM are frequent).
In the case study, a preliminary observation of the intermediate runoff model has been done using a grid DEM of a 10 m cell (scale of representation, 1:10,000, year 2005) covering the entire territory of the catchment area. It was expected that the flows of water routes down the slopes follow trajectories parallel to gravitational lines, while the model routed the flows with some perpendicular trajectories. Therefore, it has been decided to change the original DEM with an upgraded version which uses a different digital source and presents a denser pixel distribution (which was fundamental to mesh the values among the hill slopes). The geometrical precision of the LULC map was not considered because (i) it did not interact with intermediate

Model Interaction with DEM and LULC
The correction of the DEM determined a new runoff model which was validated by the group, and the nutrient retention model has been analyzed using its final output instead of the intermediate.
The DEM used in the new version is a product of the official topographical survey of the Piedmont region. It was created between 2000 and 2010 by a grid model of 25 m; thus, the pixel dimension is coarser than the previous DEM used, while the results seem to be more acceptable. The integration of the LULC map with the DEM determines the size of the output grid cells.
The evaluation of the file n_export.tif shows how the integration of the different intermediate results interact together ( Figure 5): where the run-off process in upstream areas is intense, then the load of nutrient is higher (particularly when no natural vegetation acts as a filter to the nutrients). In that case, the model shows how pollutants that are not removed by vegetation contaminates the water bodies. The use of the 25 m DEM generates a result which was closer to what was expected: the nutrient concentration happened along the buffer areas of rivers, streams, and natural water bodies following a gravitational route and infiltrates where the interaction between the DEM, soil properties, and LULC determines water sinks. Even without a quantitative assessment, the better output distribution is achieved with a coarser DEM resolution, meaning that sometimes the grid resolution does not correspond to a better model output [38] Once the distribution of the model was spatially acceptable, then a quantitative test was applied to understand how the interaction with LULC properties changes the final result. As introduced, the assessment of the load and retained absolute values of nitrogen associated with each land use type depends on the user. Therefore, to test if the model responds to a different configuration of LULC data, it has been decided to shift the LULC-associated values of the nutrient load with values of retention, thus, obtaining an output which is generated by a different configuration of LULCdependent variables (Table A2). The result (Table 3) shows a considerably different distribution of nutrients. The interaction between the DEM-dependent block (runoff) and LULC-dependent block (nutrient loads/retain) generates a new map where the contamination process happened from the upstream to the downstream area, but with a different distribution and quantity of values (as expected).
The test confirmed the response of the model to the user modification of input, generating a radically different output. Considering the "n_exp_tot" is the biophysical output that is considered during the planning process, the differences between the input configurations are outlined by the rows "diff 1-2" and "diff 2-3" (Table 4). Particularly, the correction of the DEM generates at the catchment scale an absolute augmentation of nitrogen exported to the stream equal to more than 222,658 kg, which corresponds to an augmentation of the previously-registered value of 17.9%. This measure reveals that the DEM quality influences the model workflow and sheds light on the needs of careful use of the input data, assessing both their spatial distribution and their absolute value.
While considering the "diff 2-3", it is shown that when loading and absorbing values are modified, the net absolute decrease of nitrogen exported to the stream is equal to more than 1,347,448 kg, which implies a reduction of infiltration equal to 91.8%. The decline is due to the LULC conversion which reduces the number of load areas in upslope territories and increases the filtering ones in plain zones; then interaction of runoff with LULC distribution generates an abatement of nitrogen infiltration. That configuration shows a total amount of nitrogen infiltrated of 120,150 kg in the catchment area. The reduction confirms the dependency of output from LULC associated data: the agricultural areas are double the natural and semi-natural ones (121,726.90 ha compared to 65,798.17 ha), therefore, the shift in LULC properties means that the nitrogen loads are halved while the retention areas are doubled. This result seems to confirm good feedback between the quantitative values associated with the input data and their interaction to generate the final output.
The collaboration with external expertise (Dario Masante, Joint Research Centre) confirmed that the intermediate flows of the model were analyzed by other relevant ES software (AIRES of LUCI) [5]. Outputs of the validation process were shared with the external research group, the ones that tested the nutrient retention function with different software to define a comparative assessment.

Model Sensitivity to DEM
While the feedback from LULC changes was quite intuitive and connected to a linear relationship between the quantity of a particular LULC class and its value of nutrient load or retain, the response of DEM changes still remained less understandable and requested an in-depth analysis.
Four raster DEM datasets were downloaded from different sources (national and regional web geodatabases) and used to see how the nutrient retention output changes in response to DEM substitution (Table 4). To apply the analysis, none of the other input data were changed, in order to obtain an evaluation of the DEM-related dependencies. InVEST runs four times in a sub-area characterized by a hilly and mountainous landscape to reduce the time of elaboration for each result.
The analysis of results has been conducted with a global sensitivity approach on DEM models [51,52] using a probability distribution function assigned to DEM variables. A probability distribution has been founded for DEM indices with an uncertainty degree (linear regression), and outputs were tested against the uncertainty of all input variables. DEM data were thought as a model, and probability distribution functions have been assigned to data variables. The variance between input and output was observed as an indicator of the output's sensitivity to input changes [53,54].
A global sensitivity analysis was then conducted, with particular attention to: • Analyze the variation of the model output related to a change in the raster input (DEM). A regression of DEM indices was applied and a simulation of thirteen different DEM configurations and their outputs were tested; • Analyze the correlation and the variance of each DEM index related to the output; and • Analyze the dependencies of the DEM input configuration to the nutrient retention output.
Each raster input (DEM) has been evaluated using its related statistical values. In the case study, the value "average" and "standard deviation" were inversely distributed with the DEM resolution.
As shown in Table 4, even if the large majority of indices are inversely proportional to DEM resolution (except minimum and maximum values), their interaction does not generate an output which is inversely related, too. The value of "n_exp_tot" is not linearly distributed according to the geometrical precision of pixels. Instead, the output seems to demonstrate a stochastic distribution. This condition was already explored in other DEM-dependent models (e.g., landslide prediction models, Penna et al. 39) and, what is curious is that (i) the variability of the output is wide, as values ranges from 163,392 kg to 10,398 kg, with a standard deviation of 64,271 kg; and (ii) the lower value (DEM 25) is the one that fits with the real potential distribution of nutrients in the territory, according to an expert-based analysis of the pixel distribution. The geometrical distribution of DEM 25 gives a better feedback compared to others that generate mistakes in the intermediate flow analysis (run-off index distribution).
This means that it is not true that a finer geometrical input resolution in ES models will generate a reliable output. The analysis of the different outputs shows how the change in DEM generates a water-routing intermediate model with different sizes, paths, lengths, and directions of streams in the catchment area with a better approximation to the real streams' distribution provided by DEM 25. To have a better analytical assessment of the data, the relation between the model output and DEM input configuration was analyzed with a regression analysis that simulated eleven DEMs and relative outputs. The DEM resolution of simulated values ranged from five to 75 cells, including the four cases studied (Figure 6).
Simulations were provided using the index "sum" as the regressor because it showed better correlation to the model output. They confirm that the distribution of the output is nonlinear, and not related to the DEM resolution. The "confidence range" (95th percentile) included the result of the model in just one case (DEM 10) showing, in two cases, an underestimation (DEM 20 and 75) and, in one case, an overestimation (DEM 25). The standard deviation for the simulated models was still high (23,500 kg). The greater deviation is for DEM 25 and, secondly, for DEM 75 (while DEM 20 has a lower deviation). An in-depth analysis of the values shows that the standard deviation for registered values (model output) and simulated ones remains high (respectively, 64,271.43 and 23,500.93). The confidence range is high, too. Observing how the variance (standard deviation) affects the model output it is possible to test the uncertainty of the input values.
An average influence value for the simulated case is of 26.54%, ranging from 15.27% (DEM 10) to 35.62% (DEM 70). Thus, the variability of the nutrient retention output is affected by a degree of uncertainty dependent on the DEM resolution, which should be considered high (more than 25%, on average).
In conclusion, the interaction of the nutrient retention model with different DEM resolutions is not predictable. Accordingly, it is necessary to have a visual feedback of the pixel output distribution to check how nutrients move toward the landscape.
The variability of the output, which is influenced by the DEM characteristics, is very high, thus an empirical test of the model is highly suggested. In this case, a field study measurement was not considered; nonetheless, to achieve the best input calibration, an empirical survey should be a better solution, with an analysis of the presence of nitrates into the streams in the catchment area.
Even if, for the large majority of users, the DEM resolution is considered a proxy of its reliability (here intended as a good approximation to the real morphological condition of the territory), its interaction with other elements (such as the land use configuration rather than the geological characterization of the soil) is less known and not adequately understood during the use of the InVEST software.

Research Limitations and Questions
New models, such as the nutrient retention of InVEST, are helping planners to identify areas where ES are delivered, thus obtaining real progress in the definition of planning dispositions and regulations aimed at increasing the natural capital of the territory [55,56]. Nonetheless, this study sheds light on several important limitations and new questions that arise when a new technology supports the decision-making process for planning purposes. Limitations are important considering the raising awareness of planners or scientists to ES modeling.
The main issues are: (i) the technical limitations in the use of the model; (ii) the restriction in the interpretation of the outputs. With respect to the first, the model has some limitation: (i) it is less applicable to locations where the hydrology is determined by rainfall intensity or in areas where flash rains are predominant; (ii) the model does not address any chemical or biological interactions that may occur from the point of loading to the point of interest; and (iii) the model assumes that there is continuity in the hydraulic flow path [28].
Considering the second issue, even if the study conducted resulted in a better comprehension of the model mechanisms, quantitative outputs are not evaluated by a field campaign of real measures and the model calibration is still incomplete. Such limitations are considerably negative if the model is used to suggest where new environmental compensation measures are finalized to plant new species in buffer areas. If the amount of nitrogen is just indicative than the extension of the buffering system, the selection of species, and the overall biodiversity offset assessment will be complicated [57].
Moreover, the qualitative evaluation of nitrate concentration in open water (50 mg/L, European Directive n. 676 of 12 December 1991), cannot be properly addressed. Despite that, in the near future the model should be tested extending the analysis to the wider ideological catchment area where the basin authority knows the concentrations of nutrients at pre-defined points. In that case, when the user obtains a good preliminary distribution of data, then a better calibration should be achieved.
Nevertheless, it was commonly decided to use the output just as a baseline threshold rather than an absolute value: when an alteration in land use is detected, the model ran a scenario analysis, which should be compared with the original threshold. The rate of change represents the better or worsening delivery of the specific ES.
Finally, one of the major limitations of the nutrient retention model is the complexity of designing each input to meet the data type for the entire catchment area, especially when the model is used, as in the case of SAM4CP research, at the local scale with a fine resolution. Additionally, the temporal resolution of data will affect the response in the catchment to the hydrological integration with respect to land use changes.

The Integration of the Nutrient Retention Model in Planning Activities
Even if results encourage caution and an in-depth evaluation of results, the utilization of the nutrient retention model help planners to directly identify areas where the vegetation capacity of buffering the contaminants should be improved using compensation measures finalized to abate the nutrient loads, or planting species with a specific characteristic of nutrient drainage. By this point of view, it is hard to find specific rules that directly link the use of the model with a planning prescription.
The significant benefit in applying nutrient retention during the decision-making phase is that all stakeholders involved in planning activity (technicians and politicians) are helped to understand where it is possible to retain more nutrients for different land use configurations. Moreover, the runoff distribution, as well as the flow directions, are important maps because they define the hydrological behavior of the considered landscape under normal conditions. When all the land use changes are discussed, considering the nutrient retention model, it is possible to address how the changes in soil cover and soil properties increase or decrease the resistance to flow from the surface and, consequently, diminish or augment the nutrient loading in the catchment area. Overall, the results encourage the public administrators and planners to introduce such methodologies into their traditional activities.
Nonetheless, more research is needed because the results reported here show that caution has to be used when preparing input data for an InVEST model that is DEM-dependent. The DEM configuration, associated with land use data, should hardly influence the predictive model of nutrient infiltration, in particular, as the model predictive power is shown to be DEM-configuration dependent.

Conclusions
The research conducted shows how it is important to validate quantitative analysis for specific ES in order to overcome the gaps that separate the theoretical framework on ES and its practical application.
Even if at its preliminary stage, the achieved results are encouraging and explain that planners should be aware of how to use techniques of analysis to integrate ES in their valuations. The new locations for ecological compensation zones in the three case studies were selected according to the nutrient retention capacity for specific areas. In particular, the slight reduction of arable land for new natural buffering zones and the selection of land use transformations with the lower environmental impact was supported by the use of InVEST outputs.
Increasing performances of nutrient retention, implies that the future application of planning dispositions will guarantee a better capacity of mitigating the risk of nutrient infiltration in the soil, in the subsoil, and in the water.
Although the traditional procedure for planning and its instruments of evaluation (e.g., the strategic environmental assessment) do not consider ES as quantitative values for measuring the environmental sustainability of plans and programs, the results of the study confirm that a quantification and evaluation of specific ES is, nowadays, possible and welcomed.
Acknowledgments: This study constituted an autonomous work conducted by the DIST group of research involved in the LIFE+ SAM4CP project. The group was supported in the study by Gabriele Garnero (DIST) who is an expert in geomatics and GIS applications for planning. Authors would like to also thank Dario Masante for his valuable feedback on nutrient retention validation procedures and Franco Pellerey for his willingness to work on the sensitivity analysis.
Author Contributions: All authors equally contributed to the article.

Conflicts of Interest:
The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.