Land Use Change Impact on Flooding Areas: the Case Study of Cervaro Basin (italy)

The main goal of this paper is to study the effect of the spatio-temporal changes of Land Use/Land Cover (LULC) within the hydrologic regime of the Cervaro basin in Southern Italy. selected to produce LULC maps covering a time trend of 28 years. Nine synthetic bands were processed as input data identified as the most effective for the Artificial Neural Network (ANN) classification procedure implemented in this case study. To assess the possible hydrological effects of the detected changes during rainfall events, a physically-based lumped approach for infiltration contribution was adopted within each sub-basin. The results showed an increase in flood peak and a decrease of the rangelands, forests, and bare lands between 1984 and 2011, indicating a good correlation between flooding areas and land use changes, even if it can be considered negligible in basins of large dimensions. These results showed that the impact of land use on the hydrological response is closely related to watershed scale.


Introduction
Safeguarding natural resources is essential in the preservation of the equilibrium of every ecosystem: consistent changes and the destruction of natural resources can cause many problems, including soil erosion and increased runoff.Land cover features are part of the ecosystem equilibrium, and any environmental changes that affect them is worth studying [1]-particularly due to the increasing expansion of urban areas which causes permeable land reduction and increased flooding [2,3].It is widely accepted that the magnitude and frequency of floods is currently increasing [4].This may be ascribed to two different causes: the climate change characterized by the sequence of extreme wet and dry periods, and the impact of land management upon the relationship between extreme climate events and hydrological extremes.Pattison and Lane [5] studied these hypotheses, searching for situations where land management might be influencing flood risk at a scale larger than plot or field.
In recent years, numerous researchers have shown interest in applying geographical information system (GIS) and remote sensing imagery to the extraction of land surface parameters [6], which, if applied to hydrological models, can be useful to obtain reasonable results-especially in ungauged basins [7].Incorporating processed satellite data within hydrological models has become a promising approach for better and accurate quantifications of water resources in river basins [8,9].
There are different methods used to evaluate the relation between land cover changes and runoff increase [10].Many of these use remote sensing images such as the Moderate Resolution Imaging Spectroradiometer (MODIS) [11].Van de Sande et al. [12] compared IKONOS-2 derived maps with EU CORINE Land Cover data set in a two-dimensional hydrodynamic model, simulating floodplain inundation for the southern part of The Netherlands.Panahi et al. [13], analyzing MODIS images to understand the effect of land use change on the hydrologic regime of the Madarsu basin (Iran), found that runoff was 10 times higher after 40 years, due to agricultural lands increase and a decrease in forests and rangelands.Descheemaeker et al. [14] studied the runoff-generating rainfall threshold and found a positive correlation with total vegetation cover: runoff was negligible when vegetation cover exceeds 65%, and it was dependent on soil organic matter, soil bulk density, litter cover, and slope gradient.Investigating how land use was changed, Gholami et al. [15] found that urban and road network development intensified runoff vegetation and flooding hazard within an Iranian forest area during the last 40 years.Results indicated an increased influence of urban development on the volume and peak discharge in sub-basin areas, which had considerable impervious surface development.
While it appears clear that land use is linked to flood risk, it is still to be defined how land use might impact flood trends.Pattison and Lane [5] made a detailed analysis of the connection between land-use and flood increase.They asserted that the exact relationship between land management and flooding risk is more complex, inasmuch as it depends upon local catchment characteristics, the spatio-temporal sequencing of extreme rainfall events, and the structure of the river catchment.However, it is necessary to analyze the ways in which land use management might influence flood generation [16], thought to be a consequence of the effect of land cover management on the partitioning of rainfall between runoff and subsurface flow, the watershed storage effect during a rainfall event, and the way water is transferred from hillslopes to the river network.
It is a commonplace to think that change in land use due to anthropic action results in increased flood frequency and severity.May this be true in all cases?Bloschl et al. [4] hypothesized the impact of land use and climate variability on hydrological response as a function of scale, as summarized in Figure 1.Abrupt changes in watershed response can occur as a result of land use change and can be particularly severe at small scales; however, there are examples of tremendous large-scale changes, such as the increased nitrate concentrations of the Danube as a result of political changes in eastern Europe in the 1990s [4].inundation for the southern part of The Netherlands.Panahi et al. [13], analyzing MODIS images to understand the effect of land use change on the hydrologic regime of the Madarsu basin (Iran), found that runoff was 10 times higher after 40 years, due to agricultural lands increase and a decrease in forests and rangelands.Descheemaeker et al. [14] studied the runoff-generating rainfall threshold and found a positive correlation with total vegetation cover: runoff was negligible when vegetation cover exceeds 65%, and it was dependent on soil organic matter, soil bulk density, litter cover, and slope gradient.While it appears clear that land use is linked to flood risk, it is still to be defined how land use might impact flood trends.Pattison and Lane [5] made a detailed analysis of the connection between land-use and flood increase.They asserted that the exact relationship between land management and flooding risk is more complex, inasmuch as it depends upon local catchment characteristics, the spatio-temporal sequencing of extreme rainfall events, and the structure of the river catchment.However, it is necessary to analyze the ways in which land use management might influence flood generation [16], thought to be a consequence of the effect of land cover management on the partitioning of rainfall between runoff and subsurface flow, the watershed storage effect during a rainfall event, and the way water is transferred from hillslopes to the river network.
It is a commonplace to think that change in land use due to anthropic action results in increased flood frequency and severity.May this be true in all cases?Bloschl et al. [4] hypothesized the impact of land use and climate variability on hydrological response as a function of scale, as summarized in Figure 1.Abrupt changes in watershed response can occur as a result of land use change and can be particularly severe at small scales; however, there are examples of tremendous large-scale changes, such as the increased nitrate concentrations of the Danube as a result of political changes in eastern Europe in the 1990s [4].Chen et al. [17] analyzed the effect of urbanization on flood characteristics at catchment levels for a Chinese catchment affected by an intensive urban expansion over the last twenty years.The results showed a correlation between spatial variations of urbanization and rainfall depth, but increasing urban area did not necessarily produce more runoff.This indicates that runoff generation is influenced by several factors, such as land use, soil, terrain, weather conditions, and the groundwater system [18].
There are two main approaches to analyzing the impact of land use changes on flooding [19]: an upward approach, which consists of model cascades where each model represents a sub-process such as rainfall, subsurface flow, etc.; and a downward approach, which consists of a trend analysis of long runoff data and catchment studies.
For the purposes of the case study presented in this paper, the second approach was chosen in order to investigate the Cervaro watershed in Southern Italy.Firstly, the study area is described, then the contribution of remote sensing in generating Land use/Land cover (LULC) maps from multi-date Chen et al. [17] analyzed the effect of urbanization on flood characteristics at catchment levels for a Chinese catchment affected by an intensive urban expansion over the last twenty years.The results showed a correlation between spatial variations of urbanization and rainfall depth, but increasing urban area did not necessarily produce more runoff.This indicates that runoff generation is influenced by several factors, such as land use, soil, terrain, weather conditions, and the groundwater system [18].
There are two main approaches to analyzing the impact of land use changes on flooding [19]: an upward approach, which consists of model cascades where each model represents a sub-process such as rainfall, subsurface flow, etc.; and a downward approach, which consists of a trend analysis of long runoff data and catchment studies.
For the purposes of the case study presented in this paper, the second approach was chosen in order to investigate the Cervaro watershed in Southern Italy.Firstly, the study area is described, then the contribution of remote sensing in generating Land use/Land cover (LULC) maps from multi-date LANDSAT data is presented.The methodology adopted is consequently discussed before illustrating and discussing the results and presenting the authors' conclusions.

The Study Area
The Cervaro river basin crosses the provinces of Avellino and Foggia in Campania and Apulia, respectively, two regions in Southern Italy (Figure 2).Like all the rivers crossing the Tavoliere delle Puglie plain, this basin stretches from south-west to north-east, being sourced in the mountains of Subappennino Dauno and flowing into the Adriatic Sea, between Candelaro and Carapelle outfalls.LANDSAT data is presented.The methodology adopted is consequently discussed before illustrating and discussing the results and presenting the authors' conclusions.

The Study Area
The Cervaro river basin crosses the provinces of Avellino and Foggia in Campania and Apulia, respectively, two regions in Southern Italy (Figure 2).Like all the rivers crossing the Tavoliere delle Puglie plain, this basin stretches from south-west to north-east, being sourced in the mountains of Subappennino Dauno and flowing into the Adriatic Sea, between Candelaro and Carapelle outfalls.The Cervaro river basin is one of largest in the Apulia region, with an extension of about 775 km 2 .It crosses four towns located in the province of Avellino, and fifteen towns located in the province of Foggia.Maximum and average height of the basin is 1100 and 350 m, respectively, with an average slope of 11%.The river stretches for about 110 km.
The basin is also important from a naturalistic viewpoint, as it includes a SPZ area (Special Protection Zone), a SIC area (Site of Community Interest), an Important Bird Area, and a Regional Natural Park.The Cervaro basin has a temperate maritime climate (or sub-tropical with arid summer).The average annual temperature is between 15 °C and 17 °C.The rainfall pattern is variable; in the upstream portion of the basin, the average annual precipitation is about 1000 mm, while in the The Cervaro river basin is one of largest in the Apulia region, with an extension of about 775 km 2 .It crosses four towns located in the province of Avellino, and fifteen towns located in the province of Foggia.Maximum and average height of the basin is 1100 and 350 m, respectively, with an average slope of 11%.The river stretches for about 110 km.
The basin is also important from a naturalistic viewpoint, as it includes a SPZ area (Special Protection Zone), a SIC area (Site of Community Interest), an Important Bird Area, and a Regional Natural Park.The Cervaro basin has a temperate maritime climate (or sub-tropical with arid summer).The average annual temperature is between 15 • C and 17 • C. The rainfall pattern is variable; in the upstream portion of the basin, the average annual precipitation is about 1000 mm, while in the downstream portion, it is about 400 mm.Having been the subject of previous studies, in the Cervaro basin there are seven rain gauges and a hydrometric station that is no longer functioning.
Between 1930 and 1980, the Cervaro basin was affected by many strong meteorological events.All of these events were dangerous and caused much damage to people, agricultural activities (which are a vital economic resource for the Apulia region [18,20]), and road systems-especially in the lower portion of the basin [21].Particularly, two severe floods occurred within the basin in 1956 and 1961, respectively, with peak discharge of 422 and 465 m 3 /s.

LANDSAT Data Processing
With the current availability of global spatial datasets, a common way to map land cover changes is through satellite image classifications.Different images can be processed separately, allowing variability at seasonal and directional scales to be accounted for [22][23][24][25][26].Many methods for mapping large areas with a high level of accuracy have been studied [1,27], and one of the most widespread is the Artificial Neural Network (ANN) approach, which can simulate the decision making processes of the human brain.ANNs have been used for land cover studies mainly based on raw band data, obtaining higher accuracies of classification than standard classification techniques [28,29].However, some drawbacks need to be improved in the method: some classes (e.g., vegetation) are difficult to process due to the high level of seasonal variation, and data are still lacking in many studies of historical ground reference.Spectral indices are suitable for the identification of specific features on the ground, such as water, vegetation, or imperviousness, therefore offering more information than raw band data and potentially resulting in higher classification accuracies [30].Patel et al. [31] demonstrated that-when used with spectral indices-ANNs contribute to reaching higher accuracies.
For the purposes of this study, four cloud-free LANDSAT-Thematic Mapper (TM) scenes were selected (Table 1) from the Earth Resources Observation Systems (EROS) Data Center (EDC)to perform the Land Cover (LC) classification and, consequently, the extraction of the runoff Curve Number (CN).LANDSAT-TM scenes are provided to users as Level 1 terrain corrected (L1T) data and distributed under the GEOTIFF format with the Universal Transverse Mercator (UTM) projection and the World Geodetic System 84 (WGS84) datum.A geometric correction phase is always necessary in multisensor analysis, or when the data are not previously corrected [32,33].Since Landsat-TM are provided as L1T data and come from the same sensor, a further orthorectification process was not applied to avoid a resampling with a consequent alteration of pixels Digital Numbers values.
Pre-processing of LANDSAT-TM data was aimed at enhancing the classification performance.For this purpose, nine synthetic bands were extracted through principal component analysis (PCA), tasseled cap (TC), brightness temperature (BT), and vegetation indices (Figure 3).Vegetation indices were used to enhance the detection of vegetated areas: leaf area index (LAI) and modified soil adjusted vegetation index (MSAVI).LAI is commonly used in agronomic and environmental modeling [36][37][38][39][40][41].In this study, LAI was evaluated through a semi-empirical model derived from the Lambert-Beer law (1): LANDSAT-TM digital numbers were firstly converted to Top of Atmosphere (TOA) reflectance and corrected for sun elevation angle by using gains, offsets, and brightness temperature coefficients provided in each scene metadata and related standard procedures suggested by the USGS [34].The dark object subtraction (DOS) [35] procedure was implemented to reduce atmospheric artifacts.Vegetation indices, PCA, and TC bands were computed from the DOS output of each scene.
Vegetation indices were used to enhance the detection of vegetated areas: leaf area index (LAI) and modified soil adjusted vegetation index (MSAVI).LAI is commonly used in agronomic and environmental modeling [36][37][38][39][40][41].In this study, LAI was evaluated through a semi-empirical model derived from the Lambert-Beer law (1): Where NDVI max and NDVI back are, respectively, an upper and a lower bound of the Normalized Difference Vegetation Index (NDVI), and k is an extinction coefficient.MSAVI is one of the vegetation indices that makes use of the soil line [42,43], which is the regression line of the distribution of soil pixels within the red-near-infrared (NIR)plane of the feature space with slope "a" (2)(3)(4).
The soil slope a = 1.47 was evaluated from a group of manually selected pixels (common for each considered scene) and chosen as representative of bare soils.For each subset, the WDVI (weighted difference vegetation index), the M parameter, and the MSAVI were, respectively, evaluated.
For each scene, the same training pixel locations were used.When ground reference was not available (i.e., lacking ancillary data), selected training pixels were considered as belonging to the same class of scenes with availability of ancillary data.The use of synthetic data/channels derived from the original ones has already proven its effectiveness in recent studies [44,45].In this paper, the positive effect of synthetic channels was proven by evaluating the Jeffries-Matusita Distance (JMD) [46] for each i-th and j-th class.Figure 4 shows an example of the great enhancement of spectral separability achieved by using the aforementioned combination of synthetic bands.Where NDVI max and NDVI back are, respectively, an upper and a lower bound of the Normalized Difference Vegetation Index (NDVI), and  is an extinction coefficient.MSAVI is one of the vegetation indices that makes use of the soil line [42,43], which is the regression line of the distribution of soil pixels within the red-near-infrared (NIR)plane of the feature space with slope "a" (2)(3)(4).
The soil slope a = 1.47 was evaluated from a group of manually selected pixels (common for each considered scene) and chosen as representative of bare soils.For each subset, the WDVI (weighted difference vegetation index), the M parameter, and the MSAVI were, respectively, evaluated.
For each scene, the same training pixel locations were used.When ground reference was not available (i.e., lacking ancillary data), selected training pixels were considered as belonging to the same class of scenes with availability of ancillary data.The use of synthetic data/channels derived from the original ones has already proven its effectiveness in recent studies [44,45].In this paper, the positive effect of synthetic channels was proven by evaluating the Jeffries-Matusita Distance (JMD) [46] for each i-th and j-th class.Figure 4 shows an example of the great enhancement of spectral separability achieved by using the aforementioned combination of synthetic bands.The enhancement of the classification process was also achieved by substituting the original channels with their linear transformations: tasseled cap (TC) transformation and principal component (PC).The TC transformation [47] is one of the techniques designed to enhance the informative content of multiband imagery.TC synthetic channels highlight the evolution of the vegetation spectral response in a multitemporal study.PC analysis is aimed at the reduction of the dimensionality of multiband data through the computation of uncorrelated synthetic bands [48].Only the first three PCs were used for each scene, since within them more than the 99% of the total variance of each scene is enclosed.
Training and ground reference pixels were first collected through random and stratified sampling from land use cartography and then validated with Google Earth.The totality of collected pixels was divided in two groups to verify classification accuracy without using Artificial Neural The enhancement of the classification process was also achieved by substituting the original channels with their linear transformations: tasseled cap (TC) transformation and principal component (PC).The TC transformation [47] is one of the techniques designed to enhance the informative content of multiband imagery.TC synthetic channels highlight the evolution of the vegetation spectral response in a multitemporal study.PC analysis is aimed at the reduction of the dimensionality of multiband data through the computation of uncorrelated synthetic bands [48].Only the first three PCs were used for each scene, since within them more than the 99% of the total variance of each scene is enclosed.
Training and ground reference pixels were first collected through random and stratified sampling from land use cartography and then validated with Google Earth.The totality of collected pixels was divided in two groups to verify classification accuracy without using Artificial Neural Network Classification (ANNC) training pixels.The tests carried out showed that a simple stratified sampling scheme was not adequate to achieve an elevated spectral separability.Binary masks were used to select problematic pixels from the original channels and rearrange them.The general scheme followed was the creation of a small region of interest (ROI) related to a specific category of pixels, the extraction of band statistics for ROI pixels, and the creation of the i-th mask (5).
With b i the TOA reflectance value for the i-th channel, b i,max and b i,min , respectively, the maximum and the minimum for the i-th channel of the selected ROI and n number of channels considered.This process ended with thirteen output classes (Table 2).An ANNC allows the analysis of complex relationships among variables without making a priori assumptions about data (i.e., normality or linearity) [49,50].In this study, the ANNC was implemented through the Matlab environment by Mathworks [51].Each neuron implements a logsig transfer function ( 6) and a bias.
With n ∈ (−∞; +∞) and logsig(n) ∈ (0; 1).The log-sigmoid transfer function is commonly used in multilayer networks that are trained using the family of backpropagation algorithms because it is differentiable [52].The best configuration found was a four layer ANNC with one input layer with nine neurons (one for each synthetic band), two hidden layers with respectively twelve and thirteen neurons, and one output layer with thirteen neurons (one for each output class).In order to evaluate ANNCs output accuracy, four confusion matrices were computed using 150 ground truth pixels for each class (Table 2).Final map results are shown in Figure 5.

Curve Number Extraction
The runoff curve number is a method used in hydrology as a simple watershed model, and as the runoff estimating component in more complex watershed models [53,54].The wide use of this method is not only due to its applicability in different environments, but also to its simplicity and its dependence from soil types, land use, and antecedent moisture conditions [55].

Curve Number Extraction
The runoff curve number is a method used in hydrology as a simple watershed model, and as the runoff estimating component in more complex watershed models [53,54].The wide use of this method is not only due to its applicability in different environments, but also to its simplicity and its dependence from soil types, land use, and antecedent moisture conditions [55].
In order to produce a CN extraction, the classes extracted from LANDSAT-TM imagery were merged to fit the Soil Conservation Service (SCS) hydrological cover types.The preformed association (Table 3) takes into account the great land cover heterogeneity that characterizes the study area and the fact that some ANNC output classes show a similar hydrological By using the Apulian Lithological Map [56], hydrological soil groups were set in order to extract the correct CN for each cell.

Hydrologic Modeling
The basic idea of the study was to estimate the impact of land use change on runoff within a river basin, and there were four scenarios identified for this analysis : 1984, 2003, 2009, and 2011.With the aim of estimating floods for these four scenarios, two types of data are needed: digital elevation data and land use data.As the study area relates to two regions (Apulia and Campania), two different data sets were merged: a Digital Elevation Model (DEM) of 8 m cell size of Apulia and a CTR (Regional Technical Map) of Campania.
In order to evaluate land use changes and their effects on runoff, curve number maps of each scenario were produced.Indeed, the CN method is recommended when it is necessary to accurately take into account the contribution due to land use changes over time [57].CN is a parameter that takes into account the lithology of land and the characteristics of land use and agricultural practices.The CN variation range is from 0 to 100; hence, the basin was divided into 14 sub-basins to estimate floods with a semi-distributed model.For each sub-basin, the parameters of height, slope, surface, length of main shaft, and CN values (CNI, CNII, and CNIII) were estimated [57,58].
Precipitation was calculated in each sub-basin using the VAPI Puglia method (i.e. an Italian research project, based on the use of Two-Component Extreme Value distribution probabilistic model) [59], which is useful for the estimation of rainfall within a basin when rainfall data are poor.The VAPI procedure is a regional methodology based on two-component extreme value (TCEV) distribution.It assumes that individual rainfall can be expressed as a mixture of two exponential components: (a) the first describes ordinary events (more frequent and less severe on average); (b) the second describes extraordinary extremes (more severe and less frequent).
The VAPI Puglia method [59] is a methodology performed on the Apulian territory and recommended by the competent authorities in the territory under study.
In order to estimate flow depth and the velocity of flood events, a two-dimensional numerical model called FLO-2D [60] was adopted.Numerical simulations were performed using an empirical approach, evaluating the peak flow rate and the triangular hydrograph for each sub-basin of the four scenarios of land use only for the CNII value, with a return time equal to 200 years.
The choice to use a return period equal to 200 years was due to the concept of "hydraulic safety" regulated by PAI Puglia (i.e., the hydrogeological setting plan) [61].This approach allows an understanding of how the change of land use affects the "hydraulic safety" in the study area.
Due to its large extension, the basin was divided into four domains, and flow depth and velocity was estimated for each of them.The outflow of each domain was used as inflow for the related domain.
Numerical simulations were carried out using the two-dimensional FLO-2D Model, a volume conservation model that numerically routes a flood hydrograph while predicting the area of inundation and flood wave attenuation over a system of square grid elements.The model considers river channel flow as one-dimensional, and simulates river overbank flow and floodplain flow as two-dimensional.In this case, the DEM adopted to simulate flow on a floodplain topography was provided by the Regional Authorities of Apulia.The computational domain was divided in four sub-domains due to the extension of the river and the complexity of the floodplain (Figure 6).The FLO-2D features used in this model also included bridges (in the FLO-2D software they are called "hydraulic structures"), roughness shape, and flow obstructions.Flood routing in two dimensions was accomplished through numerical integration of the motion equation (full dynamic wave momentum equation) and the conservation equation of fluid volume.
The simulation of the flooding which develops in two dimensions is obtained through the numerical integration of the constitutive equations ( 7) and (8).
The constitutive equations are [60]: Conservation equation of fluid volume: Motion equation (full dynamic wave momentum equation): Where h is the flow depth and V is the depth-averaged velocity in one of the eight flow directions x.The excess rainfall intensity (i) may be nonzero on the flow surface.The friction slope component S fx is based on Manning's equation.The other terms include the bed slope S ox pressure gradient and convective and local acceleration terms.This equation represents the one-dimensional depth averaged channel flow.
A two-dimensional model was used to provide a realistic map of flooding areas, which showed that the flow is unconfined during critical events in the Cervaro basin.
For the purposes of this study, a change of land use was simulated by changing the input in each domain because the inputs are the hydrographs, which were estimated as the output of the SCS method.
A DEM and the inflow input data given by the sum of the output of the previous domain were used for each domain, as well as the inflow of each sub-basin included in the domain (Figure 6).The simulation of the flooding which develops in two dimensions is obtained through the numerical integration of the constitutive equations ( 7) and (8).
The constitutive equations are [60]: Conservation equation of fluid volume:

Results
The CNII maps for the four scenarios extracted with the above-mentioned methodology are shown in Figure 7, Figure 8 shows the CNII differences over the whole basin experienced between two consecutive satellite data.Particularly, the analysis of CNII differences over time should take into account that the CNII inherits the accuracy of the performed classifications.The temporal CN trend analysis and CNII values can be observed in Figure 8, showing that the CN increase occurred mostly with the same spatial distribution.Overall, there was only a low average increase in CN value during the 28 years; locally, however, a non-negligible increase can be observed, especially for Domain 3 (which shows a built-up area increase), and Domain 1, in the delta area of the Cervaro basin.The observed variability can also be explained considering that the input satellite data are not n-days synthesis data (e.g., Level-3 MODIS, Spot VGT, and PROBA-V product), and that the performed land cover classifications have to be considered dependent on the specific acquisition time.Particularly, in the considered time period (summer), the basin landscape is the object of several agricultural practices (e.g., harvesting and pruning) that influence the mean spectral response at the LANDSAT-TM pixel level.This is clearly shown in Figure 7, in which most of the change is observed with reference to land covers potentially affected by seasonal and agricultural human activities.The results of the application of the two-dimensional numerical models were also analyzed for each spatial domain investigated in relation to the different land use scenarios chosen.Flooding areas for the four scenarios are shown in Figure 9, superimposing their results, with the oldest on the top.Areas characterized by a continued increase in hydraulic hazard is reported in detailed images, observing that there is a good correspondence between these areas and that with a CN increase.
In order to estimate land use change impact within the flooding areas, a comparative analysis The temporal CN trend analysis and CNII values can be observed in Figure 8, showing that the CN increase occurred mostly with the same spatial distribution.Overall, there was only a low average increase in CN value during the 28 years; locally, however, a non-negligible increase can be observed, especially for Domain 3 (which shows a built-up area increase), and Domain 1, in the delta area of the Cervaro basin.The observed variability can also be explained considering that the input satellite data are not n-days synthesis data (e.g., Level-3 MODIS, Spot VGT, and PROBA-V product), and that the performed land cover classifications have to be considered dependent on the specific acquisition time.Particularly, in the considered time period (summer), the basin landscape is the object of several agricultural practices (e.g., harvesting and pruning) that influence the mean spectral response at the LANDSAT-TM pixel level.This is clearly shown in Figure 7, in which most of the change is observed with reference to land covers potentially affected by seasonal and agricultural human activities.
The results of the application of the two-dimensional numerical models were also analyzed for each spatial domain investigated in relation to the different land use scenarios chosen.Flooding areas for the four scenarios are shown in Figure 9, superimposing their results, with the oldest on the top.Areas characterized by a continued increase in hydraulic hazard is reported in detailed images, observing that there is a good correspondence between these areas and that with a CN increase.In order to estimate land use change impact within the flooding areas, a comparative analysis between the different scenarios in each domain was also carried out.Total flooding area and volume, average flow depth and velocity for each domain were compared, showing a clear upward trend for each of these parameters (Figure 10).The remaining domains did not seem influenced by land use changes, with hydraulic parameters being almost invariable.
As further confirmation of the achieved results, it is possible to consider the areal reduction factor (ARF).The ARF is a key quantity in the design against hydrologic extremes.For a basin of area "A" and a rain duration "d" and a return period "T", ARF(A, d, T) is the ratio between the average rainfall intensity in A and d with return period T and the average rainfall intensity at a point for the same d and T. Empirical ARF charts often display scaling behavior [22,62].
In the "VAPI Puglia", the ARF can be calculated with following equation ( 9 Where the basin area A is expressed in kilometers, the duration of the rain d in hours, while the values of the three coefficients c1, c2, and c3 vary according to the compartment [59].Like with CNII values, the domains showing an increase in the extent of flooding area were Domains 1 and 3; more specifically, Domain 1 showed a significant upwards trend only in terms of flow velocity, while the trend of other variables was not clear.Domain 3 instead showed a continuous increase only in the extent of the flooding area.
The remaining domains did not seem influenced by land use changes, with hydraulic parameters being almost invariable.
As further confirmation of the achieved results, it is possible to consider the areal reduction factor (ARF).The ARF is a key quantity in the design against hydrologic extremes.For a basin of area "A" and a rain duration "d" and a return period "T", ARF(A, d, T) is the ratio between the average rainfall intensity in A and d with return period T and the average rainfall intensity at a point for the same d and T. Empirical ARF charts often display scaling behavior [22,62].
In the "VAPI Puglia", the ARF can be calculated with following equation ( 9): ARF = 1 − 1 − e −c 1 A e −c 2 d c 3 (9) Where the basin area A is expressed in kilometers, the duration of the rain d in hours, while the values of the three coefficients c 1 , c 2 , and c 3 vary according to the compartment [59].
For the total extension of the considered basin, ARF is less than one.However, the CN effect was emphasized by dividing the whole basin into several subcatchments.

Conclusions
The estimation of flooding areas for a basin depends on many factors, including hydrological parameters, morphological characteristics, land use, and adopted numerical model.Each basin is unique, and for this reason it must be analyzed in detail, taking into account all of the above-mentioned aspects.
Starting from this awareness, and by taking an integrated and multi-disciplinary approach, this paper has analyzed flooding areas for the Cervaro basin in Apulia that showed non-negligible land use changes for the 30 years under study.
In particular, the integration of scientific information and the analysis of land use changes using LANDSAT-TM scenes from 1984 to 2011 were investigated, with a two-dimensional hydraulic numerical model capable of assessing flooding area extension for each land use scene.
The analysis of the results showed a good correlation between flooding areas and land use changes, especially a slight increase in flooding areas related to a rise in the waterproofing due to urbanization growth.In fact, the areas with the greater CN value coincide with areas where urbanization or abandonment of agricultural practices has led to an increase of impervious areas.This increase justifies the increments in terms of flooding areas.
This correlation, however (as shown in Figure 1) can be considered negligible in a basin of large dimensions, such as the one analyzed in this study.
These results confirm what was previously asserted by Bloschl et al. [4]: the impact of land use on the hydrological response is a function of watershed scale.Abrupt changes in flooding can be particularly severe, especially at small scales.
Investigating how land use was changed, Gholami et al. [15] found that urban and road network development intensified runoff vegetation and flooding hazard within an Iranian forest area during the last 40 years.Results indicated an increased influence of urban development on the volume and peak discharge in sub-basin areas, which had considerable impervious surface development.

Figure 1 .
Figure 1.Hypothesized impact of land use and climate variability on hydrological response as a function of scale (Source: Bloschl et al. [4]).

Figure 1 .
Figure 1.Hypothesized impact of land use and climate variability on hydrological response as a function of scale (Source: Bloschl et al. [4]).

Figure 3 .
Figure 3. Flowchart of the methodology performed on LANDSAT data.CN: curve number; LULC: Land Use/Land Cover.

Figure 3 .
Figure 3. Flowchart of the methodology performed on LANDSAT data.CN: curve number; LULC: Land Use/Land Cover.

Figure 4 .
Figure 4. Jeffries-Matusita Distance (JMD) for different combinations of bands related to the 5th class for the 2011 scene.LAI: leaf area index; MSAVI: modified soil adjusted vegetation index; PC: principal components; TC: tasseled cap; NDVI: Normalized Difference Vegetation Index.

Figure 4 .
Figure 4. Jeffries-Matusita Distance (JMD) for different combinations of bands related to the 5th class for the 2011 scene.LAI: leaf area index; MSAVI: modified soil adjusted vegetation index; PC: principal components; TC: tasseled cap; NDVI: Normalized Difference Vegetation Index.

Figure 6 .
Figure 6.The four sub-domains along the Cervaro river.

Figure 6 .
Figure 6.The four sub-domains along the Cervaro river.

Figure 10 .
Figure 10.Total flooding area and volume, average flow depth and velocity for each domain.

Figure 10 .
Figure 10.Total flooding area and volume, average flow depth and velocity for each domain.

Table 2 .
Overall accuracies achieved for each classification output.

Table 3 .
ANNC output classes and their associated Soil Conservation Service (SCS) hydrological cover classes.