Hydrological Response to ~30 Years of Agricultural Surface Water Management

Amongst human practices, agricultural surface-water management systems represent some of the largest integrated engineering works that shaped floodplains during history, directly or indirectly affecting the landscape. As a result of changes in agricultural practices and land use, many drainage networks have changed producing a greater exposure to flooding with a broad range of impacts on society, also because of climate inputs coupling with the human drivers. This research focuses on three main questions: which kind of land use changes related to the agricultural practices have been observed in the most recent years (~30 years)? How does the influence on the watershed response to land use and land cover changes depend on the rainfall event characteristics and soil conditions, and what is their related significance? The investigation presented in this work includes modelling the water infiltration due to the soil properties and analysing the distributed water storage offered by the agricultural drainage system in a study area in Veneto (northeastern Italy). The results show that economic changes control the development of agro-industrial landscapes, with effects on the hydrological response. Key elements that can enhance or reduce differences are the antecedent soil conditions and the climate characteristics. Criticalities should be expected for intense and irregular rainfall events, and for events that recurrently happen. Agricultural areas might be perceived to be of low priority when it comes to public funding of flood protection, compared to the priority given to urban ones. These outcomes highlight the importance of understanding how agricultural practices can be the driver of or can be used to avoid, or at least mitigate, flooding. The proposed methods can be valuable tools in evaluating the costs and benefits of the management of water in agriculture to inform better policy decision-making.


Introduction
Raising interest in the importance of human drivers to understand the past and current development of our planet is of great importance [1,2].Both the landscape and the river systems in the urbanised floodplains of Europe, and elsewhere, have undergone major changes in the past, and these environmental changes have altered the flooding conditions [3,4].Degraded processes from incorrect land managements have also been underlined in many countries of the world [5,6], for different land uses [7][8][9], and at different scales [10][11][12].However, due to the complexity of the elements involved, the magnitude of the impacts is still highly uncertain [3].Amongst human practices, agricultural surface water management systems represent some of the largest integrated engineering works undertaken by humans that shaped floodplains during history, directly or indirectly affecting the landscape [13].On the one hand, in industrialised floodplains, surface water management is connected to issues arising due to changes in agricultural practices, urbanisation patterns [14,15], and emerging bank failure mechanisms [16].On the other, in developing countries, issues at stake for water management are the need to increase food production together with the increase of water shortages [17].Both examples present a broad range of impacts on society, in particular for the protection of flood prone areas, because of climate inputs coupling with human drivers [4,15,18].Several hundred years ago, small farming communities were already having a significant impact on floodplains [19,20] increasing runoff and, thus, flooding.Recently, however, the importance of cooperation of communities and landowners for the successful implementation of natural flood management projects has been underlined [21][22][23].Specifically, in the landscape of northern Italy, where industrial agro-ecosystems overlap with residential areas [24], agricultural water management is of great importance.Italy's vulnerability to floods is amplified by weak enforcement of building restrictions and low compliance with sound floodplain management principles, inevitably leading to higher flood losses [25].A common approach to flood risk management, thus, should involve an assessment of the characteristics of, and linkages between sources (e.g., precipitation and climate), preferential storage (e.g., surface water storage), and soil response to the input (areas where ponding can happen).Many countries already have some form of flood risk mapping, but these leave out flooding from malfunctioning drains, groundwater flooding, ponding, and flooding from small watercourses and culverts [26].
This research focuses on the changes of surface land drainage caused by the transformations of agricultural practices, and on the related criticalities.Agricultural land use potentially has an important role in flood risk management.Runoff and subsurface drainages from farmland act as a pathway, causing flooding in downstream receptor areas.Surface drainage can offer water storage to laminate flood waves coming from upstream.Furthermore, in relatively flat areas, rainfall is stored in the ground, depending on the soil properties, and when the ground capacity is exceeded the water leaves the system causing flooding.The analysis presented in this paper focuses on two main questions: which kind of agricultural land-use changes have been observed in the recent years (~30 years)?How does the influence of land use and land cover changes on a proxy of the watershed response depend on rainfall event characteristics and soil conditions, and what is their related significance?Research, in fact, is still needed to understand if rainfall characteristics and climate have a synergistic effect, if their interaction matters, or to understand what element has the greatest influence.

Study Area
The research focuses on the floodplain area of the municipality of Galzignano Terme (Figure 1) in the Province of Padua (in Veneto, North-Eastern Italy).The study area is located about 50 km south-west of Venice and 15 kilometres south-west of Padua.The whole municipality is an area in continuous urbanistic evolution, with a high level of hydrologic vulnerability due to flooding and ponding.Overall, the municipality is composed of a hilly side with steep slopes (>30%), and a reclaimed floodplain area with elevations ranging between 3-5 m a.s.l.The urban development in the municipality covers only a limited extent, and it is mainly composed of residential areas and hotels.The main agricultural activities include vineyards, olive orchards, crop cultivations and pastures, with more intensive crop and vines production near the urbanised areas.
The available topographic information for the study area is a LiDAR Digital Terrain Model (DTM) surveyed in 2008-2010 (horizontal accuracy of about ±0.3 m and vertical accuracy of ±0.15 m) (Figure 1a), together with aerial photographs of 1983 (Figure 1b) and orthophotos of 2006 (Figure 1c).The drainage system of the area has changed significantly over the years, but has remained mainly composed of a minor agricultural drainage network, and three major channels built between 1919 and 1930, with pumping stations controlling the water flow.The floodplain portion is highly involved in hydraulic criticalities, due to different forcing.In the higher part of the floodplain, issues arise because of change in slopes between the hills and the plain.For this area, the extent of the channel network decreased over the years, due to urbanisation.As a consequence, the storage capacity diminished, causing problems of occlusions when the streams from the hills flow into the artificial system.The plain area, instead, presents issues due to water ponding and lack of maintenance of the artificial drainage system [27].
In addition to the topographic information, local authorities provide a database of some characteristic of the soil (Table 1), which are homogeneous over the considered extent [28].The given values come from pedotransfer functions (PTFs), established to translate some measured soil matrix properties such as bulk density, organic matter and soil texture into saturation [29].The PTFs considered come from different sources [30].When texture values are within certain ranges (clay 0%-49% and sand 0%-50%) the PTFs used are those tested on data collected on benchmark soils from experimental sites in the Pianura Padano-Veneta, Northern Italy [31].Otherwise, hydraulic conductivity derives from PTFs calibrated on clay, sand and bulk density according to [32]; while the soil water retention stems from the PTFs defined by [33], based on field surveyed values of percentage of clay, sand and organic carbon.
Table 1.Main soil characteristics that are available for the study area [28].

Parameter Value
Effective hydraulic conductivity Ks 43.79 mm•h −1 Matric pressure at the wetting front Ψf 130 mm Saturated moisture content θs 0.63

Changes in the Agricultural System
In Italy, towards the end of the 20th century and at the start of the new millennium, agriculture followed a development strongly characterised by local peculiarities, with irrigation practices as a The drainage system of the area has changed significantly over the years, but has remained mainly composed of a minor agricultural drainage network, and three major channels built between 1919 and 1930, with pumping stations controlling the water flow.The floodplain portion is highly involved in hydraulic criticalities, due to different forcing.In the higher part of the floodplain, issues arise because of change in slopes between the hills and the plain.For this area, the extent of the channel network decreased over the years, due to urbanisation.As a consequence, the storage capacity diminished, causing problems of occlusions when the streams from the hills flow into the artificial system.The plain area, instead, presents issues due to water ponding and lack of maintenance of the artificial drainage system [27].
In addition to the topographic information, local authorities provide a database of some characteristic of the soil (Table 1), which are homogeneous over the considered extent [28].The given values come from pedotransfer functions (PTFs), established to translate some measured soil matrix properties such as bulk density, organic matter and soil texture into saturation [29].The PTFs considered come from different sources [30].When texture values are within certain ranges (clay 0%-49% and sand 0%-50%) the PTFs used are those tested on data collected on benchmark soils from experimental sites in the Pianura Padano-Veneta, Northern Italy [31].Otherwise, hydraulic conductivity derives from PTFs calibrated on clay, sand and bulk density according to [32]; while the soil water retention stems from the PTFs defined by [33], based on field surveyed values of percentage of clay, sand and organic carbon.
Table 1.Main soil characteristics that are available for the study area [28].

Parameter Value
Effective hydraulic conductivity K s 43.79 mm•h −1 Matric pressure at the wetting front Ψ f 130 mm Saturated moisture content θ s 0.63

Changes in the Agricultural System
In Italy, towards the end of the 20th century and at the start of the new millennium, agriculture followed a development strongly characterised by local peculiarities, with irrigation practices as a Land 2017, 6, 3 4 of 24 dominant feature [34].To describe these changes in the agricultural system and the connected changes in the water irrigation management at the local scale, the research considers different indicators taken from statistical reports [35][36][37][38][39], broadly divided in irrigation and farm indicators.
Among the considered irrigation statistics, the analysis focuses on the irrigated surface and the type of irrigation.To this latter point, non-structural irrigation refers to the use of water taken by the farmers directly from the reclamation network when needed, while organised one refers to specific amounts of water, regulated times, and specific irrigation volumes defined by official local regulations [40].
The farm system is characterised based on agricultural holdings and information about their number and type, including characteristics as management (entrepreneurial, exclusive properties, family-run business), extent, Utilized Agricultural Area (UAA), economic dimension (UDE, or gross income), and production crops.
The type of farm management is critically connected to the irrigation system: changes at the farm level imply shifts in the kind of irrigation applied (non-structural or organised).However, while the organised irrigation parameters are known, there is no accurate information about the non-structural ones at a local scale [40].Changes in the farm extent also determine variations in the water system: the use of large farm equipment demands an increasing in field size, with resulting elimination of buffer zones (areas of land, lying next to a waterway, kept in permanent vegetation), ditches and hedgerows [41].UUA and production are as well indicators of surface water management.Crops are the cultivation with the highest irrigation requirements, while the UUA extent implies information of the extent of surface water irrigation.Further, at the regional level, there is a correlation between changes in the UDE of the farms, and changes in the irrigation system [37].Activities with less than 4800 € of gross income have an irrigated surface equal to ~30% of their UAA, while the irrigated surface increases with the growth of the gross income, reaching about 60% for those activities having more than 48k € of UDE [40].

Network Extraction for 1983 and 2010
The network storage capacity has been evaluated for the reference years using two methodologies.For the current settings (year 2010), we considered a procedure already applied in [15,42].The proposed method can be easily implemented because the characteristics of the artificial channels in this context are constant over the surface, and they are strictly related to the trenchers used to build the channels.The overall aim of the approach is to estimate an average channel width and to estimate the water storage per areal units over some user-specified areas.Starting from a high-resolution LiDAR DTM, the procedure allows the automatic identification of the drainage network, and it allows to estimate some indexes of network width and length.Subsequently, associating some parameters describing the cross-section of the channels for specific channel width intervals, it allows to estimate the water storage within the channels [42].The authors refer to [15,42] for a detailed description of the procedure.For this work, the study site was divided into subsections of 50 m × 50 m.This classification is an assumption, and in this case, it is small enough to capture the geometry of the network, and wide enough to capture permanent variations in soil sealing (urban expansion).
For the past conformation of the network, the aerial photographs of 1983 (Regione del Veneto-LR n. 28/76 Formazione della Carta Tecnica Regionale) were considered to identify the channels manually.For each subarea, we assumed an average channel cross-section equal to the mean value estimated for the year 2010 [15].Water storage volumes were then estimated starting from the network length and cross-section area for each subsection.

Updated Network Saturation Index (uNSI)
For low plain areas with null slopes, the effects of network changes can be assessed using an indicator named Network Saturation Index (NSI) [15] that provide a measure of how long it takes for a designed rainfall to saturate the available storage volume.Given a designed rainfall duration and rainfall amount, and assuming that the amount of rainfall is homogeneous over the considered surface, it is possible to compute at every time step the percentage of storage volume that is filled by the rainfall.The NSI is then the first time step at which the available storage volume is 100% reached.
The original NSI approach considered a situation where there was no infiltration, and the method provided a single value of NSI over an entire area.The present research updates the index in two ways: 1.
it includes infiltration according to the Green and Ampt [43] approach; 2.
the NSI is evaluated as localised value over each user-defined subarea.
The first improvement allows considering how the hydraulic characteristics of the soil matrix influence the runoff production.The second improvement allows giving a distributed analysis of the index, thus identifying areas that can be more or less critical, or regions where changes in the land use might have improved, or decreased, the hydraulic efficiency of the landscape.
The updated NSI (uNSI) provides a local measure of how long it takes for a designed rainfall to produce a specific runoff (depending on the soil characteristics and antecedent moisture condition, see following paragraphs) that saturates the available storage volume.
The uNSI evaluation starts with synthetic rather than actual rainfall events, and it considers specific Depth-Duration Frequency curves (DDF) [15].Entering into detail on the construction of DDF curves is beyond the scope of this work, but the authors refer to [44] for a complete review of the available methods in the literature.The basic steps applied for this research follow [15] and are 1.
Plot depth versus duration for various frequencies and derive the DDF curves in the form where h is the rainfall amount in mm, t is the time in hours, and a and n are empirically derived coefficients.For the simulations, the analysis considers different return periods (from 2 up to 200 years).Originally, the NSI [15] was evaluated to create a storm producing, for the shortest return period, a volume of water about ten times larger than the total volume network storage capacity.However, the selection of the rainfall duration is an operational choice [15].For the proposed work, the uNSI is computed locally.Thus, we decided to consider a duration producing an amount of rainfall, for the shortest return period, able to fill the 90th percentile of the available volumes (~1 h for this study site).This choice is justified by the fact that the location presents ponds and retention basins built in 2010, thus choosing a rainfall duration according to [15] would produce a rainfall too high to be realistic, for the shorter return period.
Given the input rainfall, the Green-Ampt approach [43] allows computing the infiltration and, thus, it allows to evaluate runoff.The model yields through simple integration to an exact solution relating the infiltration rate, cumulative infiltration, and ponding time [51] under ponded condition into a deep homogeneous soil with uniform initial moisture content [52], giving results that match empirical observations.Conceptually, for the evaluation of the uNSI, the storage volume offered by the drainage network will start filling with water only after the soil is completely saturated.Thus, an important parameter is the ponding time (Equation (2)), or the time at which infiltration stops and runoff starts.For a rain storm of a specific duration and intensity (P), tp is given by where s K = effective hydraulic conductivity; f ψ = matric pressure at the wetting front; s θ = saturated moisture content; and i θ = initial (antecedent) moisture content.
Given the input rainfall, the Green-Ampt approach [43] allows computing the infiltration and, thus, it allows to evaluate runoff.The model yields through simple integration to an exact solution relating the infiltration rate, cumulative infiltration, and ponding time [51] under ponded condition into a deep homogeneous soil with uniform initial moisture content [52], giving results that match empirical observations.Conceptually, for the evaluation of the uNSI, the storage volume offered by the drainage network will start filling with water only after the soil is completely saturated.Thus, an important parameter is the ponding time (Equation (2)), or the time at which infiltration stops and runoff starts.For a rain storm of a specific duration and intensity (P), t p is given by where K s = effective hydraulic conductivity; ψ f = matric pressure at the wetting front; θ s = saturated moisture content; and θ i = initial (antecedent) moisture content.The local authorities [28] provided the input required for the Green-Ampt model and to solve Equation (2) (Table 1).Hydraulic conductivity of the soil matrix dynamically responds to changes in the surrounding environment.Therefore, the infiltration parameters for the Green-Ampt equation should change for each storm event [53].To account for this, different saturation conditions (0%, 1%, 10%, 30%, 50%, 70%, 99%, 100%) were simulated, thus θ i = θ s × n, with n = 0.00, 0.01, 0.10, 0.30, 0.50, 0.70, 0.99, 1.00 (3) with θ s = saturated moisture content; and θ i = initial (antecedent) moisture content.
If the area presents urban surface or water bodies, a saturation of 100% is assumed, thus providing an effective rainfall that coincides with the input rainfall.
As a result, the actual runoff to be considered in the uNSI evaluation is the effective rainfall after ponding time (Figure 3).If the area presents urban surface or water bodies, a saturation of 100% is assumed, thus providing an effective rainfall that coincides with the input rainfall.
As a result, the actual runoff to be considered in the uNSI evaluation is the effective rainfall after ponding time (Figure 3)

Statistical Influence of Climate on Ponding Time and uNSI
Hyetographs are characterised to provide quantitative information on interstorm variability to define the effect of climate on the watershed response, regarding the ponding time and differences in uNSI.
For this characterization, we considered the rainfall bursts within the storm, where a burst is an abrupt change in rainfall rate [54].Four criteria were defined (Table 2), specifically: 1. the percentage of the total storm period that has occurred at the start of the heaviest burst (BrLoc%); 2. the number of bursts (Brn); 3. the shape of the hyetograph defined as its "triangularity" (TRN).To this point, the triangularity is the ratio between the hyetograph area and the area of a triangle having as base the hyetograph duration, and as height, the height of the main rainfall burst.Values close to 1 would represent highly triangular hyetographs (e.g., Sym02 to Sym04), with rainfall having regular intensity, whereas values close to 0 would represent hyetographs whose peaks is an abrupt change in intensity (e.g., Sym10 to Sym17); 4. the asymmetry of the storm volumes (asym), defined as the ratio between the amount of rainfall fell before the maximum burst and that after the burst.Values higher than 1 imply that the rainfall amount before the burst was greater than that after the burst; values close to 1 imply a

Statistical Influence of Climate on Ponding Time and uNSI
Hyetographs are characterised to provide quantitative information on interstorm variability to define the effect of climate on the watershed response, regarding the ponding time and differences in uNSI.
For this characterization, we considered the rainfall bursts within the storm, where a burst is an abrupt change in rainfall rate [54].Four criteria were defined (Table 2), specifically: 1.
the percentage of the total storm period that has occurred at the start of the heaviest burst (BrLoc % ); 2.
the number of bursts (Br n ); 3.
the shape of the hyetograph defined as its "triangularity" (TRN).To this point, the triangularity is the ratio between the hyetograph area and the area of a triangle having as base the hyetograph duration, and as height, the height of the main rainfall burst.Values close to 1 would represent highly triangular hyetographs (e.g., Sym02 to Sym04), with rainfall having regular intensity, whereas values close to 0 would represent hyetographs whose peaks is an abrupt change in intensity (e.g., Sym10 to Sym17); Land 2017, 6, 3 8 of 24 4. the asymmetry of the storm volumes (asym), defined as the ratio between the amount of rainfall fell before the maximum burst and that after the burst.Values higher than 1 imply that the rainfall amount before the burst was greater than that after the burst; values close to 1 imply a symmetry of the rainfall volumes; values smaller than 1 imply that the rainfall volume after the burst is higher than that before.These criteria highlight some climatic characteristics.Some studies showed that changes in climate could lead to heavy rainfalls beginning earlier, lasting longer, and being more continuous within a storm [55].Thus, variations in BrLoc % can be useful in the assessment of the likely effects on the watershed response.The number of bursts can define the effects of rainfalls that are becoming more aggressive and frequent [56,57].TRN and the symmetry can highlight the effect of changes in the rainfall peak intensity and the effects of non-uniformity that can also increase flooding, particularly for short durations [58].
The relationship between hyetographs shapes and ponding time or uNSI was tested using multiple linear regression models considering each criterion as a single predictor, or implying an interaction among them.Statistically, an interaction describes a situation in which the simultaneous influence of two variables on a third is not additive.
The analysis of variance (ANOVA) defines the significance of each model.All the assumptions for the regression models were checked and, in the case of non-linearity between explanatory and response variables, a boxplot transformation was applied; further, the Cook distance [59] was used to remove outliers when needed.We tested the presence of serial correlation among the residuals, and when present, the Durbin-Watson (DW) test [60] at α = 0.05 confirmed or declined the significance of the relationship given by the ANOVA results.The significance threshold level for the ANOVA was set to p < 0.05.

Changes in the Agricultural System and Consequent Changes in Surface Water Management
Modifications of the economy and social development lead to shifts in the water surface management.In the Veneto region, agricultural holdings are mostly an exclusive property of the farmer or of close family members (grandparents, parents), but the number of exclusive properties diminished over the years, from about 80% of farms in 1982 to 69.2% in 2010.Among the family-run business, those with the greater land extend tend to move towards an organised entrepreneurial structure, where exclusive rental or flexible ownership and rental formulas are preferred [37].From 1982 to 2010, the overall UAA decreased of about 10%, but the average UAA per farm increased of about 50% [36], due mostly to the decline of the actual number of holdings.This trend of concentration has been seen in this sector for decades, in a process that dates back to 1982, but which only in the last ten years has been moving at a significantly faster pace [37].The decreasing rate of Veneto businesses was mostly evident for the smaller sized companies.Analysing the phenomenon by UAA category, the number of farms with less than one hectare halved during the years, and only in UAA categories above 20 hectares there was a positive growth trend between the two reference years.
From 1982 to 2010, production crops have remained unchanged: over two-thirds of the land owned by local agricultural holdings are dedicated to arable; ligneous crops (mostly vines) cover almost three-quarters of the total, and at the regional level they remained stable in the percentage ratio of the UAA.However, vineyard areas increased over the last decade, as a direct consequence of the doubling in size of the average land owned by each holding [37].
These changes in the agricultural system brought changes to the management of surface water.Overall, the region registered an increase in the irrigated surface, from ~230,000 ha in 1982 to ~240,000 ha in 2010.This increase is mostly due to changes from non-structural irrigation to organised irrigation with efficient irrigation networks and piping [35].
As found in other areas [15], the selected study site witnessed changes of the surface water management system (Figure 4).1982 to 2010, the overall UAA decreased of about 10%, but the average UAA per farm increased of about 50% [36], due mostly to the decline of the actual number of holdings.This trend of concentration has been seen in this sector for decades, in a process that dates back to 1982, but which only in the last ten years has been moving at a significantly faster pace [37].The decreasing rate of Veneto businesses was mostly evident for the smaller sized companies.Analysing the phenomenon by UAA category, the number of farms with less than one hectare halved during the years, and only in UAA categories above 20 hectares there was a positive growth trend between the two reference years.
From 1982 to 2010, production crops have remained unchanged: over two-thirds of the land owned by local agricultural holdings are dedicated to arable; ligneous crops (mostly vines) cover almost three-quarters of the total, and at the regional level they remained stable in the percentage ratio of the UAA.However, vineyard areas increased over the last decade, as a direct consequence of the doubling in size of the average land owned by each holding [37].
These changes in the agricultural system brought changes to the management of surface water.Overall, the region registered an increase in the irrigated surface, from ~230,000 ha in 1982 to ~240,000 ha in 2010.This increase is mostly due to changes from non-structural irrigation to organised irrigation with efficient irrigation networks and piping [35].
As found in other areas [15], the selected study site witnessed changes of the surface water management system (Figure 4).For the analysed area, the UAA decreased from 586.7 ha in 1982 to 451.9 ha in 2010 [36].In the same timeframe, the overall number of farms decreased, from 342 to 195.Confirming the regional trend, the average UAA per farm increased by about 35% from 1982 to 2010 (from 1.7 ha to 2.3 ha in average).In the study area, the overall extent of crops cultivation decreased, from 214 ha in 1982 to  For the analysed area, the UAA decreased from 586.7 ha in 1982 to 451.9 ha in 2010 [36].In the same timeframe, the overall number of farms decreased, from 342 to 195.Confirming the regional trend, the average UAA per farm increased by about 35% from 1982 to 2010 (from 1.7 ha to 2.3 ha in average).In the study area, the overall extent of crops cultivation decreased, from 214 ha in 1982 to 119 ha in 2010.By considering the regional trend, the irrigated land connected to crops might have diminished from 71 ha in 1982 to about 40 ha in 2010.This decreasing trend is confirmed by the reported value of 66 ha of irrigated land in 1990 [61].These changes in the economic sector seem to be in line with the changes highlighted by the drainage network analysis.The drainage density was 28.2 km/km 2 in 1983 (Figure 4a) in average, and declined to 21.4 km/km 2 in 2010 (Figure 4b).The overall average decrease is about −6.7 km/km 2 (Figure 4c).The storage capacity moved from 172.4 m 3 /ha in 1982 (Figure 4d) to 127.1 m 3 /ha in 2010 (Figure 4e) with an overall average decrease of −44.7 m 3 /ha (Figure 4f).
Despite the general decline, it is interesting to notice in some areas the efficiency of the network increased (blue colours in Figure 4c,f).In some locations, this depended on the creation of engineering projects, such as retention basins and some deep artificial channels.Nevertheless, on a local scale, this was due to different farmer choices, changes in the size of the farms, but also variations in the cultural systems.As noted at the regional level, there was an increase of vine cultivation, doubling from about 140 ha to 300 ha from 1990 to 2010. Figure 5 shows how some areas cultivated with crops in 1983 moved to vines cultivation in 2010.This resulted in the creation of larger channels but with less branching and being mostly located near the fields vs. the in-field network system used for the crop productions in the past (Figure 5).4e) with an overall average decrease of −44.7 m 3 /ha (Figure 4f).Despite the general decline, it is interesting to notice in some areas the efficiency of the network increased (blue colours in Figure 4c,f).In some locations, this depended on the creation of engineering projects, such as retention basins and some deep artificial channels.Nevertheless, on a local scale, this was due to different farmer choices, changes in the size of the farms, but also variations in the cultural systems.As noted at the regional level, there was an increase of vine cultivation, doubling from about 140 ha to 300 ha from 1990 to 2010. Figure 5 shows how some areas cultivated with crops in 1983 moved to vines cultivation in 2010.This resulted in the creation of larger channels but with less branching and being mostly located near the fields vs. the in-field network system used for the crop productions in the past (Figure 5).

Ponding Time
Figure 6 shows the average ponding time (Equation ( 2)) for the considered soil characteristics (Table 1).The highest ponding time (about 30 min) is related to the lowest return period (2 years) and to the dry soil condition (0%-10% saturation).At the increase of the saturation, the ponding time decreases, with the highest changes corresponding to the maximum return period (Figure 6a).

Ponding Time
Figure 6 shows the average ponding time (Equation ( 2)) for the considered soil characteristics (Table 1).The highest ponding time (about 30 min) is related to the lowest return period (2 years) and to the dry soil condition (0%-10% saturation).At the increase of the saturation, the ponding time decreases, with the highest changes corresponding to the maximum return period (Figure 6a).Note that in (a), 0% and 1% saturation overlap, presenting minimal differences.
The ponding time is only slightly influenced by the return period (Figure 6b), with variations only of few minutes depending on the magnitude of the event.Without considering the input rainfall distribution (hyetograph type) and antecedent soil moisture condition (Figure 6b), the considered soil presents a slight increase of the ponding time moving from a rainfall of a 2 years return period to a rainfall of 5 years return period ( + ~2 min).From this point on, increasing the return period results in a decrease of the ponding time: − ~0.4 min from a 5 years to a 10 years' event, or from 10 to 20 years one, and − ~0.8 min moving from a 10 to 30, 30 to 50, 50 to 100 and 100 to 200 years' event.Without considering the return period and averaging each class of antecedent soil moisture condition (Figure 6c), the changes in ponding time are higher and more marked, especially when the saturation conditions are >50%.The 100% saturated soil produces runoff immediately (as expected), while the dry soil produces a runoff about 30 min after the beginning of the rainfall.
The study area soil, for the simulated rainfall amount and duration, is characterised by a specific relationship connecting the event return period or the saturation condition, and the average ponding time, according to: with PT = Ponding Time in hours, x = either the Return period in years or the saturation condition in %, and a and b = empirically derived coefficients (Table 3).Testing the value of R 2 for this equation with a two-tailed Student's t-test for statistical significance at α = 0.05 shows that, in both cases, there is a significant relationship between the two variables for the considered datasets.When considering the saturation condition, the scaling The ponding time is only slightly influenced by the return period (Figure 6b), with variations only of few minutes depending on the magnitude of the event.Without considering the input rainfall distribution (hyetograph type) and antecedent soil moisture condition (Figure 6b), the considered soil presents a slight increase of the ponding time moving from a rainfall of a 2 years return period to a rainfall of 5 years return period ( + ~2 min).From this point on, increasing the return period results in a decrease of the ponding time: − ~0.4 min from a 5 years to a 10 years' event, or from 10 to 20 years one, and − ~0.8 min moving from a 10 to 30, 30 to 50, 50 to 100 and 100 to 200 years' event.Without considering the return period and averaging each class of antecedent soil moisture condition (Figure 6c), the changes in ponding time are higher and more marked, especially when the saturation conditions are >50%.The 100% saturated soil produces runoff immediately (as expected), while the dry soil produces a runoff about 30 min after the beginning of the rainfall.
The study area soil, for the simulated rainfall amount and duration, is characterised by a specific relationship connecting the event return period or the saturation condition, and the average ponding time, according to: with PT = Ponding Time in hours, x = either the Return period in years or the saturation condition in %, and a and b = empirically derived coefficients (Table 3).
Testing the value of R 2 for this equation with a two-tailed Student's t-test for statistical significance at α = 0.05 shows that, in both cases, there is a significant relationship between the two variables for the considered datasets.When considering the saturation condition, the scaling coefficient (a) overlaps with that of Equation ( 4) fitted for the Return period, while the exponent (b) is a larger order of magnitude.

Changes in Water Storage and uNSI
When considering the network changes between 1983 and 2010, and the network response (uNSI), the changes are mostly evident in the dry soil condition (0% saturation, Figure 7), with the highest difference displayed for the lowest return period (2 years).However, differently from the ponding time (Figure 6a), the soil conditions for this return period only matter when the soil is 100% saturated, thus when the drainage network properties completely control the response.For this saturation, there is the smallest difference in the network response between 1983 and 2010 (−7 min).For all the other saturation conditions, differences are quite similar and double the ones of the 100% saturated soil (− ~13 min in all cases).When the return period increases, differences appear more marked between the different soil moisture conditions (Figure 7a).The lowest change in the network response corresponds to the highest return period (Figure 7a) (about −9 min average).
Land 2017, 6, 3 12 of 23 coefficient (a) overlaps with that of Equation ( 4) fitted for the Return period, while the exponent (b) is a larger order of magnitude.

Changes in Water Storage and uNSI
When considering the network changes between 1983 and 2010, and the network response (uNSI), the changes are mostly evident in the dry soil condition (0% saturation, Figure 7), with the highest difference displayed for the lowest return period (2 years).However, differently from the ponding time (Figure 6a), the soil conditions for this return period only matter when the soil is 100% saturated, thus when the drainage network properties completely control the response.For this saturation, there is the smallest difference in the network response between 1983 and 2010 (−7 min).For all the other saturation conditions, differences are quite similar and double the ones of the 100% saturated soil (− ~13 min in all cases).When the return period increases, differences appear more marked between the different soil moisture conditions (Figure 7a).The lowest change in the network response corresponds to the highest return period (Figure 7a) (about −9 min average).The changes in the network response (ΔuNSI), for the considered rainfall and soil, seem to have a specific relationship with the event return period or the saturation condition, following Equation (4) as for the ponding time, but with different values of the coefficients a and b (Table 4).The changes in the network response (∆uNSI), for the considered rainfall and soil, seem to have a specific relationship with the event return period or the saturation condition, following Equation (4) as for the ponding time, but with different values of the coefficients a and b (Table 4).Testing the value of R 2 for this equation with a two-tailed Student's t-test for statistical significance at α = 0.05 shows that, in both cases, there is a significant relationship between the two variables for the considered datasets.As for the ponding time, considering the saturation condition, the scaling coefficient (a) overlaps with that of Equation ( 4) fitted for the return period, while the exponent (b) is of a larger order of magnitude.

Statistical Influence of Climate on Ponding Time and uNSI
Tables S1 and S2 in the supplementary material report the result of the ANOVA analysis about the effect of the rainfall shape, thus the possible effect of climate, on the ponding time and the ∆uNSI.On average, both parameters are influenced by the hydrograph shapes (Figure 8, Tables S1 and S2).The location of the heaviest storm burst (BrLoc % ), the regularity of the rainfall intensity (TRN), and the number of peaks (Br n ) have a statistically significant influence on both ponding time and ∆uNSI.On the other hand, and the asymmetry of the rainfall volume (asym) is significant only for the ponding time.variables for the considered datasets.As for the ponding time, considering the saturation condition, the scaling coefficient (a) overlaps with that of Equation ( 4) fitted for the return period, while the exponent (b) is of a larger order of magnitude.

Statistical Influence of Climate on Ponding Time and uNSI
Tables S1 and S2 in the supplementary material report the result of the ANOVA analysis about the effect of the rainfall shape, thus the possible effect of climate, on the ponding time and the ΔuNSI.On average, both parameters are influenced by the hydrograph shapes (Figure 9, Tables S1 and S2).The location of the heaviest storm burst (BrLoc%), the regularity of the rainfall intensity (TRN), and the number of peaks (Brn) have a statistically significant influence on both ponding time and ΔuNSI.On the other hand, and the asymmetry of the rainfall volume (asym) is significant only for the ponding time.If rainfall events evolved to delay their heaviest peak toward the end of the storm (BrLoc% = 0.9 vs. 0.3), the ponding time could increase more than half an hour (Figure 9a).The differences in the overall watershed response (ΔuNSI) would stay almost the same (main effect ~0), but with the network in 2010 having a slightly faster response than in the past, respect to a rainfall having a burst at the beginning of the storm.The higher regularity of the rainfall intensity during the storm (TRN) only slightly increases the ponding time and emphasises of few minutes the difference between 2010 and 1983.The number of peaks (Brn) has the highest effect, reducing the ponding time greatly (anticipation of about half an hour) and emphasising the difference in the network response of about 20 min).An increase in the volume of If rainfall events evolved to delay their heaviest peak toward the end of the storm (BrLoc % = 0.9 vs. 0.3), the ponding time could increase more than half an hour (Figure 9a).The differences in the overall watershed response (∆uNSI) would stay almost the same (main effect ~0), Land 2017, 6, 3 14 of 24 but with the network in 2010 having a slightly faster response than in the past, respect to a rainfall having a burst at the beginning of the storm.The higher regularity of the rainfall intensity during the storm (TRN) only slightly increases the ponding time and emphasises of few minutes the difference between 2010 and 1983.The number of peaks (Br n ) has the highest effect, reducing the ponding time greatly (anticipation of about half an hour) and emphasising the difference in the network response of about 20 min).An increase in the volume of water before the heaviest burst (asym) slightly reduces the ponding time, while it reduces the differences between the network response in 2010 respect 1983 (Figure 8).
Land 2017, 6, 3 14 of 23 of the rainfall peak (BrLoc% = 0.9 vs. 0.3), is registered if the largest volume of rainfall happens before the burst (asym = 7.4, Figure 9c).However, passing from one storm peak to multiple ones (Brn=1 to 3 in Figure 9a), the ponding time decreases the sooner the heaviest burst happens.Increasing the regularity of the rainfall (TRN = 0.04-0.96)has only slight positive effects on the ponding time, increasing of only a few minutes, with the highest effect related to rainfall having a burst at the beginning of the storm (BrLoc% = 0.3) (Figure 9b).If the largest volume of rainfall happens before the heaviest burst rather than after (asym = 0.5-7.4), and the location of the peak is at the beginning of the storm (BrLoc% = 0.3), ponding time decreases of about half an hour.On the other hand, if the location of the peak is the end of the storm (BrLoc% = 0.9), the ponding time increases of about the same amount (Figure 9c).The increase in the regularity of the rainfall intensity (TRN) (TRN = 0.04-0.96, Figure 9d,e) can have two opposite effects on the ponding time.When interacting with the increase of the peak numbers (Brn), it could delay the ponding time of half to an hour (Figure 9d).However, with a greater water volume falling before the heaviest peak (asym = 7.4 vs. 0.5), the ponding time decreases of up to half an hour (Figure 9d).
Differently, only a limited number of parameters defining the hydrograph shape influences the ΔuNSI (Table S2).Interactions are significant only concerning the asymmetry of the time to the heaviest peak and the regularity of the rainfall (BrLoc%:TRN) (Figure 10b), and the regularity of the rainfall intensity related to how much water volume falls before the heaviest burst (TRN:asym) (Figure 10e).Interactions among the predictors arise.For ponding time (Figure 9a), the changes due to the location of the heaviest burst (BrLoc % ) have the highest effect (+1 hour) related to storms having multiple peaks (Br n = 3) (Figure 9a).Increases in ponding time of about half an hour are related to regularity of the rainfall, whereas the less regular is the rainfall (TRN = 0.04, Figure 9b) the lower is the increase in ponding time (Figure 9b).The highest change in ponding time, connected to a delay of the rainfall peak (BrLoc % = 0.9 vs. 0.3), is registered if the largest volume of rainfall happens before the burst (asym = 7.4, Figure 9c).
However, passing from one storm peak to multiple ones (Br n =1 to 3 in Figure 9a), the ponding time decreases the sooner the heaviest burst happens.Increasing the regularity of the rainfall (TRN = 0.04-0.96)has only slight positive effects on the ponding time, increasing of only a few minutes, with the highest effect related to rainfall having a burst at the beginning of the storm Land 2017, 6, 3 15 of 24 (BrLoc % = 0.3) (Figure 9b).If the largest volume of rainfall happens before the heaviest burst rather than after (asym = 0.5-7.4), and the location of the peak is at the beginning of the storm (BrLoc % = 0.3), ponding time decreases of about half an hour.On the other hand, if the location of the peak is the end of the storm (BrLoc % = 0.9), the ponding time increases of about the same amount (Figure 9c).
The increase in the regularity of the rainfall intensity (TRN) (TRN = 0.04-0.96, Figure 9d,e) can have two opposite effects on the ponding time.When interacting with the increase of the peak numbers (Br n ), it could delay the ponding time of half to an hour (Figure 9d).However, with a greater water volume falling before the heaviest peak (asym = 7.4 vs. 0.5), the ponding time decreases of up to half an hour (Figure 9d).
Differently, only a limited number of parameters defining the hydrograph shape influences the ∆uNSI (Table S2).Interactions are significant only concerning the asymmetry of the time to the heaviest peak and the regularity of the rainfall (BrLoc % :TRN) (Figure 10b), and the regularity of the rainfall intensity related to how much water volume falls before the heaviest burst (TRN:asym) (Figure 10e).If the heaviest storm burst happens at the end of the storm (BrLoc% = 0.9 vs. 0.3), the recent landscape (year 2010) reacts more similarly to the past if the storm comprises multiple bursts (Brn = 3) (Figure 11a).Increasing the regularity of the rainfall (TRN = 0.04 to 0.96) enhances the difference in the network response, with the highest effect related to rainfall having a burst at the beginning of the storm (BrLoc% = 0.3) (Figure 10b).Highly asymmetric rainfalls as well increase the difference in the network response, especially if the largest volume of rainfall happens before the heaviest burst (asym = 0.5-7.4), and the location of the peak is at the end of the storm (BrLoc% = 0.9) (Figure 10c).The interaction between irregularity of the rainfall (TRN) and the number of bursts (Brn) (Figure 10d) highlights how rainfalls with more peaks, but a higher regularity in intensity, increase the difference in the network response.The interaction between the irregularity of the rainfall (TRN) and the asymmetry of the water volume (asym), on the other hand, emphasises the differences between the two years (Figure 10e).For this interaction, the increase of the rainfall intensity regularity produces If the heaviest storm burst happens at the end of the storm (BrLoc % = 0.9 vs. 0.3), the recent landscape (year 2010) reacts more similarly to the past if the storm comprises multiple bursts (Br n = 3) (Figure 11a).Increasing the regularity of the rainfall (TRN = 0.04 to 0.96) enhances the difference in the network response, with the highest effect related to rainfall having a burst at the beginning of the storm (BrLoc % = 0.3) (Figure 10b).Highly asymmetric rainfalls as well increase the difference in the network response, especially if the largest volume of rainfall happens before the heaviest burst (asym = 0.5-7.4), and the location of the peak is at the end of the storm (BrLoc % = 0.9) (Figure 10c).The interaction Land 2017, 6, 3 16 of 24 between irregularity of the rainfall (TRN) and the number of bursts (Br n ) (Figure 10d) highlights how rainfalls with more peaks, but a higher regularity in intensity, increase the difference in the network response.The interaction between the irregularity of the rainfall (TRN) and the asymmetry of the water volume (asym), on the other hand, emphasises the differences between the two years (Figure 10e).For this interaction, the increase of the rainfall intensity regularity produces the most quicker response in 2010, if the largest volume of rainfall falls before the burst itself (asym = 7.4, Figure 10e).that, therefore, can recurrently create issues for the population.The results highlighted that a greater effect of land use (even more quick response in modern times) could be seen at the increase of rainfall bursts within the storm, coupled with the heaviest peaks at the beginning of the event, or a high regularity of the rainfall intensity.As well, it could happen for events highly asymmetric with either the highest peak at the end of the storm or higher regularity of the intensity.Analogous outcomes have already been showed for artificial networks similar to that analysed in this study [15].This information should be considered in the context of changes in the climatic regime.Despite being of minor relevance on the overall climatic dynamic and usually restricted to local occurrences [3,70], high rainfall intensities are becoming more and more frequent [56,58,[90][91][92][93][94][95][96].With warming temperature, a less uniform temporal pattern of precipitation with more intense peak precipitation and weaker precipitation during less intense times can be expected [58].The proposed analysis shows that these ongoing changes might lead to increases in hydraulic criticalities.The saturation conditions and the return period of the event also influence the effects and the interaction of factors.For the ponding time, the asymmetry of rainfall volume (asym) becomes significant only if the rainfall amount exceeds a certain threshold (Return period-Rp ≥ 50 years), independently from the considered saturation (Table S1).The percentage of the total storm period that has occurred at the start of the heaviest burst (BrLoc % ), and the regularity of the rainfall intensity (TRN) and the number of bursts (Br n ) are always significant, for all the considered return periods and saturation conditions (Table S1).However, the number of bursts (Br n ) for the shortest return period (Rp = 2 years) and a dry soil does not significantly influence the ponding time (Table S1).With dry soils and short return period (Sat = 1% and Rp = 2 years), all factors interact significantly influencing the ponding time.They as well interact significantly for an average soil moisture condition (Sat = 50%) and events of a return period of 50 years (Table S1).
For the overall difference in watershed response (∆uNSI), it is important to notice that all factors become significant as a single predictor or interacting with each other when the soil is completely saturated independently from the return period (sat = 100% in Table S2).As well the interaction is always significant for medium to wet soils and short return periods.

Discussion
The proposed analysis showed how agriculture in this area, as in others [15,[62][63][64][65][66][67], has been characterised by continuing intensification resulting in changes in land management and cultivation practices.
These changes are relevant to watershed hydrology because, at larger spatiotemporal scales, they induce significant modifications in evapotranspiration [63], and because of the possible effects of different land management practices for the same land use (e.g., [7,68,69]).In addition to these direct changes, one must consider that the use of arable land and its intensification can also lead to a reduction in soil infiltration rates and available storage capacities, increasing rapid runoff in the form of overland flow [18,[70][71][72].As well, changes might happen due to the creation of subsurface drainage.This work did not account for these points, but it rather focussed on the changes in surface artificial drainage [62][63][64].
In rural floodplain catchments, the presence of linear landscape structures constituting the drainage, and their mutual interactions, strongly affect the hydraulic connectivity and the surface runoff [73][74][75][76][77].The negligible topographic gradients and possible interactions between overland and channel flows constitute the main limitation to predict flood formation, propagation, and inundation [78].As well, changes in arable production have been associated with pressures to work land unsuitable for production [79] or already involved in hydraulic criticalities.
A decrease of the agricultural areas can also be linked to soil sealing [15], and a ratio of the UAA respect to the overall surface under analysis, in this floodplain area, can offer an indirect assessment of the hydrologic vulnerability of a watershed.The variety of hydraulic paths and the availability of volumes invaded that delays the flood wave [80] influence the delays of the flood peak.By keeping the storage volume constant, the discharge peak can increase up to five times when reducing the UAA of about 20% [39].In general, a ratio of the UAA respect to the overall surface under analysis below 50% might present dangerous criticalities, and it should be carefully monitored [39,81], implying possible critical anticipation of the discharge peak and the flood volume.
For the analysed area economic changes led to a change in the farm's organisation, and while the UAA in the past covered more than 50% of the study site, in 2010 the registered UAA results in a coverage of ~25%.As well, economic changes led to a reduction of the irrigated land, mostly due to variations in the size of the fields and on the type of crops cultivated, with an overall decrease of the network density and the connected water storage of about 25%.However, for some areas, the shift from one cultivation to the other (e.g., crops to vineyards), or the need to improve the network efficiency, resulted in a decrease of the network density but a comparable increase of its storage capacity, as also seen in other contexts [82].
Analysing the actual soil response without considering the network changes offers a starting point to understand future changes.Clearly, soil properties play a fundamental role in runoff production, and different soils can result in multiple outcomes [12,83,84].The results highlighted that, for this study area, in the driest soil condition and with common events (short return periods) the overall watershed response is mostly controlled by the ground properties.At the increasing of the saturation, the network characteristics start playing an important role in controlling the hydrologic response.The analysis showed how the magnitude of the rainfall event (return period) only slightly influences the ponding time of the soil of the study area, which is in turn highly influenced by the antecedent soil moisture condition.When saturation exceeds 50%, the soil has a response about 15 min faster.
However, the results highlighted that the intrinsic characteristics of the rainfall event (rain intensity over time) influence the soil response.The proposed analysis suggests that a reduction of the ponding time (thus an expected quickening of the runoff production) is related to rainfall having a large intensity peak at the beginning of the storm, coupled with an increase in the number of storm bursts, or of the volume of water falling before the heaviest burst.As well, changes in the regularity of the intensity during the storm itself coupled with an increase in the number of bursts or the volume of water falling before the heaviest burst might produce a shortening of the ponding time.While at the pedon scale climate might not have effect on the hydrologic response and infiltration [12], these results highlight other conclusions, indicating that at the watershed scale the characteristics of the rainfall, especially the time when the heavy rainfall burst occurs during storm events, influences the soil response [85][86][87].
When considering the network changes and the network response, connected to the soil properties, in maximum soil saturation conditions, network properties completely control the watershed response, and changes are mostly evident for the lowest return period.Thus, the research confirms what found in other studies [15,88,89] that frequent events offer the worst case scenario, that, therefore, can recurrently create issues for the population.The results highlighted that a greater effect of land use (even more quick response in modern times) could be seen at the increase of rainfall bursts within the storm, coupled with the heaviest peaks at the beginning of the event, or a high regularity of the rainfall intensity.As well, it could happen for events highly asymmetric with either the highest peak at the end of the storm or higher regularity of the intensity.Analogous outcomes have already been showed for artificial networks similar to that analysed in this study [15].This information should be considered in the context of changes in the climatic regime.Despite being of minor relevance on the overall climatic dynamic and usually restricted to local occurrences [3,70], high rainfall intensities are becoming more and more frequent [56,58,[90][91][92][93][94][95][96].With warming temperature, a less uniform temporal pattern of precipitation with more intense peak precipitation and weaker precipitation during less intense times can be expected [58].The proposed analysis shows that these ongoing changes might lead to increases in hydraulic criticalities.
Changes are also related to areas where ponding is currently a pressing issue: water ponding can decrease the crop yield and result in a long-term reduction in soil fertility, leading to economic losses and exclusion of large areas from cultivation [97][98][99].The results highlight the importance of monitoring these areas, understanding if changes might imply a further negative impact of ponding.As found in other studies [3,70], the antecedent soil moisture condition is of major importance for the degree to which land-use can influence the watershed response.The return period of the event is less meaningful indicators in this respect.This probably because only storm events originated by high rainfall intensities are at least partially controlled by the conditions of the land-cover and the soil surface [3,70].Studies [100 -102] showed that future climatic conditions affect the water footprint of specific crops, which in turn can affect water management practices (e.g., irrigation techniques) inducing further land use change.As well, in urbanised floodplains, several interacting factors will have induced changes in runoff generation and its delivery to the channel network, e.g., the extent of soil compaction, the efficiency of land drains and the connectivity of flow paths [41].Recent changes in agricultural and flood defence policies create new opportunities for involving the rural land use, in particular, agriculture, in flood risk management [23].Thus, rather than continuing to rely on structural solutions for flood control, it is time to develop a comprehensive flood management strategy that includes hold precipitation where it falls and store flood waters where they occur [103].The proposed analysis can help in this field.

Conclusions
Land use changes, especially in agro-industrial ecosystems where demands for production and residential areas develop simultaneously, deserve considerable attention.These changes influence not only the landscapes in which they occurred, but also water resources at a larger scale, and areas receiving runoff from drained land.As well, climatic changes have always been a primary matter of concern, and the connection between land-use drivers and climate result in people suffering the consequences of severe floods on all continents.The results of this research highlighted that economic trends strictly control the evolvement of artificial drainage system in the considered agro-industrial landscape.These changes are related to effects on the watershed response that depends on the soil condition, network features and rainfall characteristics.In the case of the dry soil condition (0% saturation) and with frequent events (short return periods), management focus should be on the

Figure 3 .
Figure 3. Simulated hyetograph for Sym01 for an event having a return period of 2 years, duration of ~1 h and cumulative rainfall of ~32.5 mm; the effective rainfall producing runoff for different antecedent moisture conditions is shown in red.

Figure 3 .
Figure 3. Simulated hyetograph for Sym01 for an event having a return period of 2 years, duration of ~1 h and cumulative rainfall of ~32.5 mm; the effective rainfall producing runoff for different antecedent moisture conditions is shown in red.

Figure 4 .
Figure 4. Drainage density (Drain den ) and Storage Capacity (Stor cap ) in 1983 and 2010, and relative differences.

Figure 5 .
Figure 5. (a-f) Drainage density (Drainden) (left panel) and Storage Capacity (Storcap) (middle panel) in 1983 and 2010, and overview of the landscape (right panel) for the same years.The figure exemplifies changes in the sizes of the plot, and changes in the network system, and the resulting characteristics of the network.

Figure 5 .
Figure 5. (a-f) Drainage density (Drain den ) (left panel) and Storage Capacity (Stor cap ) (middle panel) in 1983 and 2010, and overview of the landscape (right panel) for the same years.The figure exemplifies changes in the sizes of the plot, and changes in the network system, and the resulting characteristics of the network.

Land 2017, 6 , 3 11 of 23 Figure 6 .
Figure 6.Average ponding time for the considered soil.Ponding time (evaluated considering all hyetograph types) as a function of the Return Period (a) averaged for each of the considered antecedent soil moisture condition (0 to 100% saturation), and (b) averaged considering all saturation conditions.(c) average ponding time considering all hyetograph types and return periods as a function of the antecedent soil moisture condition (saturation).In (b) and (c) different prediction bounds for Equation (4) are shown: BoundOB describes yn+1(x), for all x, while BoundF shows the fitted equation for all x.Note that the intervals associated with a new observation (OB) are wider than the fitted function (F) intervals because of the additional uncertainty in predicting a new response value (the curve plus random errors).Note that in (a), 0% and 1% saturation overlap, presenting minimal differences.

Figure 6 .
Figure 6.Average ponding time for the considered soil.Ponding time (evaluated considering all hyetograph types) as a function of the Return Period (a) averaged for each of the considered antecedent soil moisture condition (0 to 100% saturation), and (b) averaged considering all saturation conditions.(c) average ponding time considering all hyetograph types and return periods as a function of the antecedent soil moisture condition (saturation).In (b) and (c) different prediction bounds for Equation (4) are shown: Bound OB describes y n+1 (x), for all x, while Bound F shows the fitted equation for all x.Note that the intervals associated with a new observation (OB) are wider than the fitted function (F) intervals because of the additional uncertainty in predicting a new response value (the curve plus random errors).Note that in (a), 0% and 1% saturation overlap, presenting minimal differences.

Figure 7 .
Figure 7. Average difference in uNSI between 2010 and 1982.A negative sign indicates a faster response of the network in 2010.Differences in uNSI (evaluated considering all hyetograph types) as a function of the Return period (a) averaged for each of the considered antecedent soil moisture condition (0 to 100% saturation), and (b) averaged considering all saturation conditions.(c) Average difference between uNSI considering all hyetograph types and return periods as a function of the antecedent soil moisture condition (saturation).In (b) and (c) different prediction bounds are shown: BoundOB describes yn+1(x), for all x, while BoundF shows the fitted equation for all x.Note that the intervals associated with a new observation (OB) are wider than the fitted function (F) intervals because of the additional uncertainty in predicting a new response value (the curve plus random errors).

Figure 7 .
Figure 7. Average difference in uNSI between 2010 and 1982.A negative sign indicates a faster response of the network in 2010.Differences in uNSI (evaluated considering all hyetograph types) as a function of the Return period (a) averaged for each of the considered antecedent soil moisture condition (0 to 100% saturation), and (b) averaged considering all saturation conditions.(c) Average difference between uNSI considering all hyetograph types and return periods as a function of the antecedent soil moisture condition (saturation).In (b) and (c) different prediction bounds are shown: Bound OB describes y n+1 (x), for all x, while Bound F shows the fitted equation for all x.Note that the intervals associated with a new observation (OB) are wider than the fitted function (F) intervals because of the additional uncertainty in predicting a new response value (the curve plus random errors).

Figure 8 .
Figure 8.Estimated effects on the response (ponding time (a) and ΔuNSI (b)) from changing each predictor value (percentage of the total storm period that has occurred at the start of the heaviest burst-BrLoc%, number of bursts-Brn, hyetograph triangularity-TRN, asymmetry of the storm-asym), averaging out the effects of the other predictors.Values of the predictors are chosen automatically to produce a relatively large effect on the response.

Figure 8 .
Figure 8.Estimated effects on the response (ponding time (a) and ∆uNSI (b)) from changing each predictor value (percentage of the total storm period that has occurred at the start of the heaviest burst-BrLoc % , number of bursts-Br n , hyetograph triangularity-TRN, asymmetry of the storm-asym), averaging out the effects of the other predictors.Values of the predictors are chosen automatically to produce a relatively large effect on the response.

Figure 9 .
Figure 9.Estimated effects on ponding time of keeping one predictor fixed (percentage of the total storm period that has occurred at the start of the heaviest burst-BrLoc%, number of bursts-Brn, hyetograph triangularity-TRN, asymmetry of the storm-asym) while varying the other.The average effect is shown as a circle, while the horizontal bars are showing the confidence interval for the estimated effect.The blue symbols represent the overall average effect obtained by changing one predictor independently from the other, while the red ones represent the average effect achieved by changing one predictor over different values of the other one.

Figure 9 .
Figure 9.Estimated effects on ponding time of keeping one predictor fixed (percentage of the total storm period that has occurred at the start of the heaviest burst-BrLoc % , number of bursts-Br n , hyetograph triangularity-TRN, asymmetry of the storm-asym) while varying the other.The average effect is shown as a circle, while the horizontal bars are showing the confidence interval for the estimated effect.The blue symbols represent the overall average effect obtained by changing one predictor independently from the other, while the red ones represent the average effect achieved by changing one predictor over different values of the other one.

Land 2017, 6 , 3 15 of 23 Figure 10 .
Figure 10.Estimated effects on ΔuNSI of keeping one predictor fixed (percentage of the total storm period that has occurred at the start of the heaviest burst-BrLoc%, number of bursts-Brn, hyetograph triangularity-TRN, asymmetry of the storm-asym) while varying the other.The average effect is shown as a circle, while the horizontal bars are showing the confidence interval for the estimated effect.The blue symbols represent the overall average effect obtained by changing one predictor independently from the other, while the red ones represent the average effect achieved by changing one predictor over different values of the other one.

Figure 10 .
Figure 10.Estimated effects on ∆uNSI of keeping one predictor fixed (percentage of the total storm period that has occurred at the start of the heaviest burst-BrLoc % , number of bursts-Br n , hyetograph triangularity-TRN, asymmetry of the storm-asym) while varying the other.The average effect is shown as a circle, while the horizontal bars are showing the confidence interval for the estimated effect.The blue symbols represent the overall average effect obtained by changing one predictor independently from the other, while the red ones represent the average effect achieved by changing one predictor over different values of the other one.

Figure 11 .
Figure 11.Local difference in uNSI between 2010 and 1982.A negative sign indicates a faster response of the network in 2010.Differences are evaluated considering all hyetographs and (a-c) all return periods averaged for each specified antecedent soil moisture condition (Sat);(d-f) all saturation conditions averaged for the specified return period (Rp); and (g) displays some characteristics of the urbanisation settings and regulations for the area[27].

Figure 11 .
Figure 11.Local difference in uNSI between 2010 and 1982.A negative sign indicates a faster response of the network in 2010.Differences are evaluated considering all hyetographs and (a-c) all return periods averaged for each specified antecedent soil moisture condition (Sat);(d-f) all saturation conditions averaged for the specified return period (Rp); and (g) displays some characteristics of the urbanisation settings and regulations for the area[27].

Table 2 .
Classification of the hyetographs according to the four defined criteria: location of the heaviest burst within the storm period (BrLoc qtl ), percentage of the total storm period that has occurred at the start of the heaviest burst (BrLoc % ), number of bursts (Br n ), hyetograph triangularity (TRN), asymmetry of the storm (asym).

Table 3 .
Fitted coefficients, 95% confidence intervals (CI 95 ) and goodness of fit (R 2 and Root Mean Square Error-RMSE) for Equation (4) when considering the ponding time (PT) as a function of the return period or the saturation condition.

Table 4 .
Fitted coefficients, 95% confidence intervals (CI 95 ) and goodness of fit (R 2 and Root Mean Square Error-RMSE) for Equation (4) when considering the ∆uNSI (difference between the network response in 2010 vs. 1983) as a function of the return period or as a function of the saturation condition.