Statistical and Fractal Approaches on Long Time-Series to Surface-Water/Groundwater Relationship Assessment: A Central Italy Alluvial Plain Case Study

: In this research, univariate and bivariate statistical methods were applied to rainfall, river and piezometric level datasets belonging to 24-year time series (1986–2009). These methods, which often are used to understand the effects of precipitation on rivers and karstic springs discharge, have been used to assess piezometric level response to rainfall and river level ﬂuctuations in a porous aquifer. A rain gauge, a river level gauge and three wells, located in Central Italy along the lower Pescara River valley in correspondence of its important alluvial aquifer, provided the data. Statistical analysis has been used within a known hydrogeological framework, which has been reﬁned by mean of a photo-interpretation and a GPS survey. Water–groundwater relationships were identiﬁed following the autocorrelation and cross-correlation analyses. Spectral analysis and mono-fractal features of time series were assessed to provide information on multi-year variability, data distributions, their fractal dimension and the distribution return time within the historical time series. The statistical–mathematical results were interpreted through ﬁeldwork that identiﬁed distinct groundwater ﬂowpaths within the aquifer and enabled the implementation of a conceptual model, improving the knowledge on water resources management tools.


Introduction
The careful and correct management of water resources in densely populated alluvial areas is becoming necessary because of the groundwater depletion caused by increased exploitation associated with demographic growth. Previous scientific research on long time series was focused mainly on hydrogeological processes concerning relationships between precipitation and spring discharge or river level, quite often considering fractured aquifer systems [1][2][3][4][5][6]; these analyses provided complementary information to the hydrogeological properties of a natural system, to the transit time through the subsurface's unsaturated zone, and to the recharge of the aquifer, supporting the correct implementation of hydrogeological conceptual and numerical models [7][8][9][10]. On the other hand, the behavior analysis of long piezometric level time series was almost neglected and limited to short periods that rarely exceeded a year (e.g., contaminated site studies). Recently, Cai and Ofterdinger [11] evaluated the statistical relationships between two-year rainfall and piezometric  Table 1 for initials and abbreviations).

Study Area
The studied alluvial aquifers are located between the Apennine Mountains and the Adriatic Sea. The more elevated areas in this sector are constituted of Mio-Pliocenic terrigenous pelitic-arenaceous and pelitic-marly deposits [25][26][27]. The alluvial aquifers are sedimentary structures that include various orders of terraces (I to IV order) and the current alluvial plain. These deposits are mostly silty and sandy in nature with gravelly-sandy lenses of the period between medium-upper Pleistocene and Holocene age. The spatio-temporal evolution of the fluvial depositional environments, driven by interactions between climatic and glacio-eustatic changes and the plio-quartenary tectonic phases, has conditioned the spatial distribution of the characteristic lithotypes [28]. The presence of a mostly pelitic substrate, acting as an aquiclude (k values ranging from 10 −9 to 10 −7 m/s), below the more permeable alluvial deposits (k values ranging from 10 −5 to 10 −3 m/s) suggests the presence of shallow, single and multi-layered aquifers [29,30].
These aquifers often interact with the streams and rivers through the hyphorheic zone. Hydraulic connection between the stream and the groundwater can occur with a transfer of mass  Table 1 for initials and abbreviations).

Study Area
The studied alluvial aquifers are located between the Apennine Mountains and the Adriatic Sea. The more elevated areas in this sector are constituted of Mio-Pliocenic terrigenous pelitic-arenaceous and pelitic-marly deposits [25][26][27]. The alluvial aquifers are sedimentary structures that include various orders of terraces (I to IV order) and the current alluvial plain. These deposits are mostly silty and sandy in nature with gravelly-sandy lenses of the period between medium-upper Pleistocene and Holocene age. The spatio-temporal evolution of the fluvial depositional environments, driven by interactions between climatic and glacio-eustatic changes and the plio-quartenary tectonic phases, has conditioned the spatial distribution of the characteristic lithotypes [28]. The presence of a mostly pelitic substrate, acting as an aquiclude (k values ranging from 10 −9 to 10 −7 m/s), below the more permeable alluvial deposits (k values ranging from 10 −5 to 10 −3 m/s) suggests the presence of shallow, single and multi-layered aquifers [29,30].
These aquifers often interact with the streams and rivers through the hyphorheic zone. Hydraulic connection between the stream and the groundwater can occur with a transfer of mass and/or pressure head, based on the local hydrogeological conditions [31]. The alluvial aquifers are fed by the limestone aquifers of the Apennine Mountains in the internal areas [32][33][34][35][36], while recharge from rainfall has a major role in the total recharge to the shallow alluvial aquifers moving downstream towards the sea, particularly in correspondence of the terraced deposits that are hydraulically connected to the lower valleys [37][38][39][40]. The wells considered in this study fall in the terminal areas of the Pescara River alluvial valley, near the Adriactic Sea.
The Pescara River has a basin of 3190 km 2 (Figure 1), which can be subdivided in two main parts: the first one is the western mountainous sector, and the second one is the eastern sector, where the alluvial aquifer under study is located. The Pescara River mean discharge measured in the Pescara Santa Teresa gauge (PST) is about 50 m 3 /s (31 m 3 /s comes from the western part of its basin).
The alluvial deposits making up the aquifer present four orders of terraces, which have different elevations with respect to the base of the valley. The first and second order terraces are at higher elevations and are mainly constituted of gravelly deposits with silty-sandy lenses; the third and four order terraces are at lower elevations and are constituted of gravelly deposits with extensive silty-sandy and silty-clayey lenses. Only the fourth order terrace runs on both sides of the river, with thin gravelly-sandy deposits, considered to be the current alluvial deposits, being directly connected to this terrace ( Figure 2). The overall thickness of the quaternary alluvial deposits increases moving downstream. Thickness increases from 12 to 16 m in the middle section of the Pescara River valley to 40-50 m near the coast. As the thickness increases, the nature of the deposits varies from mainly gravelly to silty-sandy.
Water 2017, 9, 850 4 of 28 and/or pressure head, based on the local hydrogeological conditions [31]. The alluvial aquifers are fed by the limestone aquifers of the Apennine Mountains in the internal areas [32][33][34][35][36], while recharge from rainfall has a major role in the total recharge to the shallow alluvial aquifers moving downstream towards the sea, particularly in correspondence of the terraced deposits that are hydraulically connected to the lower valleys [37][38][39][40]. The wells considered in this study fall in the terminal areas of the Pescara River alluvial valley, near the Adriactic Sea.
The Pescara River has a basin of 3190 km 2 (Figure 1), which can be subdivided in two main parts: the first one is the western mountainous sector, and the second one is the eastern sector, where the alluvial aquifer under study is located. The Pescara River mean discharge measured in the Pescara Santa Teresa gauge (PST) is about 50 m 3 /s (31 m 3 /s comes from the western part of its basin).
The alluvial deposits making up the aquifer present four orders of terraces, which have different elevations with respect to the base of the valley. The first and second order terraces are at higher elevations and are mainly constituted of gravelly deposits with silty-sandy lenses; the third and four order terraces are at lower elevations and are constituted of gravelly deposits with extensive silty-sandy and silty-clayey lenses. Only the fourth order terrace runs on both sides of the river, with thin gravelly-sandy deposits, considered to be the current alluvial deposits, being directly connected to this terrace ( Figure 2). The overall thickness of the quaternary alluvial deposits increases moving downstream. Thickness increases from 12 to 16 m in the middle section of the Pescara River valley to 40-50 m near the coast. As the thickness increases, the nature of the deposits varies from mainly gravelly to silty-sandy. The groundwater flow direction follows that of the river in the most upstream areas of the valley, although groundwater flow will be dependent on the local spatial distribution of the lithotypes (presence of paleostreams, etc.) and on the morphology of the substrate, which will be connected to erosional history and tectonic deformations. The groundwater flow is toward the river moving downstream where recharge to the aquifer is clearly originating from the terraces (Figure 1). The morphology of the substrate has a strong role in controlling groundwater flow, since groundwater divides give rise to secondary flow circuits within the aquifer. The average hydraulic gradient of the shallow groundwater below the river bed varies from 4‰ to 6‰ [29,41]. A meandering hydrographic pattern is observed in the area, which, together with the changes in local river characteristics (i.e., river head slope and river width variations), controls the relationship between groundwater and surface water (Figures 1 and 3).
The geometric characteristics of the river section under investigation, assessed by means of aerial photo-interpretation and in-situ GPS surveys, indicate that the width of the river is quite variable. However, in general (Figure 3), the width increases moving downstream. There are two discontinuities along the river profile, which correspond to river dams for hydroelectric purposes, and the river head slope is very low upward of these artificial impoundments. These hydroelectric The groundwater flow direction follows that of the river in the most upstream areas of the valley, although groundwater flow will be dependent on the local spatial distribution of the lithotypes (presence of paleostreams, etc.) and on the morphology of the substrate, which will be connected to erosional history and tectonic deformations. The groundwater flow is toward the river moving downstream where recharge to the aquifer is clearly originating from the terraces ( Figure 1). The morphology of the substrate has a strong role in controlling groundwater flow, since groundwater divides give rise to secondary flow circuits within the aquifer. The average hydraulic gradient of the shallow groundwater below the river bed varies from 4 to 6 [29,41]. A meandering hydrographic pattern is observed in the area, which, together with the changes in local river characteristics (i.e., river head slope and river width variations), controls the relationship between groundwater and surface water (Figures 1 and 3).
The geometric characteristics of the river section under investigation, assessed by means of aerial photo-interpretation and in-situ GPS surveys, indicate that the width of the river is quite variable. However, in general (Figure 3), the width increases moving downstream. There are two discontinuities along the river profile, which correspond to river dams for hydroelectric purposes, and the river head slope is very low upward of these artificial impoundments. These hydroelectric plants were constructed after the period in which the data was collected and therefore had no influence on the collected datasets.
The selected wells are located on the right respect to the river bank, at variable distances from the river bed ( Figure 4, Table 1). These wells have been drilled to total depth through silty-sandy deposits and their exploitation is sporadic and limited to backyard irrigation during summer season.
Water 2017, 9, 850 5 of 28 plants were constructed after the period in which the data was collected and therefore had no influence on the collected datasets. The selected wells are located on the right respect to the river bank, at variable distances from the river bed ( Figure 4, Table 1). These wells have been drilled to total depth through silty-sandy deposits and their exploitation is sporadic and limited to backyard irrigation during summer season.   plants were constructed after the period in which the data was collected and therefore had no influence on the collected datasets. The selected wells are located on the right respect to the river bank, at variable distances from the river bed ( Figure 4, Table 1). These wells have been drilled to total depth through silty-sandy deposits and their exploitation is sporadic and limited to backyard irrigation during summer season.

Datasets
Datasets used in this study, made available by the Hydrographic Service of Abruzzo Region, refer to "Spoltore" automatic rain gauge (Sp), "Pescara S. Teresa" hydrometer (PST) and three wells (Figures 1 and 2, and Table 1), "Surricchio" (Su), "Sanità" (Sa) and "De Nicola" (DN). While rainfall data were available on a daily basis, the piezometric and hydrometric levels instead were measured every three days, which means that, for every month, 10 measures were made available. For this reason, to compare the time series, daily rainfall datasets where cumulated on a three-day basis. Since the data available are relative to the 24-year period 1986-2009, each time series consisted of 2880 records. Well characteristics and stratigraphy were identified from previous geognostic and field surveys ( Table 1, and Figure 4). The sea level time series has not been considered because of the distance of each monitoring point from the sea ( Table 1).
The three-day time series are plotted against time in Figure 5. River level and rainfall data monitored at the Pescara S. Teresa (PST) and Spoltore (Sp) sites, and piezometric heads measured from the Sa and DN wells show that the variability appears to be mainly seasonal, while at first sight a plurennial cyclic behavior is evident in the piezometric level of Su well.

Datasets
Datasets used in this study, made available by the Hydrographic Service of Abruzzo Region, refer to "Spoltore" automatic rain gauge (Sp), "Pescara S. Teresa" hydrometer (PST) and three wells (Figures 1 and 2, and Table 1), "Surricchio" (Su), "Sanità" (Sa) and "De Nicola" (DN). While rainfall data were available on a daily basis, the piezometric and hydrometric levels instead were measured every three days, which means that, for every month, 10 measures were made available. For this reason, to compare the time series, daily rainfall datasets where cumulated on a three-day basis. Since the data available are relative to the 24-year period 1986-2009, each time series consisted of 2880 records. Well characteristics and stratigraphy were identified from previous geognostic and field surveys ( Table 1, and Figure 4). The sea level time series has not been considered because of the distance of each monitoring point from the sea ( Table 1).
The three-day time series are plotted against time in Figure 5. River level and rainfall data monitored at the Pescara S. Teresa (PST) and Spoltore (Sp) sites, and piezometric heads measured from the Sa and DN wells show that the variability appears to be mainly seasonal, while at first sight a plurennial cyclic behavior is evident in the piezometric level of Su well.

Seasonal-Trend Decomposition Procedure Census I
This statistical method was applied to raw data to remove the periodical calendar related fluctuations (seasonal and cyclical) and the trend component, representing the non-casual long term

Analyses in Time Domain
With reference to the statistical theory [1][2][3][4][42][43][44][45][46], univariate autocorrelation and bivariate cross-correlation functions were applied to the time series as well as their residual component, obtained with the additive Seasonal-Trend decomposition procedure Census I. changes in the level of the time series. In this way, it was possible to extract from the original data the casual residual component describing the random short term fluctuations, neither systematic nor predictable [47][48][49].
The selected statistical approach is based on the assumption that the observed value x t at time t of the time series under consideration (or, in the case of a discrete ensemble of elements, the t-th data measured at the time t, thus 1 ≤ t ≤ N) can be thought as consisting of four different components, combined to each other in an essentially additive model: where T t is the trend component; S t is the seasonal component that occurs at regular (seasonal) intervals; C t is the cyclical component that varies from cycle to cycle and that have usually a longer duration than the seasonal component; and ε t is the random, error, or irregular component.
In the Census I method, the trend and cyclical components are customarily combined into a trend-cycle component T t × C t [47,48].

The Autocorrelation Function
The Autocorrelation Function (ACF) provides the opportunity to evaluate the linear dependency of successive values of a single parameter belonging to a time series. For this reason, the method is univariate and quantifies the "memory effect" that corresponds to the temporal reciprocal influence on subsequent data. Therefore, the ACF correlogram describes the time necessary to "forget" its own initial conditions [1,3,4,44,50,51].
In general, a correlogram with a gentle slope can be an indicator of data series persistence, while, if it drops off rapidly, it describes the random nature of the values [45]; furthermore, the slower decline of ACF function for hydrological data describes an aquifer with low draining properties, low permeability or major groundwater storage, conversely a faster decline of the autocorrelation function indicates a more rapid flow of water through the aquifer and/or its limited storage capacity [1,2,5,11,18,52,53].
To compare the "inertia" of signals, a delay (lag) is identified in correspondence of the decorrelation coefficient 0.2, since events with a lower ACF coefficient can be considered quasi-independent [1,6,54].

The Cross-Correlation Function
The Cross-Correlation Function (CCF), bivariate analysis performed between rainfall and river level input signals and between rainfall and piezometric output signals, describes the time lag necessary for the latter to reach the maximum CCF value in response of the input signals. Furthermore, the function's slope characterizes an aquifer in terms of infiltration rate, in terms of draining capacity and storage and, by doing so, it describes the system's modulating power and consequently its inertial nature [1][2][3]13,42].
The cross-correlation function, to be acceptable, has to be characterized by a significant correlation at the 95% confidence level, which means a correlation coefficient superior to the standard error ∼ 2 N = 0.037, where N is 2880 and corresponds to the number of values in the selected datasets [55,56]. Applying this function to the raw data, we obtain information on both the casual and non-casual components relationships of input and output variables, as well as their importance [3]. On the other hand, cross-correlation applied only to residuals describes the relationship between the random components of the data, giving an idea about the response of the system against a unitary impulse, hence the cross-correlogram can be used as the "unitary hydrograph" in surface hydrology [2,3]. Furthermore, the delay of the response gives an estimation of the pressure pulse transfer times (the piezometric level increase due to hydrostatic pressure increase in the aquifer) and of the particle travel times through the aquifer [5,6,9,44,[48][49][50]54,57]. In our study, because of the small surface area considered, the input rainfall signal, referred to SP rain gauge, is common to all wells and consequently the variability of the transfer function can be related only to site-specific structural and functional variability of the hydrological processes [57].

Analyses in Frequency Domain
The frequency domain is useful for data manipulation to extract frequency signals and provide a picture of the frequency contents that are unclear in the time domain data [58]. Furthermore, this analysis allows detecting average periodicities in long term series, to see if there are any dominants oscillation modes and to identify the frequency bands of interest [42,59,60].
Assuming that the time series under study {x t } is a linear combination of sinusoidal functions, each representing the multiple periodical nature of the phenomena, the discrete Fourier's analysis, also called harmonic analysis, can be employed. This univariate analysis decomposes the time series with cyclical components, into few underlying co-sinusoidal and sinusoidal functions of all possible (discrete) frequencies, each characterized by an amplitude related with the variance of the time series. Accordingly, when N is odd: where x is the mean of the time series, ν n = n/N is the n-th frequency (harmonics) multiple of the fundamental frequency 1/N; and ε t is the random component associated to the t-th data and q = (N − 1)/2. If N is even, q = N/2, a q = ∑ N t = 1 (−1) t · x t /N and b q = 0. It worth noting that the cosine and the sine parameters, respectively, a n and b n in Equation (2), are obtained applying the least squares estimation method and indicate the degree to which the respective cosine and sine functions are correlated with data. To explore how much each cyclic phenomenon affects the overall time series trend, the amplitude square of each component, function of the related frequency (spectral Analysis), was identified as an "importance measure" parameter.
The Fourier Transform, usually employed through its fast computational version (FFT; a particular gain is obtained if the number of data is equal to N = 2 n 1 × 3 n 2 × 5 n 3 × · · · , in this paper N = 2880 = 2 6 × 3 2 × 5, is a common approach aimed at performing univariate spectral analysis, providing the feasibility to switch from the temporal to the frequency domain. Accordingly, the following expression: , called "periodogram" or the "sample spectrum", was introduced. After some straightforward manipulations, a link between the sample spectrum and the estimate of the autocovariance, called also "power spectrum" [42,44] emerges. Moreover, in order to compare time series obtained at different scales of measurements, it is useful to introduce the estimated (for this reason characterized by a circumflex accent) "spectral density function" SDF(ν), normalizing the power spectrum respect to the varianceσ 2 x . Finally, the following expression was exploited: whereĉ xx (k) is the estimate of the k-th lag auto-covariance andĉ xx (0) ≡σ 2 x is the estimate of the variance (for more details see [42,44]). That is, SDF(ν) could be interpreted as the Fourier transform of the autocorrelation AĈF(k) of the time series. The physical meaning of the SDF(ν) is that the area under a portion of the total plot, related to a selected frequency range, represents the contribution to variance of components with frequencies within the range. When the entire spectrum is drawn, the equation indicates that the area underneath the curve is equal to the variance of the process. However, the simple Fourier analysis is based on the assumption of fixed amplitudes, frequencies and phases, while time series are affected by random changes of these kind of parameters.
A consequence is that Expression (3) fluctuates uncontrollably and, for this reason, is not useful for any meaningful interpretations. Accordingly, the SDF(ν) expressed by Equation (3) is not a consistent estimator [43] for an ensemble of variables taken from a discrete purely random process (de-trended time series). A way to solve this kind of issue is to smooth the periodogram, exploiting a sort of truncation of the autocorrelation function, giving less weight to the values of the AĈF(k) as k increases. As discussed in [42,43] , the following truncated SDF(ν) T appears to be a consistent estimator to consider properly the randomness nature of the time series under study: where {λ(k)} is an ensemble of weights called the "lag window" and M < N is called the "truncation point" [43]. In this paper, we adopted, among many others, the Blackman-Tukey window expression: The value of M, following the detailed discussion reported in [16,43], was selected equal to 5. Hereafter, to simplify the symbology, ACF, CCF and SDF stand for AĈF(k), CĈF(k), and SDF(ν) T .
Summarizing, the power spectrum of a time series is the simple translation of its autocovariance function (i.e., the ratio of autocorrelation and variance) in the frequency; complementary to correlation analysis, it may reveal certain characteristics present in a time series that are buried in noise and not discernable otherwise [61,62].

Fractal Analysis
The data clustering of time series, expressed by only the number of events occurring within a period, cannot fully describe the possible large time irregularities of events [24]. For this reason, a preliminary fractal approach application was proposed in this study, in addition to the time series statistical approach described above. The purpose of the fractal analysis was to gain more information in relation to the variability of the data, and therefore provide insights into this issue.
One of the most important properties of a fractal object, time series included, is that it can be subdivided into parts (or intervals of time), each of which is (usually only approximately) a reduced-size copy of the whole. This approximate property is termed "self-affine" and implies that properties such as frequencies of selected events are almost invariant across spatial and temporal scales. Their irregular or complex occurrences, in general, can be "measured" to some extent, by the distance of the Hausdorff-Besicovitch dimension D, called hereafter the "fractal dimension", from an integer number [19,20,[63][64][65]. However, not all time series data are truly fractal [66]. Therefore, a preliminary analysis of the assessed data is recommended in order to verify the occurrence of fractal behavior. A self-affine performance of a time series could be detected, for example, by long-term correlations resulting from the analysis of the ACF and by the evidence of a power-law relation between the power spectrum S(f) and the time frequency of the occurrence of a selected event: S( f ) ≈ f −β , resulting from the Spectral Analysis [66]. The exponent β is linked to the Hurst's coefficient H and to the fractal dimension D [67][68][69].
The fractal features of different hydrological time series, including rainfall and river dynamics, were highlighted in previous papers [19,20,22,23,63,[65][66][67][68][69][70][71][72]. To have a fast response, we carried out a fast fractal analyses as a first step. Nevertheless, as we demonstrated in the following sections, some interesting results were obtained. The feasibility of the fractal analysis application to the time series studied in this paper emerge from the inspection of the results obtained by both Autocorrelation Function and Spectral Analysis. Different types of fractal dimension and the methods of estimating these have been proposed in literature [65,73].
In this paper, an algorithm based on the Box-Counting and Cantor-Dust methods, successfully applied in [22,23,70,71], was developed. As a first attempt, only raw data were considered. Subsequently, the values of the data belonging to each time series, i.e. rainfall intensities, river hydrometry and three wells piezometric levels, were normalized with respect to the full range of their variability within their own time series. The resulting percentages (accordingly from 0 to 100%, to be considered as fluctuations above the minimum values) were considered as the "events" that can be compared to each other, although originating from different sources. In the following the related analytical expression is given: where f i,j (%) ("fluctuation") is the "event" related to the i-th data x i,j belonging to the j-th time series (j = 1,2,3,4,5, since the time series under consideration are 5); and x min,j and x max,j are, respectively, the minimum and the maximum numerical values of the j-th time series data. According to the algorithm procedures, different values of f i,j (%), hereafter "F_threshold", should be introduced.
No particular constraints are required, thus, in this paper 7-9 values between 0% up to 90% for both the river hydrometry and wells, while just only up to 40% for Spoltore rainfall were selected. Then the time range was divided into increasing integer numbers (1, 2, 3, . . . ) of equal-time interval groups, Ng. Each interval group contains the same integer number ε of measured data. Ng values that did not satisfy the Ng × ε = 2880 requirement was excluded. Therefore, the resulting Ng numbers are all the integer divisors of 2880. For each accepted number Ng of intervals, the number of interval time N(ε), in which at least one event above the selected F_threshold occurred, was counted. We introduced the parameter "n_events" as the number of the fluctuations f i,j (%) that, in the selected interval time, satisfy a selected comparison relationship (≤ or ≥) with a selected F_threshold. To be more clear, if over the interval time under consideration (i.e., from 9-th May 1990 to 10-th June 1990) the rainfall fluctuations f i,j (%), measured at Spoltore station, were equal to or above the selected F_threshold (i.e., 60%, thus F_threshold = 0.6) for 5 times, accordingly it results that n_events = 5 and the value of N(ε) is increased by just only 1. The following power law relation [22,23,65] was assumed: where m is a constant (equal to 1 in this paper) and D is a parameter that is identified as the fractal dimension if: ln N(ε) was plotted versus ln(1/ε) and D was estimated from the slope of the best fit straight line, obtained by means of the Minimum Square Root Method. Generally, the D parameter is not constant and depends on the selected F_threshold values. Accordingly, the value of the D parameter calculated in this paper should be intended, in general, as an average fractal D dimension. Moreover, the comparison of the fractal analyses to the results obtained in the previous statistical approaches should be considered as qualitative, rather than quantitative, although some quantitative considerations are also made.
To better clarify the physical meaning of the variability of the fractal dimension D by the mean of the Cantor Dust approach (D F in Figures 13 and 15), the following example can be useful. Two rainfall time series are analyzed. The data were obtained by daily measurements over 30 years (10,920 values), at two gauges, A and B, representative of two different sites (A and B). We are interested to study how the occurrences of daily rainfall above or equal to 50 mm (F_threshold) are distributed over the selected interval of time. If we consider a time scale of one year, the entire interval is subdivided into 30 subintervals, each lasting one year. Accordingly, N g = 30 and ε = 365, as resulting from the application of the Cantor Dust method. Suppose that for 40 times (n-events = 40) the intensity of the rainfall, measured at A and B stations exceeded the value of the selected F_threshold over, respectively, 2 and 22 different years, even non-continuously. Consequently, N(ε) = 2 and N(ε) = 22, respectively, for time series A and B.
It could be reasonably inferred that, within the site where Station A is located, rainfall values exceeding the considered F_threshold are rather uncommon events over the selected time scale of 30 years. Their occurrences are concentrated during just two years. On the other hand, within the site where Station B is located, this kind of events is quite common and the expected return time would be reasonable equal to about one year. Then, consider a further partition of the whole time range into 60 subintervals (N g = 60 and ε = 182). Moreover, assume that, as the time scale decreases, N(ε) = 2 remains the same for Station A, while N(ε) = 37, increases (thus almost 1 for each semester) for Station B. Repeating this operation by the means of the Cantor Dust approach, the resulting value of the D fractal dimension related to the time series A would be low (in this case close to 0), while for the other time series, D would be close to 1.
As highlighted by the previous example and according to Mazzarella [23,24], we can infer that, for the time series considered in this paper, the D Cantor Dust fractal dimension can range from 0, corresponding to randomly distributed events, just as isolated points to the lim ε→0 D = 0, to 1 (lim ε→0 D = 1: a line by a geometrical point of view) corresponding to a distribution of the selected events occurrence that remains almost the same within all the considered time scales.

Descriptive Statistics
The variability of the data, described by the standard deviation analysis (Table 2), points out that the piezometric level fluctuations within the Su and Sa wells are higher than in the DN well. Their average water depth (unsaturated zone thickness) is calculated from the difference between the altitude and the average piezometric level for the 1986-2009 period. Averaging for each month the 10 available piezometry, hydrometry and cumulated rainfall data, over the 24-year period of 1986-2009, it is possible to represent the monthly average of each variable ( Figure 6). River level at the PST site presents the maximum value in April and the minimum value in July, while the Sp rain gauge presents a maximum in December and a minimum in July. Furthermore, the standardization of the computed piezometric levels through the normalization, based on the mean and standard deviation, gives us the chance to compare at the same scale the entire 1986-2009 group of data, showing that all wells present the lowest piezometric level in August and the highest in January, February and March, respectively, for the DN, Su and Sa wells. Averaging for each month the 10 available piezometry, hydrometry and cumulated rainfall data, over the 24-year period of 1986-2009, it is possible to represent the monthly average of each variable ( Figure 6). River level at the PST site presents the maximum value in April and the minimum value in July, while the Sp rain gauge presents a maximum in December and a minimum in July. Furthermore, the standardization of the computed piezometric levels through the normalization, based on the mean and standard deviation, gives us the chance to compare at the same scale the entire 1986-2009 group of data, showing that all wells present the lowest piezometric level in August and the highest in January, February and March, respectively, for the DN, Su and Sa wells.

The Residual Analysis
The comparison at a three-day resolution between the 2880 piezometric residuals data for the whole 1986-2009 period identifies in Su site the records with the greatest piezometric variability; furthermore, a similar residual distribution especially for the Su well, the Sp rain gauge and the PST river level can be identified, since the peaks fall together on several dates (Figure 7, left side).
The periods during which the casual components are concurrently more influent can be identified by calculating the average annual series on the basis of a monthly arrangement of the data, which consisted in averaging the 10 monthly three-day resolution records, over the 24 years of 1986-2009 period (Figure 7, right side); it is evident that the months of March-April, October-November and June-July, in correspondence of the spring, autumn and summer periods, respectively, present the greatest random variability.
Some exceptions can be highlighted: for example, a greater residual concentration during the 2007 summer is self-evident at the DN location and is characterized by a strong variability of the raw data. In December 1994, an intense rainfall event was observed at the Spoltore rain gauge and a correspondent peak was identified in the rainfall residual. However, it was not identified in the datasets of river stage and piezometric head residuals, probably because of the local character of this event.

The Residual Analysis
The comparison at a three-day resolution between the 2880 piezometric residuals data for the whole 1986-2009 period identifies in Su site the records with the greatest piezometric variability; furthermore, a similar residual distribution especially for the Su well, the Sp rain gauge and the PST river level can be identified, since the peaks fall together on several dates (Figure 7, left side).
The periods during which the casual components are concurrently more influent can be identified by calculating the average annual series on the basis of a monthly arrangement of the data, which consisted in averaging the 10 monthly three-day resolution records, over the 24 years of 1986-2009 period (Figure 7, right side); it is evident that the months of March-April, October-November and June-July, in correspondence of the spring, autumn and summer periods, respectively, present the greatest random variability.
Some exceptions can be highlighted: for example, a greater residual concentration during the 2007 summer is self-evident at the DN location and is characterized by a strong variability of the raw data. In December 1994, an intense rainfall event was observed at the Spoltore rain gauge and a correspondent peak was identified in the rainfall residual. However, it was not identified in the datasets of river stage and piezometric head residuals, probably because of the local character of this event.

Autocorrelation
The correlogram relative to rainfall observed at the Sp rain gauge ( Figure 8A) shows how, after a three-day time lag, the ACF function immediately decreases to uncorrelated values. This confirms its random behavior because events with duration longer than three days are rare and, therefore, the correlation between the data for intervals longer than the reference lag is nil.

Autocorrelation
The correlogram relative to rainfall observed at the Sp rain gauge ( Figure 8A) shows how, after a three-day time lag, the ACF function immediately decreases to uncorrelated values. This confirms its random behavior because events with duration longer than three days are rare and, therefore, the correlation between the data for intervals longer than the reference lag is nil. The river level reaches the decorrelation value threshold ACF = 0.2 after 81 days, highlighting a strong memory effect ( Figure 8B). However, the two factors influencing the ACF function shape are the runoff velocity and the catchment area; indeed, the curve shape in the first three days, when the ACF value decreases to 0.7, is characterized by a greater slope and a more rapid decorrelation rate down, meaning a stronger independence of the data. In fact, in this short period, the river level ACF curve has the same shape of the rainfall one because it describes the flooding action caused by direct runoff, while after three days the more autocorrelated shape of the river level ACF is controlled by the baseflow due to the large recharge area and groundwater contribution [74].
The hydrogeological time series as the piezometric levels of our study are strongly autocorrelated in all sites because of inertia and carryover processes in the physical system, that impart a correlation to successive piezometry levels. In detail, Su reaches 0.2 threshold value after 880 days, while the ACF curves for Sa and DN reach that level after respectively 105 and 90 days ( Figure 8). Su well stronger autocorrelation might be related to site-specific physical and structural conditions of the aquifer, such as the hydraulic conductivity influenced by changes in the grain size of the alluvial aquifer or its degree of compaction [18,75], because, of course, as suggested by Imagawa et al. [53], persistence can be an indicator of the local storage capacity of alluvial aquifers, provided that the assumption of homogeneity in the hydraulic conductivity throughout the aquifer is respected.

Raw Data Cross-Correlation
The piezometric level fluctuation in answer to rainfall (Figure 9) highlights low correlations value for all the wells. The Su and DN wells behave the same, since the shape of the curves are similar, with a maximum CCF value around 0.16, reached at both sites after, respectively, a time lag of nine and six days. The river level reaches the decorrelation value threshold ACF = 0.2 after 81 days, highlighting a strong memory effect ( Figure 8B). However, the two factors influencing the ACF function shape are the runoff velocity and the catchment area; indeed, the curve shape in the first three days, when the ACF value decreases to 0.7, is characterized by a greater slope and a more rapid decorrelation rate down, meaning a stronger independence of the data. In fact, in this short period, the river level ACF curve has the same shape of the rainfall one because it describes the flooding action caused by direct runoff, while after three days the more autocorrelated shape of the river level ACF is controlled by the baseflow due to the large recharge area and groundwater contribution [74].
The hydrogeological time series as the piezometric levels of our study are strongly autocorrelated in all sites because of inertia and carryover processes in the physical system, that impart a correlation to successive piezometry levels. In detail, Su reaches 0.2 threshold value after 880 days, while the ACF curves for Sa and DN reach that level after respectively 105 and 90 days (Figure 8). Su well stronger autocorrelation might be related to site-specific physical and structural conditions of the aquifer, such as the hydraulic conductivity influenced by changes in the grain size of the alluvial aquifer or its degree of compaction [18,75], because, of course, as suggested by Imagawa et al. [53], persistence can be an indicator of the local storage capacity of alluvial aquifers, provided that the assumption of homogeneity in the hydraulic conductivity throughout the aquifer is respected.

Raw Data Cross-Correlation
The piezometric level fluctuation in answer to rainfall (Figure 9) highlights low correlations value for all the wells. The Su and DN wells behave the same, since the shape of the curves are similar, with a maximum CCF value around 0.16, reached at both sites after, respectively, a time lag of nine and six days. The cross-correlation analysis between groundwater levels and river level highlights a major dependency of groundwater level fluctuations to the river stage, rather than to rainfall. In fact, the maximum CCF values vary from 0.20 in Su, 0.25 in DN and 0.37 in Sa with corresponding time lags of three days (Su, DN) and six days in Sa (Figure 9). The DN well, which is 1960 m from Pescara River, presents the intermediate maximum cross-correlation value.

Residual Cross-Correlation
The CCF between rainfall and river level residual components shows an initial high positive correlation value (0.34), indicating that an immediate positive fluctuation of river stage corresponds to the three-day precipitation impulse. After three days, this correlation becomes negative, highlighting that river level continues to rise positively for three days after the end of the rainfall event ( Figure 10). It is important to remember that the river level is affected by both the local precipitation events (monitored at the Sp rain gauge) and by those events happened upward in the river drainage basin, far from the selected rain gauge. The rainfall-piezometric and river-piezometric levels residual cross-correlation are similarly shaped, even though the river impulse affects the well fluctuations more than rainfall, especially for the Su and Sa wells, which are closer to the river. In fact, in Su well, the residuals CCF precipitation-piezometric level and river level-piezometric level are, respectively, 0.28 and 0.36 ( Figure 11); these maximum values are reached immediately as a consequence of the first impulse The cross-correlation analysis between groundwater levels and river level highlights a major dependency of groundwater level fluctuations to the river stage, rather than to rainfall. In fact, the maximum CCF values vary from 0.20 in Su, 0.25 in DN and 0.37 in Sa with corresponding time lags of three days (Su, DN) and six days in Sa (Figure 9). The DN well, which is 1960 m from Pescara River, presents the intermediate maximum cross-correlation value.

Residual Cross-Correlation
The CCF between rainfall and river level residual components shows an initial high positive correlation value (0.34), indicating that an immediate positive fluctuation of river stage corresponds to the three-day precipitation impulse. After three days, this correlation becomes negative, highlighting that river level continues to rise positively for three days after the end of the rainfall event ( Figure 10). It is important to remember that the river level is affected by both the local precipitation events (monitored at the Sp rain gauge) and by those events happened upward in the river drainage basin, far from the selected rain gauge. The cross-correlation analysis between groundwater levels and river level highlights a major dependency of groundwater level fluctuations to the river stage, rather than to rainfall. In fact, the maximum CCF values vary from 0.20 in Su, 0.25 in DN and 0.37 in Sa with corresponding time lags of three days (Su, DN) and six days in Sa (Figure 9). The DN well, which is 1960 m from Pescara River, presents the intermediate maximum cross-correlation value.

Residual Cross-Correlation
The CCF between rainfall and river level residual components shows an initial high positive correlation value (0.34), indicating that an immediate positive fluctuation of river stage corresponds to the three-day precipitation impulse. After three days, this correlation becomes negative, highlighting that river level continues to rise positively for three days after the end of the rainfall event ( Figure 10). It is important to remember that the river level is affected by both the local precipitation events (monitored at the Sp rain gauge) and by those events happened upward in the river drainage basin, far from the selected rain gauge. The rainfall-piezometric and river-piezometric levels residual cross-correlation are similarly shaped, even though the river impulse affects the well fluctuations more than rainfall, especially for the Su and Sa wells, which are closer to the river. In fact, in Su well, the residuals CCF precipitation-piezometric level and river level-piezometric level are, respectively, 0.28 and 0.36 ( Figure 11); these maximum values are reached immediately as a consequence of the first impulse The rainfall-piezometric and river-piezometric levels residual cross-correlation are similarly shaped, even though the river impulse affects the well fluctuations more than rainfall, especially for the Su and Sa wells, which are closer to the river. In fact, in Su well, the residuals CCF precipitation-piezometric level and river level-piezometric level are, respectively, 0.28 and 0.36 ( Figure 11); these maximum values are reached immediately as a consequence of the first impulse event, while after six days the correlation becomes negative (Figure 11) pointing out that the piezometric level continues to rise for six days after the rainfall-river stage pulse.
Water 2017, 9,850 16 of 28 event, while after six days the correlation becomes negative ( Figure 11) pointing out that the piezometric level continues to rise for six days after the rainfall-river stage pulse. At Sa well, the maximum CCF value is lower (0.05-0.1), while DN well level presents slight correlations to rainfall and river stage (0.06-0.05).

Spectral Analysis
The power spectrum of precipitation, containing homogeneously unimodal distributed amplitudes, is dominated only by an annual cycle representing the wet/dry seasonal cycle and is void of other prevailing periodic component. Hydrometric and piezometric levels spectrum, instead, are not homogeneously distributed, and present a bimodal behavior, because in addition to the seasonal cycle present high amplitude oscillations at lower-frequencies (higher periods) concentrated particularly between five and eight years, with broader peaks more evident in correspondence of the interannual 72-month cycle (≈6 years) and of the interdecadal 144-month cycle (≈12 years), especially in Su and Sa wells located closer to the river ( Figure 12).  At Sa well, the maximum CCF value is lower (0.05-0.1), while DN well level presents slight correlations to rainfall and river stage (0.06-0.05).

Spectral Analysis
The power spectrum of precipitation, containing homogeneously unimodal distributed amplitudes, is dominated only by an annual cycle representing the wet/dry seasonal cycle and is void of other prevailing periodic component. Hydrometric and piezometric levels spectrum, instead, are not homogeneously distributed, and present a bimodal behavior, because in addition to the seasonal cycle present high amplitude oscillations at lower-frequencies (higher periods) concentrated particularly between five and eight years, with broader peaks more evident in correspondence of the interannual 72-month cycle (≈6 years) and of the interdecadal 144-month cycle (≈12 years), especially in Su and Sa wells located closer to the river ( Figure 12).
Water 2017, 9,850 16 of 28 event, while after six days the correlation becomes negative ( Figure 11) pointing out that the piezometric level continues to rise for six days after the rainfall-river stage pulse. At Sa well, the maximum CCF value is lower (0.05-0.1), while DN well level presents slight correlations to rainfall and river stage (0.06-0.05).

Spectral Analysis
The power spectrum of precipitation, containing homogeneously unimodal distributed amplitudes, is dominated only by an annual cycle representing the wet/dry seasonal cycle and is void of other prevailing periodic component. Hydrometric and piezometric levels spectrum, instead, are not homogeneously distributed, and present a bimodal behavior, because in addition to the seasonal cycle present high amplitude oscillations at lower-frequencies (higher periods) concentrated particularly between five and eight years, with broader peaks more evident in correspondence of the interannual 72-month cycle (≈6 years) and of the interdecadal 144-month cycle (≈12 years), especially in Su and Sa wells located closer to the river ( Figure 12).  Hydrometric and piezometric level response to the rainfall main driver results from the passage of infiltrated water through several subsurface filters that, depending on the scale range, have the effect to either amplify or attenuate the precipitation variability modes [76,77]. Accordingly, the hydrometric and piezometric level output signals can have significantly different spectral characteristics than the input precipitation signal; in our study comparing their power spectrum in a log-log scale graph it is evident how strongly is related the frequency of piezometry and hydrometry, while this relationship is much weaker with precipitation frequency. Consequently, the unimodal spectrum of the latter is not affected by a power law type, whilst the previous are characterized by a bimodal power spectra shifted towards low frequency scales (Figure 12, in red).
Spectral density and amplitude in Figure 12 indicate that piezometric trends seem to follow approximately a decaying exponential function S( f ) ≈ f −β = T β (T is the period and f is the frequency). They are localized in two intervals for Su and Sa, 12 months < T < (72~144) months and 0 months < T < 12 months, and only in the latter interval for DN. These occurrences further support the applicability of fractal analysis, as discussed above. Moreover, according to many authors [78][79][80] due to the variation of the trends observed in Figure 12 (related to the piezometry Spectral densities, before and after one year), it could be argued that the time series of the piezometric levels and just weakly for the hydrometry, show a multifractal behaviors.

Cantor-Dust Fractal Analysis
The characteristic of the selected F_threshold events and their occurrences are detailed in Table 3 for each time series. More details regarding the meaning of parameters are discussed in the text, in particular in the Discussion section (Section 6.3. Fractal analysis).
In Figure 13, for Sp rainfall and PST hydrometry respectively, the highest values of the average fractal dimension, ranging from D = 1.000 to D = 0.851 and D = 0.906 were reached for level values up to 5% and 10%, corresponding to ∆L = 100 mm and ∆L = 0.52 m fluctuations ( Table 3).
The outcomes of piezometric level variations are plotted in Figure 14. Irregular distributions of points occur for all the three wells at high values of F_threshold events with normalized level variations higher than 80-90%. This kind of drawback, due to the sensitivity of the method to small differences in the box size, can be solved using the Shifted Box Counting Method [22]. For the DN time series (Table 3), normalized events, related to values ranging from 1% to 20% (Figure 14 (Table 3). a fractal dimension that is very close to unity. As discussed above, if the fractal dimension of a dataset is close to unity for a specific time scale, the data show similar distributions within time scales higher than the identified one. Accordingly, this time-scale could be considered for the introduction of a "distribution return time" (DRT) of the selected event, which should not be confused with the commonly used "event return time", because the last parameter is calculated by means of different approaches. values are reported in Figure 15.  It is roughly possible to identify several "breaking points" in Figures 13 and 14 and the underlying data on which the figures are based. Such breaking time scales are specific for each F_threshold value and for each time series source, for which the relation between the number of covering boxes and the interval time (box size), shows a departure from a linear trend. The latter has a fractal dimension that is very close to unity. As discussed above, if the fractal dimension of a dataset is close to unity for a specific time scale, the data show similar distributions within time scales higher than the identified one. Accordingly, this time-scale could be considered for the introduction of a "distribution return time" (DRT) of the selected event, which should not be confused with the commonly used "event return time", because the last parameter is calculated by means of different approaches. Therefore, DRT values (estimated through inspection of Figures 13 and 14 and the underlying data) of different f i,j (%) values are reported in Figure 15.   It is interesting to note that the Sp Rainfall and the DN estimated distribution curve of DRT show a similar sigmoidal regression trend of the following type: where T is the DRT (days); Nvp is the Normalized variability percentage (that is f i,j (%)); and "a", "b", and "c" are three parameters, the values of which are different rainfall and DN data. The regression of DRT related to the other sources follows exponential trends: where "a", "b", and "c" are the three regression parameters.

Discussion
The relationship between rainfall, river level and piezometric levels was assessed by means of statistical-mathematical analyses and a conceptual model of the local hydrogeology.

Analyses in Time Domain
The DN well is the furthest from the main stream, therefore recharge effects from the Pescara River are considered to be negligible, while rainfall infiltration represents the main local recharge; this hypothesis is confirmed by the CCF analysis of the rainfall and piezometry data ( Figure 9A), which shows that DN reaches the highest CCF value. An additional confirmation is provided by the CCF of the residuals, which indicates the absence of a relationship between the DN well and the river ( Figure 11B) meaning that extreme rainfall events appear to have less influence on the piezometry, as suggested by its low residual component. The higher CCF value between the well and the river level ( Figure 9B) can be explained by the simultaneous changes of piezometric and river levels, due to the seasonal cycle, even if the two oscillations are not physically correlated. This is confirmed by the hydrogeology of the area. In fact, the water table depth in the DN well is quite low (Figure 4; Table 2). Similarly, the difference between the head in the DN well (mean 10.46 m a.s.l.) and the elevation of the Pescara River (mean 3.33 m a.s.l.) is significant. In addition, the groundwater hydrodynamic indicates the presence of a groundwater divide, which has a direction that is approximately parallel to the direction of the river (Figure 1). This excludes that the water level changes in the well can be induced by the river.
Conversely, groundwater and surface water interaction is likely in the area near the Su and Sa wells, since they are located near the Pescara River ( Figure 16). this hypothesis is confirmed by the CCF analysis of the rainfall and piezometry data ( Figure 9A), which shows that DN reaches the highest CCF value. An additional confirmation is provided by the CCF of the residuals, which indicates the absence of a relationship between the DN well and the river ( Figure 11B) meaning that extreme rainfall events appear to have less influence on the piezometry, as suggested by its low residual component. The higher CCF value between the well and the river level ( Figure 9B) can be explained by the simultaneous changes of piezometric and river levels, due to the seasonal cycle, even if the two oscillations are not physically correlated. This is confirmed by the hydrogeology of the area. In fact, the water table depth in the DN well is quite low (Figure 4; Table 2). Similarly, the difference between the head in the DN well (mean 10.46 m a.s.l.) and the elevation of the Pescara River (mean 3.33 m a.s.l.) is significant. In addition, the groundwater hydrodynamic indicates the presence of a groundwater divide, which has a direction that is approximately parallel to the direction of the river (Figure 1). This excludes that the water level changes in the well can be induced by the river.
Conversely, groundwater and surface water interaction is likely in the area near the Su and Sa wells, since they are located near the Pescara River ( Figure 16).  The Su well is located near the river (Figure 4, Table 1) and has piezometric heads that are comparable to the elevations of the river section that falls just upstream of the well itself. The piezometric level spatial distribution indicates that groundwater flows from the river to the aquifer, thus recharging the latter ( Figure 16). The analysis of the raw data CCF ( Figure 9A) indicates that there is a low correlation between rainfall and piezometric level fluctuations. This is in agreement with both the high average depth of the water table (6.70 m; Table 2) and with the idea of aquifer recharge from the river. The value of the river/well raw data CCF ( Figure 9B) is generally lower compared to the corresponding residual component CCF ( Figure 11B). A possible explanation is the position of the gauging Station compared to the well location; the analysis of the hydraulic conditions ( Figure 3) indicates that the river section in correspondence of the gauging Station is wider compared to the section just upstream of the well. Therefore, considering the same discharge, a wider section should imply smaller river level changes. The higher cross-correlation of the residuals confirms this. In fact, when high discharges occur during flood events, the higher river level recorded downstream of the well is indicative of even higher changes upstream where the river section is narrower. However, the strong impulsive changes in the piezometric level within the well, clearly represented by the strong oscillation of the residual component ( Figure 11B), is due to a pressure transfer rather than a mass transfer, given the low hydraulic conductivity of the aquifer (about k = 10-5 m/s).
The Sa well is also located near the Pescara River ( Figure 4) and has piezometric head that is constantly higher than the river level upstream of the well itself ( Figure 16). This indicates a constant discharge from the aquifer to the river. The piezometric level supports this idea and the flow lines are constantly towards the river. The cross-correlation with the rainfall shows a very low value ( Figure 9A) that, together with the high water depth (mean 3.89 m; Table 2), allows assessing the low contribution of the direct recharge from local rainfall. The CCF with the river level ( Figure 9B) shows a relatively high value. This apparent anomaly can be explained with the presence of a significant annual cycle in both datasets, based on the hydrogeological understanding of the area. A further confirmation is provided by the CCF function of the residuals, which shows very low values for both the rainfall and the river level fluctuation (Figure 11).

Analyses in Frequency Domain
The spectral analysis has highlighted a clear seasonal cycle (one year) in piezometry, hydrometry and rainfall spectral densities distribution. Furthermore, evident multi-year periodicities have been detected in the piezometric and river levels fluctuations (pointed out in Figure 12). The rainfall spectrum is characterized by a unimodal and homogeneous power spectrum distribution and, even though multi-year cycles are present, they are not dominant compared to the annual cycle. On the other hand, piezometric and river levels present a clearly bimodal and not homogeneous spectral density distribution, characterized by a sharp spectral power decrease with frequency increase (Figure 12).
This effect is well known in the literature and is due to low-frequency pass filter effect of the aquifer and of the river basin, on the piezometric and river levels, respectively [53,59,76,78,81]. Thus, these two measures modulate the rainfall input signal, lowering the short period fluctuations and highlighting the long period ones (seasonal and multi-year) that are in association to large scale climatic oscillations, as many researchers have reported [59,79,[82][83][84][85][86][87].
The water volumes that are at stake proportionally influence the weight of this filtering effect [16,[75][76][77], therefore extended aquifers and river basins can be considered as excellent climate signal proxies, at regional scale [18,88,89]. Similar results are described by Luque-Espinar [16] and Fendeková et al. [90], who found that piezometers, located close to stream with leakage to the aquifer, or along the main river in those areas where it recharges the aquifer, present long period cycles. In addition, Chiaudani et al. [91] found out the presence of multi-year and decadal cycles in Abruzzo Region in 11 wells belonging to alluvial aquifers located near to their main stream.

Fractal Analysis
It was possible to estimate a sort of "pseudo" fractal dimension through the fractal analysis discussed in this paper. This fractal dimension can be regarded as average values of trends, within the appropriate log-log domain for this kind of investigation. As expected, the rainfall plot (left of Figure 13) shows that these types of events are more "isolated" and random than "events" associated to river level (right of Figure 13) and piezometric levels ( Figure 14). This can be explained by the fact that the river level and piezometric variations are more continuous in time. They are the result of a much more complex phenomena related to not only rainfall at one site, but also to the hydrological dynamics of the entire basin.
The river level plot shows that almost the same event distribution is related to variations of levels lower than 10% (<0.52 m). This distribution of events, which corresponds to F_threshold = 3.26 m level (see Table 3) occurs within all time scales. It is interesting to note that the qualitative trends, shown in Figure 14, are in agreement with the trends obtained for the Warta River flows in [22], even if the phenomena dynamics are different. The occurrence of the same fractal dimension equal to 1 means that each distribution of such level fluctuations remains similar in all the selected time interval scales; however, distributions related to different sources of data are not necessary equal. Consequently, it would be hard to identify any correlation between different time series, where only the value of the fractal dimension obtained by the Cantor Dust approach is observed. Figure 15 and Table 3 show that the DRT distributions of both rainfall and the DN piezometric data are well interpolated by sigmoidal curves, while Su, Sa and river stage are better interpolated by exponential ones. This result is in agreement with the statistical analysis that identifies two distinct groups of data with respect to well defined hydrological relationships. The shift of the sigmoidal curve of DN well respect to the rainfall one indicates that the piezometric levels shows a random behavior for higher L_thresholds, as the breaking points in Figures 13 and 15 and the previous discussion highlight. This could be related to the modulator effect of the aquifer which gives higher stationarity to DN time series. Moreover, the interpolating sigmoidal and exponential expressions could be compared to the classical return-time laws proposed in literature ( [92], and references therein).

Conceptual Model
All the applied statistical-mathematical analyses, considering hydrogeological knowledge, allowed defining a final conceptual model of the system hydrodynamic ( Figure 17).

Conclusions
The assessment of the relationships between groundwater and both rainfall and river recharge influencing the alluvial aquifer hydrodynamic needs a detailed hydrogeological framework definition. Long time series are useful for this purpose, and, in this study, hydrological data series, monitored for the 24-year period 1986-2009, have been analyzed.
The statistical-mathematical analysis of rainfall, river level and groundwater level time series data has produced useful information for the understanding of their relationships existing within the studied hydrogeological context. To better understand the hydraulic conditions and to make statistical results coherent, aerial photo-interpretation analyses and in-situ GPS surveys have been realized.
The "memory effect", representing the self-coherency indicator of each time series, was assessed with the Auto-Correlation Function and with the fractal analysis pointing out how stationary behaviors were higher for piezometric and river levels than for rainfall. Interdependency between different hydrological parameters was found by means of the Cross-Correlation Function, highlighting a strong groundwater/surface water interaction monitored between the Pescara River and the nearest wells. On the other hand, the farthest well showed a clear correlation with local rainfall. This occurrence could also be observed by fractal analysis from the regression trends of the "distribution return times", for high time scales.
Moreover, the Cross-Correlation Function applied to residual data, after seasonal cycle and trend removal, points out a strong pressure transfer from the river to groundwater during flooding In particular, the DN well is located in a Pescara aquifer part where a secondary flow path has been recognized. Here, the recharge is completely related with local rainfall, described by the Sp rain gauge, and with the runoff infiltration at the clayey slope base.
On the other hand, the main flow path has been identified near the Pescara River. In this area, there are clear surface-water/groundwater relationships, considering either an aquifer recharge from the main stream or the opposite. In detail, upstream, where the Su well is located, the main river recharges the aquifer. Furthermore, when huge flooding events occur the Su well detects outright the increase of the hydrometric level, probably due to a pressure transfer. Contrariwise, downstream near the Sa well, the piezometric head is constantly higher than river level, even during huge flooding events. Thus, it is clear that the aquifer discharges groundwater into the Pescara River.
Finally, the local rainfall monitored in Sp rain gauge, near the Pescara River, does not seem to influence substantially the piezometric fluctuation, due to the importance of aquifer baseflow.

Conclusions
The assessment of the relationships between groundwater and both rainfall and river recharge influencing the alluvial aquifer hydrodynamic needs a detailed hydrogeological framework definition. Long time series are useful for this purpose, and, in this study, hydrological data series, monitored for the 24-year period 1986-2009, have been analyzed.
The statistical-mathematical analysis of rainfall, river level and groundwater level time series data has produced useful information for the understanding of their relationships existing within the studied hydrogeological context. To better understand the hydraulic conditions and to make statistical results coherent, aerial photo-interpretation analyses and in-situ GPS surveys have been realized.
The "memory effect", representing the self-coherency indicator of each time series, was assessed with the Auto-Correlation Function and with the fractal analysis pointing out how stationary behaviors were higher for piezometric and river levels than for rainfall. Interdependency between different hydrological parameters was found by means of the Cross-Correlation Function, highlighting a strong groundwater/surface water interaction monitored between the Pescara River and the nearest wells. On the other hand, the farthest well showed a clear correlation with local rainfall. This occurrence could also be observed by fractal analysis from the regression trends of the "distribution return times", for high time scales.
Moreover, the Cross-Correlation Function applied to residual data, after seasonal cycle and trend removal, points out a strong pressure transfer from the river to groundwater during flooding events, when the river leakage recharges the aquifer.
The spectral analysis identified in all data series a predominant annual cycle, linked to seasonal fluctuations, while multi-year climatic related regimen are detected in the hydrometry and in piezometry affected respectively from the low-pass filter action of the river basin and the aquifer. These results are supported also from autocorrelation stationary analysis and the fractal approach that shows a more systematic behavior (high D values) both in hydrometry and in piezometry than in rainfall (low D values).
These research results highlighted that statistical-mathematical analyses performed on long time series, when coupled with detailed assessments of the local hydrogeological conditions, can provide a quantitative understanding of the spatial-temporal relationships that exist between input and output parameters. In this respect, the application of such tools can be also useful for numerical modeling and in all those cases (e.g., water resource management, groundwater vulnerability, etc.), where it is important to understand the type and the entity of the relationships between surface-water and groundwater.