A methodology for assessing spatio-temporal dynamics of flood regulating services

The effects of land use alteration, migration and urbanization are key aspects in flood management, as human activities can strongly influence the capacity of ecosystems to provide flood regulating ecosystem services and determine their demand. This study analyzes spatio – temporal dynamics of flood regulating ecosystem services to support watershed management planning. A methodology for mapping the supply and demand of flood regulation is proposed and applied to the Arno River basin, in central Italy. The spatial explicit analysis of flood regulating ecosystem services supply is carried out with SWAT - Soil and Water Assessment Tool, whose outputs are synthetized by two indicators to evaluate the retention capacity of each land use class originating from CORINE data sets. Quantification of demand for flood regulating ecosystem services is based on flood hazard classes derived from the existing local flood management plans (i.e., PAI-Piano per l ’ Assetto Idrogeologico and PGRA-Piano di Gestione del Rischio Alluvioni). Supply and demand data are then combined to obtain budget maps of flood regulating ecosystem services and their evolution, between 1990 and 2018. The results show how both demand and supply of ecosystem services have changed in the last decades, highlighting the main hotspots at the catchment and sub-catchment scales. With the increasing urbanization, the demand values have grown in the Arno floodplains, where residential, industrial and commercial zones are located. At the same time, land use changes have altered the water regulation supply, resulting in a generalized decrease of the basin capacity to provide flood regulation services. The maps and tables obtained show the fundamental role of forest and other vegetated areas whose protection is a priority to assure future flood regulation and associated co-benefits (e.g., regulation of air quality, reduction of erosion, improvement of water quality, wood fuel). The assessment of flood regulating here proposed is a powerful tool for decision makers to improve flood regulation and provides a sound base of knowledge to identify and locate flood prevention and mitigation measures.


Introduction
Floods are among the most affecting natural hazards worldwide, causing loss of lives and damage to buildings, industrial settlements, communication routes and agricultural areas (CRED and UNISDR, 2015;Schanze et al., 2007;Dottori et al., 2018).According to worldwide statistics, flood damages are continuously increasing (Schanze et al., 2007;CRED and UNISDR, 2015) and this trend is exacerbated by anthropogenic climate change and human interventions on river morphology and land use (Bronstert et al., 2002;Pall et al., 2011;Ekka et al., 2020).These effects, combined with population migration and urbanization, continually alter flood exposure and vulnerability putting extra pressure on physical and social infrastructure (Eigenbrod et al., 2011;DePaul, 2012;Mazzoleni et al., 2020;Akhter et al., 2021), as confirmed by numerous disastrous flood events that have occurred in recent decades despite huge flood risk reduction efforts (Kundzewicz et al., 2018).
In this context, ecosystem-based adaptation approaches can provide an interesting perspective on flood management, aiming at promoting a more sustainable interaction between human and natural systems and enhancing climate change coping capacity (Chong, 2014;Doswald et al., 2014).Indeed, well-functioning ecosystems have the potential to buffer the impacts of climate change and to provide multiple societal benefits, namely ecosystem services (Millennium Ecosystem Assessment, 2005;Iacob et al., 2014).In particular, ecosystems can provide flood regulating services, preventing the rapid runoff of surface water and reducing peak discharges (Bayley, 1995;Vallecillo et al., 2020).This capacity is strongly influenced by land use management that affects both the supply (altering the hydrological cycle) and the demand, driven by the socioeconomic system (Vandecasteele et al., 2018;Stammel et al., 2021).Therefore, it becomes crucial to set up management strategies that can positively influence the environment capacity to contribute to flood mitigation and prevention.This includes the assessment of flood regulating services that can be analyzed both in terms of supply (capacity of the environment to provide the service) and demand (human driven request of that service) (Stürck et al., 2015;Wamsler et al., 2016;Halbe et al., 2018).
The concept of ecosystem services (ES) has gained extraordinary popularity in the recent years (Seppelt et al., 2011;Zhang et al., 2019) and has been accompanied by the development of different methodologies to support the understanding of the complex relationships between ES and the socio-ecological processes underlying them (Egoh et al., 2008;Crossman et al., 2013;Burkhard & Maes, 2017;Santos-Martín et al., 2019;Pacetti et al., 2020a).Several studies have focused specifically on flood regulating ES, trying to value and map the biophysical features that determine flood regulation at different scales (continental, national, watershed or city scales) and to estimate the associated service demand (Egoh et al., 2008;Syrbe & Walz, 2012).Stürck et al. (2014) performed a qualitatively assessment of spatial patterns of ES supply and demand.Shen et al. (2019) focused on the integration of supply and demand analysis at a single city scale.Only a few studies have analyzed the evolution in time of ES: Stürck et al. (2015) using intermediate-complexity ES models determined the supply and demand evolution at the continental scale, while Li et al. (2019) used the high-complexity model CLM-GBHM to investigate flood regulating service supply under different management scenarios, without evaluating the associated ES demand.Nedkov and Burkhard (2012) based their analysis on the ES matrix approach proposed by Burkhard et al. (2009) that provides a flexible approach based on lookup table that are used to link ecosystem types or other geospatial units (e.g., Land Use and Land Cover types-LULC) with ES.In their study, they introduced the use of the basin-scale hydrological model (i.e., KINEROS and AGWA tool) for the quantification of ES supply, coupling the analysis with demand mapping to determine the budgets for flood regulation.Their analysis highlighted the advantages of using hydrological models to investigate the link between hydrological processes and the ES production but also stressed the limitation of using an event-based model that represents mainly runoff (overlooking important element as evapotranspiration or infiltration).The oversimplification in representing catchment characteristics (e.g., flood frequencies, peak flows or flood durations using statistical methods or simplified hydrological models) can strongly affect the final ES assessment, as highlighted by studies focused on specific aspects, such as the effects of land use change or climate variability (Bagstad et al., 2013;Stürck et al., 2014;Stürck et al., 2015).Boyanova et al. (2016), Boyanova et al. (2017) used the Soil Water Assessment Tool (SWAT, Arnold et al., 1998) hydrological model for the quantitative assessment of flood regulating ES showing its suitability to overcome shortcomings of other models (Nedkov and Burkhard 2012;Boyanova et al. 2014), but also stressing the model dependency on a large amount of data, whose quality is fundamental.
Building up from previous studies, this research contributes to improving understanding and spotting limitations of flood regulating ES by offering spatio-temporal dynamics of ES.This approach provides a useful insight on flood regulation as one of multiple ecosystem services that can operationalize ecosystem-based approach and support decision making in flood management (Iacob et al., 2014).
In particular, taking advantage of the SWAT model capability to provide a spatially explicit analysis of the hydrological process underpinning regulating ES, this study focuses on the spatio-temporal dynamics of flood regulating ES supply and demand over the past thirty years (from 1990 to 2018) in the Arno River basin (central Italy).

Case study
Italy is characterized by frequent flooding occurring due to intense and localized rainfall combined with complex characteristics of the territory.Major flood events affect many people and cause serious economic damage (Ministry for the Environment, 2007).In this context, the Arno River basin, located in Central Italy, appears as an exemplary case study, with a long track of major recorded flooding events: in 1333, 1557, 1589, 1844and 1966(Galloway, et al., 2020)).After the dramatic flood of 1966, several interventions have been implemented to improve flood regulation with reforestation, stabilization of slopes, and construction of weirs (Caporali et al., 2005).However, the land use changes over the past decades caused environmental degradation and negative effects on water regulation due to several factors, among which: 1) increase of urban areas, with the reduction of permeable surfaces; 2) reduction of the agriculture productive potential, resulting in degradation of the territory, and 3) internal migration from rural to urban areas (IRPET, 2010).Indeed, the exposure to flooding, the complex and varied lithology (Baiocco et al., 2003), the active land use change and recent population dynamics make the Arno River the ideal case for analyzing spatio-temporal analysis of supply and demand budgets of regulating flood services.
The Arno River basin, according to the Köppen and Geiger climate classification (Rubel et al., 2017) is located in the Csa category, i.e., temperate climate zone, characterized by hot and dry summers.The analysis proposed here focuses on the portion of the river basin upstream of the city of Florence (closing section corresponding to the Ponte a Signa gauge station), including 7 river sub-basins (i.e., Chiana, Casentino, Upper Valdarno, Greve, Medium Valdarno, Sieve and Bisenzio; Fig. 1) and covers an area of approximately 4350 km 2 .

Materials and method
Supply, demand and budgets of flood regulating ES are quantified and mapped following the conceptual scheme shown in Fig. 2. The SWAT model is calibrated and validated over the years 2003-2014 and used to simulate the evolution of catchment behavior due to LULC changes from 1990 to 2018.The model output is converted into two indicators to quantify ES supply, resulting in four maps (one per year, in 1990, 2000, 2012 and 2018) for each indicator.
Two demand maps are derived from existing flood management plans, which contain the identification of flood hazard levels.The Hydrogeological Plan (PAI -Piano per l'Assetto Idrogeologico) is used to obtain flood regulating demand quantification for the years 1990 and 2000 and the Flood Risk Management Plan (PGRA -Piano di Gestione del Rischio Alluvioni), which replaced the previous plan in compliance with the Flood Directive 2007/60/EU, is used for the years 2012 and 2018.Finally, supply and demand data are combined to obtain budget maps of flood regulating ES and assess their spatio-temporal analysis.

Hydrological modelling
The assessment of flood regulating services supply is based on the evaluation of the hydrological processes occurring in the catchment.The hydrological understanding of the area is processed with SWAT, a physically based, semi-distributed and continuous time model that allows several different physical processes to be simulated in a watershed (e.g., hydrology, climate, soils, land management, plant growth, pesticides and the nutrients cycle).Water balance equation is the driving force behind the hydrological processes occurring in the basin: where all the quantities on day i are in mm of water; SW t is the final soil water content; SW o is the soil water content at the beginning of the analysis period t (days); R day is the amount of precipitation; Q surf is the amount of surface runoff; E a is the amount of evapotranspiration; W seep is the amount of water entering the vadose zone from the soil profile; and Q gw is the amount of return flow.
The equation is applied at the level of Hydrological Response Units (HRUs), the smallest spatial units of the model, representing the unique combinations of land use, soil and morphology characteristics.Model

Table 1
Input variables, source and relative scale used to set up the hydrological model.input data (Table 1) were mainly derived from the Tuscany Region (Tuscany Region, 2020).A DEM with 10 m resolution is used identify the drainage patterns and 180 soil types (see Supplementary material) are identified based on their composition, while 29 land cover classes (see Supplementary material) are derived from the CORINE land cover dataset and then associated with corresponding SWAT land use/plant classes.The soil parameters are calculated using the Pedo Transfer Function (PTF), developed by Saxton and Rawls (2006) while the land cover classes are associated with the similar SWAT database land use/ plant classes.Rainfall and temperature data are extracted from the Tuscany Region dataset (see Supplementary material).While relative humidity, solar radiation and wind speed data of 15 points (see Supplementary material) are acquired from the global climate database of the National Centers for Environmental Prediction (NCEP) -Climate Forecast System Reanalysis (CFSR).
The model is calibrated and validated for the period from 1 January 2000 to 31 July 2014, adopting the CORINE Land Cover map of the year 2012 as baseline.Then, keeping all the other variables unchanged, the model is run three more times using the land cover maps of 1990, 2000 and 2018 to evaluate the effects of land use change.

Capacity assessment of flood regulating ecosystem service supply
Two flood regulating supply indicators are elaborated based on the output provided by the SWAT model for the years 1990, 2000, 2012 and 2018.The first indicator (I 1 ) is built according to the methodology proposed by Boyanova et al. (2017), extracting the average annual value of 5 main hydrological variables (i.e., Evapotranspiration (ET), surface runoff (SURFQ), lateral flow (LATQ), percolation (PERC), groundwater flow (GWQ)) for any land use scenario and interpreting them in terms of their contribution to flood mitigation (Fig. 3).Higher values of ET, PERC, GWQ, LATQ correspond to higher values of supply as they contribute to generate a delay on the peak discharge.Indeed, ET and PERC represent water flows that do not contribute to the streamflow, while GWQ and LATQ represent the portion of streamflow that slow down the peak propagation.Higher values of SURFQ correspond to lower level of service supply, as it represents the portion of water that directly contributes to the streamflow.Then, the variable distribution among the different land uses is rescaled based on the quartile function assigning semi-quantitative scores, ranging from 0 to 3 (0 = no relevant supply, 1 = low supply, 2 = moderate supply, 3 = high supply).The final land use class is obtained by averaging all values again and then applying the quartile function.
The second indicator (I 2 ) is obtained computing the ratio between the surface runoff and the relative precipitation variables for each HRU; the average annual value is assessed for each land-use class, following Castelli et al. (2017).
The range of the four classes, from 0 to 3, are again obtained with the quartile function, by arranging the I 2 values from maximum to minimum.
By linking the 0-3 values obtained with the two previous indicators to spatial data in GIS, estimates of ecosystem service supply can be assessed and transferred to different spatial and temporal scales.

Assessment of demands for flood regulating ecosystem services
The assessment of flood regulating ES is based on hazard maps.Flood hazard maps identify areas potentially affected by floods, according to pre-established scenarios, indicating, where possible, information on discharge, water level and flow velocity.Hazard is defined as the likelihood of occurrence of a flood with specific attribute at a given location in a reference interval of time (Wright, 2015).Flood hazard magnitude can be expressed in probabilistic recurrence interval, known as return period (T) which is an average time or an estimated average time between the occurrence of two events of similar magnitude.The highest hazard areas will be in the floodplains, where urban, industrial and commercial areas are located.Therefore, they can be associated to the highest demand values for flood regulating ES.While lower demand rates can be likely expected in areas of lower population density or industrial assets.In this study, demand quantification is derived by existing flood management plans.The enactment of Law n. 183 in 1989 reordered the legislation governing the management and functions of soil protection and in 1996 the Arno River Basin Authority adopted a Provisional Plan (Piano Stralcio) for hydraulic risk mitigation.In 2000, in the accomplishment of the Italian Law n. 180 (1998), called Sarno Law, published after the catastrophic mudflows that occurred in southern Italy on May 5, 1998, the Arno River Basin Authority adopted an Extraordinary Plan (Piano Straordinario).Then, the Hydrogeological Plan (PAI -Piano per l'Assetto Idrogeologico) was approved in 2002.It contains the identification and the perimeter of four hydraulic hazard classes (PI1, PI2, PI3 and PI4) for the Arno River basin.From the experience gained with the PAI, the Arno River Basin Authority, within the application of the Flood Directive 2007/60/EU produced the hazard maps of the Flood Risk Management Plan (PGRA -Piano di Gestione del Rischio Alluvioni), accounting for landscape evolution and effects of implemented interventions.The PGRA of the Arno River Basin Authority was adopted by the deliberations n. 231 and 232 taken on December 17, 2015 and has been finally approved by the deliberation n. 235 (March 3, 2016).The last two plans have used same assessment criteria; however, the PGRA hazard maps present a greater level of detail thanks to more precise input data and more powerful and new calculation tools.Moreover, PGRA defines three, rather than four, hazard levels (P1, P2 and P3) using different return period thresholds.
In this study, the PAI hazard maps is used for assessing demand for flood regulating ES in the years 1990 and 2000, while the PGRA maps is used for the years 2012 and 2018.In order to obtain homogenous results, as can be seen in the Table 2, maps from the two management plans are compared and four classes identified, assuming that the highest hazard level areas will have the highest demand, according to the following scale: 0 = no relevant demand, 1 = low demand, 2 = moderate demand, 3 = high demand.

Model calibration
SWAT model is calibrated, from 01/01/2003 to 12/31/2011, using manual calibration procedures and validated, from 01/01/2012 to 07/ 31/2014, using statistical and graphical techniques (Pacetti et al., 2020b).According to Moriasi, et al. (2007), streamflow model simulation can be judged as satisfactory if the Nash-Sutcliffe efficiency is NSE > 0.50, the ratio of the root mean square error to the standard deviation of measured data is RSR ≤ 0.70 and if percent bias is PBIAS ± 25%.Daily discharge data at the Ponte a Signa gauge station (Fig. 1) are retrieved form Tuscany Region database, from 1 January 2003 to 31 July 2014.These measured data are compared with the corresponding simulated daily streamflow (Fig. 4).
Following the procedure described in Shrestha (2017) and Krause et al. (2005) five parameters are selected and modified for calibration (Table 3).
The statistics (NSE RSR, PBIAS) show that there is a good agreement between the measured and simulated flows both for calibration and validation periods (Table 4).

Mapping of flood regulation supply
The supply capacities of different CORINE land cover classes are defined (Table 5) with the two indicators described above and mapped with GIS tools for four different years, i.e.,1990years, i.e., , 2000years, i.e., , 2012years, i.e., and 2018 (Fig. 5 (Fig. 5).Then, starting from the 1990 map, the variations over time at the sub-basins scale are evaluated (Fig. 6).
Both indicators show that forests (FRST, FRSE) and densely vegetated areas (SHRB) have high flood regulating capacities.While urban (URBN) areas due to their impervious surfaces and barren (BARR) areas without vegetation show low supply capacity.Even the agricultural (AGRL) areas due to the many are characterized by sparse vegetation and a low soil retention capacity.The major differences between the two indicators are for two land uses: Orchard (ORCD) and Cropland/ Woodland Mosaic (CRWO).However, at the basin scale the results are similar.There is a continuous reduction in the flood regulation capacity from 1990 to 2012 due to the expansion of urbanization and intensive agriculture.While the 2018 map shows a slight increase in ES supply capacity with the increase in wooded and vegetated areas in the countryside (Table 6).These results are confirmed by the fact that, using the calibrated hydrological model, at the Ponte a Signa river stage section, there is a continuous increase in the simulated daily water flow for the first three CLC maps and a slight reduction with the 2018 map.Between the seven sub-basins there is a trend that varies over time (Table 6).With both indicators, the Chiana River sub-basin has shown the best percentage growth in supply since 1990, while the situation worsens more in the Bisenzio and Sieve River sub-basins.

Mapping demands for flood regulating ecosystem services
As described in the methodology section, two flood regulating services maps are obtained according to PAI (1990 and2000) and PGRA (2012 and 2018) flood management plans available for the Arno River basin.
Both maps, apparently similar as can be seen in the Fig. 7, show that the areas with high and moderate hazard are predominantly along the river streams and in the floodplains where historically the population is concentrated.Consequently, the low hazard is mostly associated with landscapes dominated by agricultural use or large areas of natural vegetation.Compared to the PAI map, in the map obtained from the PGRA there is an increase in demand of about 6% at the river basin scale.In all the river sub-basins there is an increase in demand over time.

Mapping the budget between flood regulation supply and demand
Flood regulating ecosystem service supply and demand data are merged with a spatial overlay to produce a map showing demand-supply balances over the time in the Arno River basin (Fig. 8).Demand classification values are considered negative, thus obtaining a relative scale from 0 (no relevant demand) to − 3 (high demand).While the values of the four supply classes are considered positive from 0 (no relevant supply) to + 3 (high supply).
Then, with a spatial addition operation the final maps are obtained.The supply maps obtained for the years 1990 and 2000 with both indicators are compared with the demand map obtained from the PAI classification.While supply maps obtained for 2012 and 2018 are compared with the demand map obtained from the PGRA.Starting from the 1990 map, the variations over time at the sub-basins scale are evaluated (Fig. 9).As seen in the previous paragraphs, at the sub-basin scale there is a continuous reduction in the water regulation capacity from 1990 to 2012 with a slight increase in 2018 and an increase in the ES demand over the time.Consequently, as shown in Table 7, in the entire study area there is a reduction of the final budget from 1990 to 2012 with a slight trend inversion for 2018.Both indicators show that Chiana and Greve river sub-basins are characterized by increasing final budgets, while other sub-basins show a decreasing percentage trend.In the sub-basin of the Middle Valdarno, which is the most urbanized one, there is the maximum reduction of the budget over time.

Spatio-temporal evolution of flood regulating services
The application of the proposed methodology to the Arno River basin has shown that the growing urbanization over time in the Arno River floodplains has generated an increase in flood regulating service demand.It increased of about 6% at the river basin scale by comparing PAI and PGRA maps.The greatest flood damage occurs in economic centers and urban agglomerations that are concentrated in a very restricted area.The impermeable surfaces of the structures, infrastructures, and services in the urban area with the intensification of agricultural practices in the most fertile and accessible lands simultaneously generated negative effects on the water retention capacity of the soil.However, the supply increased in the most marginal areas not suited to agriculture which over time have been abandoned and where forests and densely vegetated areas are dominant.The spread of woodland systems increased the potential supply of regulating services.At the same time, the decline in semi-natural open landscapes negatively affected the supply of services linked to traditional uses.Both indicators, the I 1 , obtained extracting the average annual value of 5 main hydrological variables for each land use, and the I 2 , obtained computing the ratio between the surface runoff and the relative precipitation for each HRU, show that regions with a high capacity to provide flood regulation are mainly characterized by large patches of natural vegetation (Table 5).At the basin scale, a decreasing trend in the final ES budget is observed (Table 7) from 1990 to 2012 equal to − 2.28% for the I 1 and equal to − 1.30% for the I 2 .However, there is a trend inversion in the most recent budget (2018) equal to 0.98% for the I 1 and equal to + 0.19% for the I 2 also due to a significant change in land management.With the regional law 65/2014 the Tuscany Region has limited new buildings outside the already urbanized areas.This approach allowed limiting the increase in  waterproof surfaces, reducing soil consumption and slightly increasing the supply capacity of the basin.Even at the sub-basin level, the two indicators have a similar trend.Since 1990, the Chiana river sub-basin has shown the highest percentage growth, while the Greve and Upper Valdarno river sub-basins have shown a slight increase (Table 6).While the Sieve and Bisenzio River sub-basins show the largest percentage decrease in water retention capacity.Similar considerations can be made for the positive values of the budget, but the greatest reduction occurs in the sub-basin of the Medium Valdarno where the city of Florence is located.

Strengths, limitations and outlook of this study
Flood prevention and mitigation measures are matters of growing  importance.It is possible to mitigate the risk of flooding both through structural interventions such as levees and retention basins, and through non-structural interventions such as land planning and management (Petry, 2002).The construction of the first type of measures is frequently associated with detrimental effects on biodiversity and ecosystems (Birkland et al., 2003), therefore the interest in non-structural mitigation measures and the assessment of their flood-regulation capacity has increased (Kundzewicz et al., 2002).Informed land use modifications and changes in land management intensity can enhance flood regulation capacity of catchments areas.Recognizing spatial landscape patterns and interactions between adjacent ecosystems can be a powerful tool to support decision makers in reducing flood risk.Moreover, the ES assessment opens the way to the experimentation of alternative flood management interventions such as nature-based  solutions that can increase the water retention capacity of the soil and then the resistance of environment against extreme rainfall events while supporting a wider set of ecosystem services (Keesstra et al., 2018).Building up from previous studies, this research takes advantage of the SWAT model in order to support watershed management planning and to provide the flood regulating ES.The use of hydrological modelling, which represents the hydrological processes of the environment through simplified mathematical relationships, offers the opportunity to determine the water-related ecosystem services based on the water balance of the study area.The SWAT model was applied in QGIS using the QSWAT extension because the GIS-based models provide the possibility for spatially explicit analyses of output variables previously processed.With this approach, flood regulation demand is assessed and compared to the spatial patterns of a supply indicator.The indicators can also be used to analyze the consequences of historical and projected land use changes on flood regulation services.
However, this approach has some limitations.It quantifies demands solely taking geomorphological properties into account and the land's hazard to floods.While, regulating capacity consider primarily land cover variation without considering, for example, agricultural intensity.The CORINE data set provides satisfactory results, but for more detailed analyses, further data with higher spatial resolution could be used.Moreover, the uncertainties due to anthropogenic climate change were neglected by keeping constant the meteorological data for the period 2000-2014.Climate change involves complex interactions and changing likelihoods of diverse impacts.In recent decades, they have caused impacts on natural and human systems on all continents and across the oceans (IPCC, 2014).
By assessing quality of the input data used for this application, and by looking at further improvement of the methodology here proposed, we suggest to: • Include land management intensity, which could alter the water retention capacity of land cover over time; • Account for potential climate change that can influence the natural processes and the hydrological cycle; • Extend the quantification and mapping approach to other study area and assess multiple ecosystem services (e.g., fresh water, wood production, nutrient loss) to provide a more comprehensive view of the landscape functionality; • Obtain demand maps with other assessment tools such as that of economic damage.In the application shown here, the demand maps accounted only for topographical and hydrological characteristics of the areas; • Use finer input data (e.g., maps with a higher resolution, locally measured meteorological data) to improve the potential of spatial ecosystem service mapping approaches.

Conclusions
Informed land use modifications and changes in land management intensity can enhance flood regulation capacity of water catchments areas.Recognizing spatial landscape patterns and interactions between adjacent ecosystems can be a powerful tool to support decision makers in reducing flood risk.In this study, we presented a methodology to assess budgets of supply and demand of flood regulation ES and to analyze their evolution over the last three decades.The aim of the study was to highlight how landscape management can influence the capacity of the ecosystems to provide regulating services.The analysis was performed for the Arno River basin in central Italy, an area highly exposed to flooding that has faced several Land Use and Land Cover changes over the past years.
The results allowed for the spatial explicit quantification of supply and demand of regulating ES, showing the fundamental role of forest and other vegetated areas whose protection is a priority to assure future flood regulation and associated co-benefits (e.g., regulation of air quality, reduction of erosion, improvement of water quality, wood fuel).Moreover, the results show how socioeconomic development can influence the capacity of the ecosystems to regulate hydrological processes.Land cover changes are linked to socio-economic transformations and agricultural policies.Often economic resources and technical knowledge are not enough to mitigate various risks, whether natural or generated by activities human.
Sustainable territorial planning and integrated resources management are the key elements to sustain human well-being while preserving the ecosystems and biodiversity of the territory.
The present approach allows for relatively simple quantification and mapping of the ES distribution in the watershed and facilitates the economic, environmental and social management of the area.Maps and tables are useful for supporting the understanding of information which otherwise might be difficult to interpret, thus potentially contributing to inform a wide range of decision makers and of stakeholders.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Map of the basin and the sub-basins of the Arno River, central Italy.

Fig. 3 .
Fig. 3. Hydrological variables considered to define the first indicator I 1.

Fig. 4 .
Fig. 4. Observed and simulated discharge data for the calibration period at the Ponte a Signa gauge station.

Fig. 5 .
Fig. 5. ES supply maps obtained with the two indicators for the year 1990.

Fig. 6 .
Fig. 6.ES supply variation maps over the time at the sub-basin level obtained with indicator I 1 in the upper row and indicator I 2 in the lower row.

Fig. 7 .
Fig. 7. ES demand maps obtained from the Hydrogeological Plan -PAI and the Flood Risk Management Plan -PGRA.

Fig. 8 .
Fig. 8. ES budget maps obtained with the two indicators for the year 1990.

Fig. 9 .
Fig. 9. ES budget variation maps over the time at the sub-basin level obtained with indicator I 1 in the upper row and indicator I 2 in the lower row.

Table 2
Demand classification obtained from the Hydrogeological Plan -PAI and the Flood Risk Management Plan -PGRA.

Table 3
Modified parameters for the manual calibration of the model.

Table 5
Supply classes of the different land use types per indicator.

Table 6
Percentage variation over the time of supply in the river sub-basins per indicator.

Table 7
Percentage variation over the time of the budget in the river sub-basins per indicator.
S.Mori et al.