Space–time Characterization of Rainfall Field in Tuscany

Precipitation during the period 2001–2016 over the northern and central part of Tuscany was studied in order to characterize the rainfall regime. The dataset consisted of hourly cumulative rainfall series recorded by a network of 801 rain gauges. The territory was divided into 30 × 30 km 2 square areas where the annual and seasonal Average Cumulative Rainfall (ACR) and its uncertainty were estimated using the Non-Parametric Ordinary Block Kriging (NPOBK) technique. The choice of area size was a compromise that allows a satisfactory spatial resolution and an acceptable uncertainty of ACR estimates. The daily ACR was estimated using a less computationally expensive technique, averaging the cumulative rainfall measurements in the area. The trend analysis of annual and seasonal ACR time series was performed by means of the Mann–Kendall test. Four climatic zones were identified: the northwestern was the rainiest, followed by the northeastern , north-central and south-central. An overall increase in precipitation was identified, more intense in the northwest , and determined mostly by the increase in winter precipitation. On the entire territory, the number of rainy days, mean precipitation intensity and sum of daily ACR in four intensity groups were evaluated at annual and seasonal scale. The main result was a magnitude of the ACR trend evaluated as 35 mm/year, due mainly to an increase in light and extreme precipitations. This result is in contrast with the decreasing rainfall detected in the past decades.


Introduction
In recent decades, there has been a gradual increase in the average temperature of the atmosphere.Global warming has determined a higher water vapor concentration and thus an increase in global precipitation [1].Despite this overall scenario, a decrease in rainfall was found in the Mediterranean area; for example, experimental results for Italy and Spain indicated a reduction in precipitation [2,3].Along with the decrease in the average rainfall amount, there has been a change in rainfall patterns, with extreme events becoming more frequent [4]; in particular, Brunetti [5] and Martinez [6] found similar results for Italy and northern Spain (Catalonia).The pattern was not homogeneous on the Italian territory, being more pronounced in the north [7].
Climate change may cause increased drought and extreme events, therefore the study of the evolution of the rainfall regime is a very important task in order to construct a hydrological model, prevent flooding and hydraulic risk [8,9].It is thus necessary to study rainfall phenomena in different time periods.Evaluation of the rainfall regime for relatively long time periods, annual or seasonal, is essential for the study of climate change and hydrological models.Prevention of hydraulic risks requires an assessment of rainfall phenomena for relatively short time periods, daily or hourly.In this perspective, the amount of rain falling on the territory has to be estimated and its uncertainty evaluated in order to establish the reliability of the results.
In this study, the rainfall regime in the Tuscany region was analyzed during the years 2001-2016.A rain gauge network was used that provides measurements of the cumulative rainfall relative to a small area (in the order of hundreds of square centimeters).There were only a limited number of rain gauges in the territory, the typical distance between them was a few km in more populated areas and more than ten km in hilly areas.
The spatial variability of the cumulative rainfall field implies that for a distance of a few km, the values of the correlation coefficient of such an observable is very low [10][11][12].Therefore, by means of a rain gauge network with such a low density, it was not possible to reconstruct the details of the rainfall field using some spatialization techniques.Given the limitation of the instrumentation, an observable was chosen that evaluates the overall rainfall amount on a selected area: the Average Cumulative Rainfall (ACR) i.e., the cumulative rainfall averaged over the area.The spatialization technique chosen to estimate the ACR was Non-Parametric Ordinary Block Kriging (NPOBK) [13][14][15][16][17][18].It allows the best value to be estimated and the standard deviation of the ACR, which is a measure of the uncertainty.
An observable closely related to the ACR is the Rain Volume (RV) i.e., the ACR multiplied by the area; RV can be estimated during thunderstorms by means of the area time integral method, a technique based on non-quantitative rainfall detection [19,20].Griffith et al. [21] used geosynchronous visible or infrared satellite imagery to estimate rainfall over large space and time scales.All these techniques are characterized by lower spatial and temporal resolution.
The estimation of ACR degrades the spatial resolution compared to rain gauge measurements but it allows the amount of rainfall to be assessed on a selected area.In this context, the Tuscany territory was divided into square areas where the ACR was estimated.The uncertainty of the ACR decreases by increasing the area under consideration, but the increase in the area deteriorates the spatial resolution.A compromise was therefore found for the area size that allows a satisfactory spatial resolution and an acceptable uncertainty of ACR estimates.
The annual, seasonal and daily ACR were estimated for all the areas, highlighting the difference in the precipitation amount in different parts of the territory.The trend detection of annual and seasonal ACR time series was performed by means of the Mann-Kendall test.The analysis showed higher rainfall in the north-western area, and an overall increase during the period, more pronounced on the Tyrrhenian coast.This increase occurred mainly in winter.The precipitation increase in Tuscany is contrary to that of previous decades, when a null or negative trend was recorded [22,23].The analysis of daily ACR time series showed an increase in extreme events, therefore in this case, in agreement with the results obtained for the previous years [5,7].

Data
The instrumentation used to evaluate the rainfall field was a network of weather stations equipped with rain gauges and other meteorological sensors.The rain gauges measure cumulative rainfall with an hourly accumulation time.The network is managed by the Hydrological Service of Tuscany Region that provides the quality control of the data.The area covered extends over the Tuscany region and neighboring areas of Emilia-Romagna, Lazio and Umbria.Figure 1  The number of rain gauges increased during the period considered (2001-2016) as new gauges were installed: 325 in 2001 and 801 in 2016.The cumulative rainfall for different accumulation times was evaluated for each rain gauge by adding the hourly measurements.The accumulation times considered were annual (from 1 September to 31 August); seasonal (from 1 September to 30 November in autumn; from 1 December to 28 or 29 February in winter; from 1 March to 31 May in spring; and from 1 June to 31 August in summer); and daily.
Only cumulative rainfall measurements of rain gauges that had recorded more than 90 percent of the hourly measurements during the accumulation time were considered.The results obtained were multiplied by the ratio between the number of hours of the accumulation time and the number of hourly measurements recorded in order not to underestimate the result.This procedure eliminates the bias due to a loss of hourly measurements, but introduces an increase in the uncertainty of the cumulative rainfall.

The NPOBK Technique
The complete characterization of the rainfall regime on the territory can be obtained by evaluating the amount that falls on an arbitrary area for an arbitrary time period; this can be obtained by means of ACR, which can be expressed in terms of rain rate as: where t is an arbitrary instant of time period T, x an arbitrary point inside area A, r the rain rate evaluated in point x and time t, C(x) the cumulative rainfall in point x.Determination of the rain rate field for every instant and every point is a necessary and sufficient condition to determine the ACR for any space-time domain.This is obviously impossible given the limitations of the measurements provided by the available instrumentation; the rain gauges perform cumulative rainfall measurements that represent an integral of the rain rate, furthermore they are only installed at a limited number of points.The first limitation implies that 1 h is the shortest accumulation time to evaluate the ACR.The second implies that the ACR estimates should be made by means of a spatialization of the measurements.The spatialization technique used was the NPOBK method; it is based on the assumption that the cumulative rainfall field C(x) is statistically describable by means of a second-order stationary random function [13][14][15][16][17][18] i.e., the first and second moments are invariant under translations; they are described by means of the following equation: The number of rain gauges increased during the period considered (2001-2016) as new gauges were installed: 325 in 2001 and 801 in 2016.The cumulative rainfall for different accumulation times was evaluated for each rain gauge by adding the hourly measurements.The accumulation times considered were annual (from 1 September to 31 August); seasonal (from 1 September to 30 November in autumn; from 1 December to 28 or 29 February in winter; from 1 March to 31 May in spring; and from 1 June to 31 August in summer); and daily.
Only cumulative rainfall measurements of rain gauges that had recorded more than 90 percent of the hourly measurements during the accumulation time were considered.The results obtained were multiplied by the ratio between the number of hours of the accumulation time and the number of hourly measurements recorded in order not to underestimate the result.This procedure eliminates the bias due to a loss of hourly measurements, but introduces an increase in the uncertainty of the cumulative rainfall.

The NPOBK Technique
The complete characterization of the rainfall regime on the territory can be obtained by evaluating the amount that falls on an arbitrary area for an arbitrary time period; this can be obtained by means of ACR, which can be expressed in terms of rain rate as: where t is an arbitrary instant of time period T, x an arbitrary point inside area A, r the rain rate evaluated in point x and time t, C(x) the cumulative rainfall in point x.
Determination of the rain rate field for every instant and every point is a necessary and sufficient condition to determine the ACR for any space-time domain.This is obviously impossible given the limitations of the measurements provided by the available instrumentation; the rain gauges perform cumulative rainfall measurements that represent an integral of the rain rate, furthermore they are only installed at a limited number of points.The first limitation implies that 1 h is the shortest accumulation time to evaluate the ACR.The second implies that the ACR estimates should be made by means of a spatialization of the measurements.The spatialization technique used was the NPOBK method; it is based on the assumption that the cumulative rainfall field C(x) is statistically describable by means of a second-order stationary random function [13][14][15][16][17][18] i.e., the first and second moments are invariant under translations; they are described by means of the following equation: where E[] is the mean value operator, m the mean value of cumulative rainfall independent of point x, F(h) the correlation function that depends only on the distance h of the two points as the isotropy of the field was also assumed [13].In this perspective, ACR is described by means of a random variable (R_ACR), its mean was considered the best estimate of ACR and its variance a measure of uncertainty.
The technique produces an interpolation function that gives the best unbiased linear estimate of the mean of R_ACR (Equation ( 3)) and its variance (σ 2 ACR ) (Equation ( 4)): where C(x α ) is the rain gauge measurement value on the points x α , A the area, T the accumulation time, γ(x, y) the variogram relative to the pair of points x and y.The parameters λ α and µ are the solution of the following linear equations system: The uncertainty of the ACR estimate was a decreasing function of the number of rain gauges installed within and near the analyzed area, so the size of the area had to be large enough to contain a sufficient number of these in order to perform acceptably accurate estimates.It is worth noting that to estimate ACR, only rain gauge measurements placed within or in proximity to area A were used.The use of measurements away from the area decreases the value of the variance very little at the expense of a substantial increase in the calculation time required.
The variogram, for the isotropy of the model, is a function that depends only on the distance of the points and it can be evaluated by the following equation: C(x) is second-order stationary therefore the second term of the second member is equal to zero (Equation ( 2)).The sample values of the variogram were estimated nonparametrically, i.e., without supposing the shape a priori, evaluating approximately the second member of Equation ( 6) by means of the following: where C(x i ) is the cumulative rainfall measurement of the rain gauge installed in point x i .The sample values of the variogram were determined accurately by considering large values of n(h), therefore using the measurements of rain gauges installed in a large area N. To determine the largest area N where the field was stationary, the variogram was evaluated using measurements of the rain gauges located within increasingly larger areas centered on the same point.When increasing the size of the area, a significant increase was observed in the sill evaluation of the variogram, so the field cannot be considered stationary.Indeed, the apparent growth of the sill is due to the non-stationarity of the field, so in this case, the second term of the second member of Equation ( 6) is not negligible and therefore Equation (7) does not accurately determine the sample values of the variogram.In all the cases studied (Section 3.1), the largest area N where the cumulative rainfall field was stationary contains the area A where the ACR was estimated (a 30 × 30 km 2 square area), therefore the cumulative rainfall is also stationary in area A and the NPOBK method can be applied to estimate the ACR in this area.
To avoid negative values of the R_ACR variance, a conditionally negative defined variogram [15][16][17][18] was assumed.For this reason, the variogram was determined by fitting the sample values with the following exponential function [18]: where Q and r are the fitting parameters.
The rain rate field was always different for each area and in each time as every rainfall event has its own peculiarities, so the characteristics of the accumulation rainfall field differed for each analyzed case.For this reason, for each ACR estimate, the procedure described above was repeated; first the maximum area was determined where the accumulation rainfall field was stationary, then the variogram was calculated using Equations ( 7) and (8).

ACR Estimation
The territory was divided into square areas, the annual, seasonal and daily ACR(A i ,T j ) were evaluated for each area during the period 1 March 2001-31 May 2016.A i denotes an arbitrary area in Figure 2, T j an arbitrary year, season or day.Each area was identified by means of its center, latitude and longitude, expressed in decimal degrees.The annual and seasonal ACR were estimated by means of the NPOBK technique described in Section 2.2.The relative error, defined by 3σ ACR /ACR, was evaluated.Only ACR estimates with relative error lower than 90 percent were considered sufficiently accurate.
Water 2017, 9, 86 5 of 18 cases studied (Section 3.1), the largest area N where the cumulative rainfall field was stationary contains the area A where the ACR was estimated (a 30 × 30 km 2 square area), therefore the cumulative rainfall is also stationary in area A and the NPOBK method can be applied to estimate the ACR in this area.To avoid negative values of the R_ACR variance, a conditionally negative defined variogram [15][16][17][18] was assumed.For this reason, the variogram was determined by fitting the sample values with the following exponential function [18]: where Q and r are the fitting parameters.
The rain rate field was always different for each area and in each time as every rainfall event has its own peculiarities, so the characteristics of the accumulation rainfall field differed for each analyzed case.For this reason, for each ACR estimate, the procedure described above was repeated; first the maximum area was determined where the accumulation rainfall field was stationary, then the variogram was calculated using Equations ( 7) and (8).

ACR Estimation
The territory was divided into square areas, the annual, seasonal and daily ACR(Ai,Tj) were evaluated for each area during the period 1 March 2001-31 May 2016.Ai denotes an arbitrary area in Figure 2, Tj an arbitrary year, season or day.Each area was identified by means of its center, latitude and longitude, expressed in decimal degrees.The annual and seasonal ACR were estimated by means of the NPOBK technique described in Section 2.2.The relative error, defined by 3σACR/ACR, was evaluated.Only ACR estimates with relative error lower than 90 percent were considered sufficiently accurate.The size of the areas had to be large enough to contain a sufficient number of rain gauges in order to perform an acceptably accurate ACR estimate.Conversely, the spatial resolution of the  The size of the areas had to be large enough to contain a sufficient number of rain gauges in order to perform an acceptably accurate ACR estimate.Conversely, the spatial resolution of the accumulation rainfall field is better for small areas.Based on these statements the choice of the size of square areas was 30 × 30 km 2 , which is the best compromise between a good resolution and an acceptable uncertainty of the ACR estimates.Only the northern and central part of Tuscany was examined, as the low density of rain gauges installed in the southern part of the region does not allow the ACR to be estimated with sufficient accuracy (Figure 2).
Over a hundred thousand daily ACR had to be estimated, which is an enormous number.The algorithm to perform the NPOBK is computationally expensive so it was not possible to estimate the daily ACR in a reasonable amount of time.For this reason, the NPOBK technique was not used and another computationally less expensive technique had to be considered.The daily ACR were therefore estimated averaging the cumulative rainfall measurements recorded by rain gauges within the study area.This technique, identified as AM, provides less accurate estimates than the NPOBK technique, and also does not allow the uncertainty to be evaluated.
To establish whether the estimates obtained by means of the AM technique were sufficiently accurate for the aims of this study, the daily ACR estimated by AM and NPOBK techniques was compared for some special cases.The discrepancy of the values in all cases was in the order of a few tens of percentage points, a quantity considered small enough.

Time Average of Local Rainfall
For each area (Figure 2), in order to establish the difference of the typical rainfall amount falling in different zones of Tuscany, the means of the annual and seasonal ACR over the entire period (ACR_AV) were considered: where A i is an arbitrary area (Figure 2), T 1 , . . . .Tn, all the years or all autumn, winter, spring and summer seasons of the considered period, f the joint probability density function of the random variables R_ACR(A i ,T 1 ), . . ., R_ACR(A i ,T N ).
The exact value of ACR_AV could not be established as the joint probability density function was unknown.ACR_AV was therefore evaluated approximately averaging the values of the ACR estimated for each year or season:

Trend Analysis of Local Rainfall
To establish the rainfall regime evolution in each area, in particular if an increase or a decrease in ACR occurred, the trend analysis of annual and seasonal ACR(A i ,T j ) time series was conducted together with an evaluation of the magnitude of the trend.
The trend analysis of the time series was performed by means of the non-parametric Mann-Kendall test for autocorrelated data [24,25].The Z statistic used to perform the trend analysis is a random variable whose probability density function is known and can be approximated to a Gaussian with zero mean and standard deviation equal to 1 for time series with a number of elements greater than or equal to 10.The significance level chosen to reject the null hypothesis, i.e., the absence of a trend, was 0.05, corresponding to Z > 1.96 (positive trend) or Z < −1.96 (negative trend).
The magnitude of the trends was computed using the Theil-Sen estimator (TSA), which is less sensitive to outliers than the least square linear fitting [26,27].

Global Rainfall
Annual, seasonal and daily ACR evaluated over the entire territory (ACR_TOT) was analyzed in order to study the overall evolution of the rainfall regime.The cumulative rainfall field was not second-order stationary as the area was too large, so the NPOBK technique was not applicable.It was therefore chosen to evaluate ACR_TOT by means of the mean of ACR in all areas (Figure 2) corresponding to the same accumulation time: where T j is an arbitrary year for annual ACR_TOT, or an arbitrary autumn, winter, spring, summer season for seasonal ACR_TOT, or an arbitrary day for daily ACR_TOT; A 1 . . ... A m are all the areas in Figure 2, f the joint probability density function of random variables R_ACR(A 1 ,T j ), . . ., R_ACR(A M ,T j ).
Once again, the joint probability density function of the R_ACR was unknown so the ACR_TOT was evaluated, approximately averaging the ACR estimates;

Observables for the Study of the Global Rainfall Regime
The following observables were analyzed on annual and seasonal time scales, [27,28]: ACR_TOT.Number of rainy days (NRD): the number of days during which daily ACR was greater than 1 mm.If the ACR was less than 1 mm, the day was considered not rainy.
The time scale will be indicated in subscript to the name of the observable (ann, aut, win, spr, sum indicate the annual and seasonal time scales respectively).
It should be noted that the daily ACR evaluated over relatively large areas never presents the peak of daily cumulative rainfall values observed by some rain gauges, which take measurements on a much smaller area.The correlation coefficient of daily cumulative rainfall presents values significantly lower than 1, even for distances of a few km [10][11][12].This implies that the peak values may not remain on the whole area.For this reason, the threshold value (40 mm) chosen for the last intensity group is significantly lower than the cumulative rainfall measurements of some rain gauges during an extreme event, which may exceed 100 mm.
The trend analysis was performed for all time series and the magnitude of the trends was computed.

Evaluation of the Variogram
The sample values of variograms for all the years and seasons were derived by means of Equation (7), using the cumulative rainfall measurements inside 30 × 30, 40 × 40, 50 × 50, 60 × 60, 70 × 70 km 2 square areas.The variograms were evaluated by fitting the sample values obtained by the measurements within the largest areas where the condition of stationarity occurred.In almost every case, the stationarity of the field is checked inside areas of 60 × 60 km 2 or smaller.In all cases, the value of the SILL (Q parameter of Equation ( 8)) was between 40,000 and 600,000 mm 2 for the annual ACR and between 3000 and 100,000 mm 2 for the seasonal ACR.The r parameter (Equation ( 8)) was between 7000 and 30,000 m for the annual ACR and 15,000 and 70,000 for the seasonal ACR.As an example, Figure 3 reports the variograms of annual and seasonal ACR for the year 2009-2010 (areas centered at lat. 43.75 • , lon.11.25 • ).In all cases, the 60 × 60 km 2 area was the largest where the stationarity condition was verified, except in spring 2010, where this condition was also verified for the 70 × 70 km 2 area.

Evaluation of the Variogram
The sample values of variograms for all the years and seasons were derived by means of Equation (7), using the cumulative rainfall measurements inside 30 × 30, 40 × 40, 50 × 50, 60 × 60, 70 × 70 km 2 square areas.The variograms were evaluated by fitting the sample values obtained by the measurements within the largest areas where the condition of stationarity occurred.In almost every case, the stationarity of the field is checked inside areas of 60 × 60 km 2 or smaller.In all cases, the value of the SILL (Q parameter of Equation ( 8)) was between 40,000 and 600,000 mm 2 for the annual ACR and between 3000 and 100,000 mm 2 for the seasonal ACR.The r parameter (Equation ( 8)) was between 7000 and 30,000 m for the annual ACR and 15,000 and 70,000 for the seasonal ACR.As an example, Figure 3 reports the variograms of annual and seasonal ACR for the year 2009-2010 (areas centered at lat. 43.75°, lon.11.25°).In all cases, the 60 × 60 km 2 area was the largest where the stationarity condition was verified, except in spring 2010, where this condition was also verified for the 70 × 70 km 2 area.

Estimation of Rainfall Amount
The annual and seasonal ACR were estimated for each area (Figure 2) by means of the NPOBK technique during the considered period.Examples are given in Table 1 (for some areas and all years and seasons) and in Figures 4 and 5 (for all areas and some years and seasons).The typical values of the relative error were some tens of percentage points; in all cases, it was lower than 75 percent.The density of rain gauges was not homogeneous, it was greater near the main towns and in the most populated areas in general (Figure 1), therefore the relative error decreases accordingly.New rain gauges were installed during the considered period, so in some areas the uncertainty was lower when evaluated more recently.For some seasons and areas, it was not possible to estimate the ACR since the measurements acquired by rain gauges were less than 90 percent of the total.It should be noted that, for each year, the sum of seasonal ACR estimates is not exactly equal to the annual one.This discrepancy is due to two distinct causes.First, the data quality control sometimes excluded the measurements of some rain gauges; the rain gauge measurements excluded by the annual estimates may differ from the seasonal ones.Therefore, the set of data in the interpolation equations (Equation ( 3)) to estimate annual and seasonal ACR were different.Second, the annual and seasonal cumulative rainfall fields and thus the shape of the related variograms did not coincide, so the coefficients of the NPOBK linear systems (Equation ( 5)) related to the same rain gauges were different.
The daily ACR were estimated by the AM technique.Table 2 reports examples of estimates obtained by means of the NPOBK and AM techniques.The results show a discrepancy in the order of some tens of percentage points, which is sufficiently accurate.The annual and seasonal ACR_AV of each area was calculated (Equation ( 10)); the results are reported in Table 3 and Figure 6.The annual values of ACR_AV highlight different rainfall regimes for different areas.In central-southern areas (centered at lat. 43.21 • and 43.48 • ) ACR_AV were between 800 and 900 mm, in central-northern areas (centered at lat. 43.75 • ) they were about 1000 mm.In north-western areas (centered at lat. 44.01 • , 44.28 • and lon.10.12 • , 10.87 • ) ACR_AV were about 1600 mm, in north-eastern areas (lat.44.01 • , lon.11.25 • , 11.62 • ) 1100 mm.
In spring, the ACR_AV values in central-southern areas were about 210 mm.In central-northern areas a gradient in the east-west direction was observed; in the east, the ACR values were about 200 mm, while in the west 280 mm.In north-western areas, they were about 350 mm, while in north-eastern areas 300 mm.In the summer, the ACR_AV values in central-southern and central-northern areas were about 130 mm, and in northern areas 200 mm.In autumn, in central-southern areas, a gradient in the east-west direction was observed; in the west, the ACR_AV values were about 350 mm, in the east 250 mm.In central-northern and northern areas, values were about 350 mm and 500 mm respectively.In the winter, the ACR_AV values in central-southern areas were about 250 mm, in central-northern areas about 300 mm and 500 mm in northern areas.
Based on these results, it may be noted that the annual volume of rainfall decreases moving in the north-south direction; furthermore, the north-western part of the territory was rainier than the north-east.The rainfall regime, on a seasonal basis, reproduces the same pattern; if a region is rainier annually, the same pattern occurs for all seasons.In spring, the ACR_AV values in central-southern areas were about 210 mm.In central-northern areas a gradient in the east-west direction was observed; in the east, the ACR values were about 200 mm, while in the west 280 mm.In north-western areas, they were about 350 mm, while in northeastern areas 300 mm.In the summer, the ACR_AV values in central-southern and central-northern areas were about 130 mm, and in northern areas 200 mm.In autumn, in central-southern areas, a gradient in the east-west direction was observed; in the west, the ACR_AV values were about

Evolution of the Local Rainfall Regime
The trend analysis was performed for the time series of annual and seasonal ACR of each area (Figures 2 and 6, Table 3).For almost all areas, the value of the Z statistic of the Mann-Kendall test was positive, which is a clue to a positive trend.The Z value of many time series was lower than 1.96, the threshold chosen to reject the null hypothesis, therefore it is not possible to state with sufficient confidence the presence of a positive trend in many of the areas considered.
The magnitude of trends of annual ACR was higher in the north-western coastal areas; in general the gradient of the trend magnitude was directed towards the north-west; in other words, there is a Water 2017, 9, 86 14 of 18 trend increase by moving in the western and northern directions.The typical trend magnitude was some tens of mm/year.In autumn, winter and summer, the trend analysis highlights a similar pattern to the annual one; in spring, the magnitude of the trend was overall homogeneous.

Characterization of the Global Rainfall Regime
The annual and seasonal values of ACR_TOT, NRD, MIP, SD0, SD10, SD25, SD40 were evaluated (Figures 7 and 8).The Mann-Kendall test was performed for all time series (Table 4).
The trend analysis was performed for the time series of annual and seasonal ACR of each area (Figures 2 and 6, Table 3).For almost all areas, the value of the Z statistic of the Mann-Kendall test was positive, which is a clue to a positive trend.The Z value of many time series was lower than 1.96, the threshold chosen to reject the null hypothesis, therefore it is not possible to state with sufficient confidence the presence of a positive trend in many of the areas considered.
The magnitude of trends of annual ACR was higher in the north-western coastal areas; in general the gradient of the trend magnitude was directed towards the north-west; in other words, there is a trend increase by moving in the western and northern directions.The typical trend magnitude was some tens of mm/year.In autumn, winter and summer, the trend analysis highlights a similar pattern to the annual one; in spring, the magnitude of the trend was overall homogeneous.

Characterization of the Global Rainfall Regime
The annual and seasonal values of ACR_TOT, NRD, MIP, SD0, SD10, SD25, SD40 were evaluated (Figures 7 and 8).The Mann-Kendall test was performed for all time series (Table 4).ACR_TOT ann highlights an increase in precipitation; the Z statistic value of 1.86 was very close to the threshold value, suggesting the presence of a growing trend.The trend magnitude of 35.4 mm/year indicates a significant increase in the rainfall amount.The Z values of ACR_TOT aut , ACR_TOT win , ACR_TOT spr , ACR_TOT sum show an increase in precipitation only in winter.
The trend analysis of the NRD ann time series does not clearly highlight an increase in this observable.However, the Z values of seasonal time series show that a positive trend was detected only in winter.
The Z value of MIP ann was 1.86; also in this case, the presence of a positive trend is very likely.The analysis of the Z statistic of the seasonal time series did not show an increase in this observable; although the Z value of MIP win suggests the possible presence of a positive trend in winter.
The study of these observables highlights that one of the causes of the increase in annual precipitation was the growth of MIP ann .The Mann-Kendall test of NRD ann does not clearly show a positive trend, although the Z value might suggest it; therefore, it cannot be excluded that this phenomenon could also have caused the increase in annual rainfall.For all the considered time series, the positive trend was always attributable to a change of the rainfall regime in winter; a clear growth of the observables was not highlighted in the other seasons.
The trend analysis of SD0 ann shows a positive trend; the value of the Z statistic was 2.36.The magnitude of trend equal to 13.3 mm/year was very high, so the rainfall amount caused by low-intensity events has increased a lot during the considered period.SD10 ann and SD25 ann did not show any positive trend, so the rainfall amount due to medium intensity events can be considered stationary.The value of the Z statistic of SD40 ann was 2.58, so the presence of a positive trend was established.The magnitude of the trend equal to 10.5 mm/year was very high; it shows a significant increase in the contribution of extreme events to the rainfall amount.
The analysis of seasonal time series shows that the rainfall regime in autumn, winter and summer was similar to the annual one; there was an increase in rainfall amount for the light and extreme precipitations and there was substantial stability for the intermediate ones.In the spring, no trend was observed in any of the intensity groups.

Discussion and Conclusions
The estimation of rainfall amount is fundamental to evaluate the rainfall regime on the territory.The observable introduced to perform this task was the ACR [13,14], evaluated on some square areas in the studied territory.In addition to the best value of the observable, it is very important to quantitatively evaluate the uncertainty to establish the validity of the results.The annual and seasonal ACR and its uncertainty was therefore estimated, spatializing the measurements of a rain gauge network by means of the NPOBK technique [13][14][15][16][17][18].The daily measurements were evaluated without an assessment of the error.
The uncertainty was due mainly to the limited density of rain gauges on the territory; indeed, the standard deviation of the ACR estimate depends essentially on the number of rain gauges in the area.In this context, a compromise was reached about the size of the square areas (30 × 30 km 2 ) that Water 2017, 9, 86 16 of 18 allows a satisfactory spatial resolution and an acceptable uncertainty.Only the northern and central part of Tuscany was studied as the low density of rain gauges in the south does not allow acceptable estimates to be obtained.The period considered was 1 March 2001-31 May 2016.
The results can be summarized as follows: the territory can be divided into four areas distinguished by the ACR_AV values during the considered period.The north-western area is the rainiest (ACR_AV about 1600 mm) and experienced the highest increase in precipitation.The north-eastern area (ACR_AV about 1200 mm), central area (ACR_AV about 1000 mm) and southern-central area (ACR_AV about 800 mm) had the most limited increase in precipitation.For all areas, the precipitation in autumn and winter was more abundant than in spring; summer precipitation was the least abundant.
The ACR_TOT time series showed a very strong increase in precipitation; this was assessed as 35 mm/year and was mainly due to winter rainfall.It was caused by the positive trend in MIP and probably also in NRD.The ACR_TOT positive trend was not homogenous on the territory but the magnitude of the trend was higher in the north-western area.Analysis of the distribution of daily precipitation shows an increase in the contribution of low and extreme intensity events, while the contribution of medium intensity events can be considered stationary.
The increase in rainfall that occurred during the considered period is in contrast with the overall trend of the last decades.Bartolini et al. [23] showed a tendency to a decrease in rainfall during the period 1948-2009 in two distinct sites in Tuscany, while Fatichi et al. [22], during the period 1916-2003, showed an absence of any trends in precipitation amount or intensity of extreme events.
The results of this paper are only apparently surprising.Within a generally decreasing rainfall trend, a period of precipitation increase may occur without affecting the overall trend.For example, Romano et al. [28] showed an overall decreasing trend in annual precipitation in the Tiber river basin, an area close to Tuscany.The details of the pattern of the rainfall time series showed a succession of dry and wet periods when the rainfall amount increased.
The values of SD40 have increased, so extreme rainfall events have become more frequent.This result is in agreement with Crisci et al. [29] who reported an increase of extreme events in Tuscany during the period 1973-1994, while Fatichi et al. [22] did not detect any trend on the basis of a longer time series.The increase in extreme precipitation aggravates the risk of flooding; in this context, the monitoring of these events, also with the technique described in this paper, becomes essential.
In the future, the ACR estimates will be improved by means of the measurements from radar operating on the territory.The radar can measure the average rain rate on 1 × 1 km square pixels, so can improve the spatial resolution of ACR estimates.The radar sweeps all Tuscany, so the measurements will allow the rainfall regime to be studied over the whole region.

Figure 1 .
Figure 1.The Tuscany rain gauge network.(a) shows the position of the Tuscany region on the Italian peninsula; (b) the rain gauges installed in the region (March 2015).

Figure 1 .
Figure 1.The Tuscany rain gauge network.(a) shows the position of the Tuscany region on the Italian peninsula; (b) the rain gauges installed in the region (March 2015).

Figure 2 .
Figure 2. Areas where the ACR was estimated.Latitude and longitude of the centers of the areas expressed in decimal degrees are indicated below and above them.

Figure 2 .
Figure 2. Areas where the ACR was estimated.Latitude and longitude of the centers of the areas expressed in decimal degrees are indicated below and above them.

Figure 3 .
Figure 3. Variograms of annual (a) and seasonal (b,c,d,e) ACR (year: 2009-2010, areas centered at lat. 43.75°, lon.11.25°).The sample values obtained with the measurements inside 50 × 50, 60 × 60, and 70 × 70 km 2 areas are reported.The distances of the points (meters) are reported in the abscissa, the values of the variogram (mm 2 ) in the ordinate.The fitting function is also reported.

Figure 3 .
Figure 3. Variograms of annual (a) and seasonal (b-e) ACR (year: 2009-2010, areas centered at lat. 43.75 • , lon.11.25 • ).The sample values obtained with the measurements inside 50 × 50, 60 × 60, and 70 × 70 km 2 areas are reported.The distances of the points (meters) are reported in the abscissa, the values of the variogram (mm 2 ) in the ordinate.The fitting function is also reported.

Figure 6 .
Figure 6.Annual and seasonal ACR_AV estimates (a,d,g,j,m); values of Z statistic (b,e,h,k,n); and magnitudes of the trend (c,f,i,l,o).

Figure 6 .
Figure 6.Annual and seasonal ACR_AV estimates (a,d,g,j,m); values of Z statistic (b,e,h,k,n); and magnitudes of the trend (c,f,i,l,o).

Table 1 .
Examples of annual and seasonal ACR (mm) and their standard deviations (σ) and relative errors (ERR).M denotes missing values.Latitude and longitude are expressed in decimal degrees.

Table 2 .
Comparison of daily ACR estimates obtained by means of AM and NPOBK techniques.
3.3.Characterization of the Local Rainfall Regime3.3.1.Estimation of the Time Average of the Local Rainfall

Table 3 .
Annual and seasonal ACR_AV estimates, values of Z statistic and magnitudes of the trend.The Z values >1.96 are reported in bold.

Table 4 .
Values of Z statistic and trend magnitudes of annual and seasonal observables.The Z values > 1.96 are reported in bold.ANNUAL AUTUMN WINTER SPRING SUMMER Z T Z T Z T Z T Z T ACR_TOT 1.86 35.4 0.49 2.9 3.57 14.8 0.20 1.1 1.48 3.9

Table 4 .
Values of Z statistic and trend magnitudes of annual and seasonal observables.The Z values > 1.96 are reported in bold.