Application of short time series analysis for the hydrodynamic characterization of a coastal karst aquifer: the Salento aquifer (Southern Italy)

Water scarcity together with its ever going deterioration is a crucial issue experienced across the world, and its solution is a top priority for an increasing number of countries. Especially in the Mediterranean area, groundwater in coastal karst aquifers is essential to supply fresh water for drinking and irrigation purposes. Such resources, being vulnerable and highly exploited for many decades, are subject to pollution and saltwater intrusion to a variable space and time extent (Leduc et al. 2016). Because of the complex nature of coastal aquifers, their management and projection of mitigation measures should be based on a clear understanding of their hydrodynamic behaviour. A coastal aquifer is complex in itself from the hydrogeological point of view, and this complexity can only increase if it is of karst type. The organized heterogeneity (Kiraly 1998), which is the distinctive characteristic of a karst aquifer, refers to the spatial distribution of intrinsic permeability (and, therefore, of hydraulic conductivity) and reflects the degree of evolution of karst processes within the carbonate masses.

In general, the intrinsic permeability of a karst aquifer is described following a ‘triple porosity’ model, where the overall permeability is composed of various essentially independent proportions, of matrix, fractures, and conduits. However, there are karst aquifers in which an element can be dominant. The hydraulic fields connected to each of these elements of intrinsic permeability dominate at different scales (Jeannin 1998).

The first element refers to the rock matrix is characterized by low specific discharge and hydraulic conductivity. The second element is the system of small and medium fractures in the rock matrix. The third relates to the major conduits and fractures.

In general, while the transmissivity develops in the network of conduits, the role of storage capacity is played by the associated system of fractures and the less permeable rock surrounding the network itself. The conduit component is the most essential feature of a karst system and is the principal source of uncertainty in its study.

Ultimately, a karst aquifer is among the types of aquifers where the water flow shapes the aquifer voids; this induces specific characteristics to the aquifer itself, which condition, according to a feedback loop, water circulation. Moreover, during high-intensity precipitation, the hydraulic head in the conduits is larger than in the rock matrix and fractures (Mangiarotti et al. 2019): the flow is either driven from the conduits to the matrix or directly discharges from springs. During low or no precipitation periods, the matrix slowly releases water into conduits producing long recessions in spring hydrographs. Technical literature reports several approaches and methods to study the water exchanges between conduits and surrounding matrix (Martin & Dean 2001; Peterson & Wicks 2005; Bailly-Comte et al. 2009; Li et al. 2020); the hydraulic relationships between the conduit networks and set of fissures and fractures are a factor of complexity (Jenkins & Watts 1968) influenced by the geometry of conduits, fractures, and rock matrix, and by water pressure (Li et al. 2020). The above leads to the conclusion that the karst systems do not allow for a simple description. As a result, it is necessary to decipher the intrinsic characteristics of the aquifer indirectly, studying the effects of permeability structure on the hydraulic response of the system, rather than directly assessing its values and spatial and temporal distribution.

Because of aquifer anisotropy and heterogeneity, such schemes consider karst as a filter that transforms the input signal represented by the precipitation into an output signal (Larocque et al. 1998), as the spring discharge (Benavente & Pulido-Bosch 1985; Angelini 1997; Larocque et al. 1998; Amraoui et al. 2003; Fiorillo & Doglioni 2010; An et al. 2019), the groundwater level (Imagawa et al. 2013; Delbart et al. 2014; Cai & Ofterdinger 2016; Chiaudani et al. 2017), the river flow rate (Bailly-Comte et al. 2008; Chiaudani et al. 2017) or other physical or chemical parameters characterizing the aquifer (Bailly-Comte et al. 2011).

The uncertainties about the karst structure become even higher when considering extensive coastal karst aquifers, because of the spatial variability of permeability and fluid density variation.

The Salento aquifer is a coastal karst aquifer in Puglia Southern Italy: saltwater represents the lower hydraulic limit of the regional groundwater flow system and groundwater outflow only occurs to the sea through diffuse flow and focused brackish karst coastal springs (Fidelibus & Pulido-Bosch 2019). The Salento karst terrain does not favour surface water accumulation. As there is also a lack of freshwater springs and extra-regional water resources, Salento relies on groundwater as the only source of supply for drinking, irrigation, and industrial purposes. Groundwater over-exploitation started in the 1960s, consequent to the demographical and economic development, thus causing a relentless groundwater quantity reduction and quality deterioration (Cotecchia et al. 2005; Polemio 2005; De Filippis et al. 2019; Alfio et al. 2020; Di Nunno & Granata 2020). Because of the current exploitation regime and the projected climate change induced impacts, knowledge of the hydrodynamic response of the aquifer to precipitation is a crucial aspect for present and future water management. Water balance of the Salento aquifer is difficult to estimate, especially because of the uncertainties in the evaluation of the outflows, comprised of the discharge into the sea (diffuse flow and focused springs), and the exploitation, considering that thousands of not licensed wells abstract large volumes of water. Considering the hydrogeological complexity as a whole, we selected the time series analysis, which is an indirect method, as an effective approach to qualify the intrinsic permeability structure that drives the recharge–discharge processes of the Salento coastal karst aquifer.

This study presents the results of a short time series analysis based on four-, three-, and two-year time series of precipitation and groundwater level (GWL), the latter obtained by pressure transducers with hourly monitoring frequency in 11 monitoring wells. Section 2 describes the study area and the Salento aquifer; Section 3 presents the considered dataset, and the methods; Section 4 summarizes the results, discussing the key aspects of the short time series analysis; Section 5 shares thoughts on the adopted methods and discusses the corresponding hydrogeological implications. Finally, Section 6 outlines the major conclusions of the study.

The study area is located in the south-eastern part of Puglia (Salento Peninsula, Southern Italy), which extends from the Ionian to the Adriatic Sea (Figure 1). Its borders roughly coincide with those of the Lecce province, covering 2,760 km². The Salento Peninsula is part of the Apulia carbonate platform, composed of a well-bedded succession of Jurassic–Cretaceous carbonate rocks, the thickness of which varies from about 3 to 5 km. Besides Salento, this platform includes other three hydrogeological and structural units, namely Gargano, Tavoliere, and Murgia (Polemio & Limoni 2006), which determine a hydrogeological domain hydraulically connected by a deep circulation (Maggiore & Pagliarulo 2004).

The geological basement of the Salento Peninsula, which crops out in extensive areas, comprises layers and banks of fractured and karstified limestone and dolomitic limestone of Cretaceous age; sand, silty or sandy clay, and calcarenite of Miocene to Pleistocene age cover the basement (Ciaranfi et al. 1988). The system is multi-layered comprising seven litho-stratigraphic units (De Filippis et al. 2019), which are spatially distributed in different thicknesses and order. Tectonics involved a combination of gentle folds, E-W dextral and sinistral strike-slip faults, and NW-SE sub-vertical normal faults that dislocated the basement forming Horst and Graben structures. Poor incised ephemeral streams, which are normally dry but discharge significant flow rates and flash floods generated by heavy or intensive rainfalls, cross the karst surface. The networks of the ephemeral streams drain the outflow of exoreic basins along the coasts and converge inland in hundreds of endorheic basins, which cover over 40% of the total area of Salento (Forte & Pennetta 2007). The endorheic area includes several dolines with more or less hydraulically active shallow holes at their bottom. Horizontal and vertical discontinuities strongly controlled in the geological time the development of dolines and other subsurface karst forms such as conduits and horizontal caves in the carbonate rock mass. At the regional scale, the features of the discontinuity network vary from zone to zone. The rainfall depth is low (600–700 mm/year) and precipitation mainly occurs in the autumn–winter period (October–March). Alfio et al. (2020) analysed the Standardized Precipitation Index (SPI) for the Salento aquifer from 1949 to 2011 on the basis of monthly precipitation. They verified the presence of several drought periods during this time, including a variable severity one, from mild to severe that occurred between 2006 and 2008.

The carbonate basement makes up the main coastal aquifer (identified as the ‘deep’ aquifer). The element that marks the hydrogeological and morpho-structural boundary between the Murgia and Salento aquifers is classically represented by the Soglia Messapica, a paleostructure reactivated by tectonics and today buried by sediments of the Plio-Pleistocene cycle (Cotecchia 2014). The transition from Murgia to Salento aquifer occurs through an intermediate buffer zone, located at a lower altitude and exhibiting lower hydraulic heads than both the southern part of the Salento Peninsula and the southern Murgia aquifer. The current NW border of the Salento aquifer is set SE of the (buffer) zone. The highest hydraulic heads in the Salento aquifer (around 3 m a.m.s.l.) occur at its southernmost part and suggest the existence of a regional flow system. Groundwater is mainly in phreatic conditions and flows towards the coastline and the buffer zone, reaching water levels of less than 1 m a.m.s.l. with a mean hydraulic gradient of the order of a few tenths per thousand.

Permeability is locally variable also because of different lithological characteristics of the carbonate rocks making up the basement. The complex structural setting also suggests that lateral groundwater exchanges between the shallow aquifers and the deep aquifer may occur.

Groundwater of the deep aquifer represents the primary water resource for the Salento territory and its recharge depends on direct infiltration of precipitation. The intense groundwater withdrawals have compromised the freshwater–saltwater equilibria, causing groundwater shortage and salinization.

The hydrogeological balance of the Salento aquifer is difficult to define because of the uncertainties in evaluating the coastal discharge and the amount of groundwater abstractions, especially from the unauthorized wells, whose ratio to the legal ones is estimated to 10/1 (Parisi et al. 2018).

The study is based on statistical analyses of precipitation time series originating from 18 rain gauge stations and GWL registrations from 11 monitoring wells in the Salento Peninsula, with the principal purpose to investigate the effects of precipitation on groundwater levels.

Input series comprised daily precipitation records of the 18 rain gauge stations for the period 1 October 2007 to 30 September 2011, provided by the Civil Protection of the Puglia Government (https://protezionecivile.puglia.it/centro-funzionale-decentrato/rete-di-monitoraggio/). Figure 2 illustrates that the spatial distribution of the rain gauges and the corresponding Thiessen polygons, while Table 1 their salient features. Univariate analyses were applied to the 7, out of the 18, rain gauges for which the corresponding Thiessen polygons contained the monitoring wells considered, i.e. Collepasso, Copertino, Galatina, Lecce, Maglie, Nardò, and Ruffano stations (Table 2).

Comparing the precipitation records of the seven selected rain gauges for the defined hydrological period, it was found that the precipitation behaviour is generally homogeneous except for 2007–2008, the driest hydrological year with annual precipitation varying from about 393 to 715 mm. Notably, the hydrological year 2007–2008 shows an average annual value of 549 mm/year, which is significantly lower than that reported by Portoghese et al. (2013) equal to 638 mm/year, for the period 1951–2002. Hydrological year 2008–2009 was the wettest with annual precipitation varying from 857 to 1,144 mm.

The output data are the daily GWLs (average of hourly registrations) registered at 11 monitoring wells in the Salento Peninsula (Figure 1). Data are provided by the Geographical Information System of the Puglia Region (http://www.sit.puglia.it/portal/portale_cis/Corpi Idrici Sotterranei/Dati del Monitoraggio). The selected monitoring wells, which are in a static condition, were equipped with pressure sensors from July 2007 to December 2011 during the development of the ‘Tiziano’ groundwater monitoring project conducted by the regional government to assess the quality and quantity status of the Salento groundwater. The continuous monitoring of GWL concerned 28 wells: administrative inconveniences and changes in technical equipment and personnel caused time gaps in the water level registrations, leading to the selection of 11 out of 28 wells presenting continuous recordings. The depth to groundwater measured with hand-held water level meters at each well every 4 months allowed validating and controlling the accuracy of automatic registrations by comparing measures in the same period. A relative and acceptable difference between measurements confirmed the quality of the continuous recordings. Table 3 provides the location of the 11 wells and their main technical characteristics obtained from historical sources, where available. Table 3 also reports the average groundwater level detected during the observation period. Three wells find groundwater in phreatic condition (LE_12/IIIS, LE_LS21LE, and LE_1/LR); the wells LE_NC4, LE_PS24LE, and LE_PS17LE intercept groundwater in confined condition because of low permeability layers above the carbonate basement. Groundwater at the LE_2/BS, LE_19/IIS, and LE_4/IIS wells results confined at variable elevations below sea level because of basement dislocation under it. Table 3 reports the ranges of elevation of well-casing and screens as well. Any technical information is available for wells LE_SG3 and LE_P1TAU.

Based on data availability, we divided the 11 wells into three groups according to the number of considered hydrological years. The study period is from 1 October 2007 to 30 September 2011 for Group 1, from 1 October 2008 to 30 September 2011 for Group 2, and from 1 October 2009 to 30 September 2011 for Group 3.

For the study aims, we used daily means of hourly water levels. Figure 3 shows the hydrograph of the 11 wells over the studied hydrological years. The evolution of groundwater levels clearly indicates a smooth seasonal change between recharge and recession periods (Cai & Ofterdinger 2016). During the summer dry season, water levels occasionally exhibit a considerable decline in response to excessive abstractions to satisfy drinking and irrigation demand.

Methods based on time series are useful to represent the internal structure of a karst system, the karstification degree, the prevailing karst permeability structure and its functioning.

Time series studies involve both univariate and bivariate analysis in the time and frequency domains (Pulido-Bosch 2021). Mangin (1984) and Box et al. (2013) were the first to implement time series analysis on karst studies.

The autocorrelation and cross-correlation functions have been implemented in a code, which has been written by the authors of this paper through the software R-Studio (R Core Team 2019) according to the formulations detailed below and subdivided into two subsections relating respectively to the time and frequency domain analysis.

The shape of the correlogram provides useful information on a karst system; a gentle slope shows the resilience of the aquifer in terms of infiltration time and/or major groundwater storage. In contrast, a steep slope highlights a more rapid flow through the aquifer or a poor storage capacity (Padilla & Pulido-Bosch 1995; Panagopoulos & Lambrakis 2006; Duvert et al. 2015; Chiaudani et al. 2017).

The time required before the autocorrelation reaches a value of 0.2 is known as memory effect (Mangin 1984). For a random variable, the ACF decreases quickly and reaches zero for very short time lags (Delbart et al. 2016); instead a slow decrease shows an idle nature of the karst, characterized by large storage (Sağir et al. 2020).

Mangin (1984) suggests Cross-Correlation Function (CCF) as useful to analyse the rainfall and discharge relationship to understand the impulse response of the karst system to the recharge process because of rainfall events normally known as a random process.

Defining the two discretized chronological series of length n, ⁠, and so that the first causes the second, the CCF with a truncation point m (k = 0, 1, 2, … m) is not symmetrical.

A significant correlation between input and output time-series at 95% confidence interval is reached when the CCF exhibits a correlation coefficient greater than the standard error 2/n^0.5 (Diggle 1993; Lee et al. 2006; Cai & Ofterdinger 2016).

The CCF indicates a random process with an impulse response: when the CCF shows a maximum for k>0 it indicates an impact of the input signal on the output signal; on the contrary (for k < 0), the output influences the input (Larocque et al. 1998). The delay calculated from a time lag equal to 0, to a time lag corresponding to the maximum value of CCF represents the transfer velocity of the aquifer, which is known as response time (Cai & Ofterdinger 2016).

Spectral analysis allows to break down the total variance of a variable based on the frequency of the events that compose it and therefore to highlight its structure. The spectral approach compared to the correlation approach is twofold; it expresses in the frequency domain what the correlogram reflects in the time domain. Mangin (1984) pointed out that the analysis of a phenomenon in the frequency domain can often facilitate its interpretation. The spectral analysis allows detecting the randomness of a time series since the absence of peaks of the spectral density function characterizes a purely random phenomenon for characteristic frequencies. In the frequency domain, we could obtain different functions to characterize aquifer behaviour, particularly with complex systems where the quick flow superimposes on the baseflow. A short sign of the primary expressions used for conducting the univariate and bivariate analyses in the frequency domain is presented in the following, based on a few studies (Padilla & Pulido-Bosch 1995; Larocque et al. 1998; Panagopoulos & Lambrakis 2006).

The above expression is known as raw spectral estimate (Chow 1969) because it is not a consistent evaluation of the spectral density (Molénat et al. 1999). Thus, a smoothing procedure is performed in order to determine the smoothed spectral estimate (Chow 1969). This procedure permits to give more importance to the periodogram ordinates in correspondence to the interested frequencies, and progressively less weight to the periodogram ordinates at increased frequencies (Chatfield & Xing 2013). In this study, the Daniell or ‘rectangular’ window (Priestley 1981; Molénat et al. 1999) was applied as smoothing procedure.

The spectral density allows defining the duration of the impulse response of the system (Larocque et al. 1998) in terms of regulation time, which is a measure of the inertia of the system; it gives a sign of the length of the impulse response of the system (Lee & Lee 2000). Differently from the memory effect, it is less sensitive both to the sampling interval of the time series and the correlation between distant precipitation events (Zhang et al. 2013).

The CAF indicates the duration of the impulse response function and the filtering of the periodic components of input data, permitting the characterization of the modulating effect of the aquifer in the short, medium, and long term.

The PHF represents the mean delay, for different frequencies, between the input and output signal according to the expression of Padilla & Pulido-Bosch (1995).

COF indicates whether certain variations in the output signal can be attributed to the input one. A loss of coherency at lower frequencies indicates a scarcely karstified aquifer where baseflow dominates; whereas, on the contrary, a high COF also to greater frequencies denotes a highly karstified system. GAF identifies the type of flow and expresses the amplification/attenuation of the input signal: the frequency associated to (i) GAF = 1 represents the duration of baseflow, (ii) GAF = 0.4 indicates the duration of quick flow, and (iii) 0.4 < GAF < 1 characterizes an intermediate flow.

The short-term spectral analysis aimed at investigating the role of precipitation, in terms of intensity and frequency for a few hydrological years, on GWLs measured in the Salento karst coastal aquifer. The analysis used daily precipitation as input and GWL as output considering 11 monitoring wells for a minimum of two up to four hydrological years depending on the well.

In the following, we will show the results of ACF and CCF (time-domain), CAF, PHF, COF, and GAF (frequency) domain analyses. Coherently with Mangin (1984), we fixed a step of 1 day and a truncation point at 125 days for each series.

The ACF of daily precipitation for every considered rain gauge station shows a very fast decrease in 2 days. It lies into a specific confidence interval from (+0.05 to −0.05) for every station and does not show oscillation, because of the variability of the time in which minimum and maximum rainfall occur, even if precipitation shows a seasonal pattern (Angelini 1997). This behaviour confirms that precipitation events are uncorrelated random processes, which can be considered as a pure aleatory function (Panagopoulos & Lambrakis 2006). As an example, Figure 4(a) shows the ACF of precipitation for the Ruffano rain gauge station.

In contrast, the ACFs of GWLs decrease gradually. For some wells, the ACF analysis suggests a very different response if referred to a less rainy hydrological year or a much wetter one (e.g. 2008/09, which is wetter than the average year). Table 4 summarizes the memory effect for every considered hydrological year and well, where the evaluation of the time lag corresponds to the ACF value of 0.2 (Mangin 1984). Calculated time lags are sometimes varying considerably from each other, even for the same well in different hydrological years (e.g. LE_19/II S, LE_PS24LE, and LE_SG 3).

Correlograms show a gentle slope (Figure 5). Figure 5 also provides information about the dependency of correlograms on the intensity and rainfall frequency. The decrease of ACF is even more rapid as the periodicity of the time series becomes shorter or as the frequency increases.

The ACFs of most wells in different hydrological years show negative values (Figure 5): for the Salento aquifer, these values indicate an annual cycle that varies from 30 to 125 days (Larocque et al. 1998).

Considering the frequency domain, the Spectral Density Function highlights aspects of the output signal that ACF cannot clearly distinguish. A comparison of the spectral density of precipitation of the Ruffano rain gauge station (Figure 4(b)) to the spectral density of the GWL of LE_P1TAU allows us verifying that oscillations of the input signal are absent at the GWL spectrum (Figure 6).

Table 4 shows the regulation times evaluated for each well and hydrological year. They range from 50 to 100 days, highlighting the noticeable regulation capacity of the system.

Figure 7 shows the CCF between precipitation and GWL for each well and hydrological year; we can observe that the shape of CCF is dissimilar on a case-by-case basis, with variable amplification. Table 4 shows the cross-correlation coefficients and the response times for every hydrological year and considered well; the evaluation of the time lag corresponds to the maximum values of CCFs (Cai & Ofterdinger 2016).

Furthermore, the CCFs for the Salento wells indicate a bimodal hydrodynamic behaviour (Panagopoulos & Lambrakis 2006), where some of them show a series of peaks after the first peak. Padilla & Bosch (1995) noted a similar aspect, which was attributed to additional flow components arriving in the aquifer after the first peak. All the cross-correlograms show small CCF values (never over 0.24) and a variable response time for the same well in the different hydrological years.

Some wells LE_P1TAU, LE_4/IIS, and LE_1/LR show a negative response time, indicating that precipitation does not have a direct influence on GWL (Cai & Ofterdinger 2016).

During the driest year (2007–08) wells of Group 1 show a fast response time ranging from 1 to 14 days, except for LE_NC4 (56 days). During the wettest hydrological year (2008–09), the wells of Groups 1 and 2 present different behaviour in terms of response time. LE_LS21LE and LE_P1TAU wells exhibit a very fast response time, equal to 3 and 28 days, respectively; during the same year, LE_19/IIS, LE_12/IIIS, and LE_SG3 show a time lag of about 56–57 days, while LE_NC4, LE_PS24LE, and LE_PS17LE have a response time lag of 71 days. For 2009–10 and 2010–11, more similar to the average hydrological year described by Portoghese et al. (2013), the response time varies depending on the well and hydrological year considered.

Wells of Group 3 exhibit different behaviour: LE_4/IIS and LE_1/LR show a response time of 28 days for 2009–10, the wettest of the two considered hydrological years, and 103 and 113 days, respectively, for the second year. On the contrary, LE_2/BS shows an opposite response, i.e. 69 days for 2009–10 and 8 days for 2010–11. This faster response time during a drier year could be associated with a concentrated period of rainfall occurred from October to November 2010, i.e. at the beginning of 2010–11 hydrological year. Figure 3 shows the hyetograph of the large central endoreic area, Maglie gauge station compared to the LE_2/BS GWL hydrograph: precipitation generates a piston flow that causes water to quickly drain towards the deep aquifer.

The CAF, for each well and hydrological year, shows different results. In general, for the hydrological year 2007–08, which is the driest among the considered, the CAF reaches a maximum value of 1, in contrast to the rest of the studied years that exhibit highly variable amplitudes greater than 1. However, for 2008–09, 2009–10, and 2010–11 hydrological years, the CAF reaches the value of 1 at frequencies between 0.05 and 0.1; subsequently, the CAF tends to zero demonstrating a common slow hydraulic behaviour of the Salento aquifer. Therefore, the CAF shows a series of peaks that progressively decrease for high frequencies, reflecting the same shape of the precipitation spectral density, and recognizing and conserving the main frequencies associated with the rainfall events. Figure 8 shows a comparison between the CAFs of LE_P1TAU well for the wettest (Figure 8(a)) and driest (Figure 8(b)) hydrological year and the corresponding spectral density of precipitation of the Ruffano rain gauge station. Both during the 2008–09 and 2010–11 hydrological cycles, the CAF tends to zero long before the frequency of 0.1.

Regarding the PHF, not in all cases there is a linear trend at low frequencies that allows us to calculate the average time lag (Figure 9(a)–9(c)). Results (Table 4) show that this aspect depends on the hydrological year and not on the monitoring well (e.g. for LE_NC4 and LE_P1TAU, the average lag cannot be calculated for all hydrological years).

The COF is a measure of the input–output linearity. The average COF for each well and hydrological year shows a fluctuating behaviour indicating poor linearity between precipitation and groundwater levels; this poor linearity is more accentuated in the drier hydrological years. In general, for the wells of Group 1, the coherency is more variable and insufficient for the 2007–08 hydrological year, while better values are obtained for the rest of the considered cycles, with values ranging from about 0.2 to 0.7 at frequencies not greater than 0.15. The same results are obtained for the wells of Group 2 in relation to the hydrological year 2008–09; scarce and in some cases undistinguishable coherency values are found for the 2009–10 and 2010–11 hydrological years.

Figure 10(a) illustrates as an example the COF for the LE_NC4 well, corresponding to 2007–08 and 2009–10 hydrological years; Figure 10(b) shows the result for the LE_P1TAU regarding 2008–09 and 2009–10 periods, demonstrating the fall of the coherency according to the precipitation in the selected time window. As for the previous cases, Group 3 shows the same trend of COF with values oscillating between 0.2 and 0.7 for 2009–10 and 2010–11, respectively and for a frequency of 0.2 and 0.1.

The GAFs for each well and hydrological year yield values that are never higher than 0.1, thus indicating the attenuation of the input signal at low frequencies (Zhang et al. 2013) independently from the precipitation trend of every hydrological year. Figure 11 illustrates the GAF only for the wells LE_NC4 and LE_P1TAU, based on the available data over the considered hydrological years.

Lastly, we verified whether by using longer time-series there was a more evident link between precipitation and aquifer response. To this end, we carried out the long-term ACF analyses for each well considering the entire time series with a cutting point of 125 days. It was revealed that, in most cases, the response of GWLs in terms of memory effect is close to, or longer than 100 days (Figure 12(a)). Only wells LE_1/LR and LE_P1TAU exhibit shorter lags, from 57 up to 82 days.

The CCFs of the long-term analyses show a statistically poor correlation between precipitation and GLW for all wells: they manifest a similar general behaviour and correlation coefficients not greater than 0.1 (Figure 12(b)). Although, in some cases, the correlation peak corresponds to a small response time of less than 25 days (e.g. LE_LS21LE, LE_4/IIS, and LE_12/IIIS), all CCFs slowly goes to zero with numerous supplementary peaks, confirming the long-time lags with which the aquifer responds to precipitation events. Contrary to the previous ones, wells LE_19/IIS, LE_NC4, LE_PS24LE, LE_PS17LE, and LE_P1TAU show a response time of about 75 days.

We first consider the results of the short term series analysis in the time domain. The results of ACF allow asserting that baseflow is the major water transfer process that plays the key role because of a modest connection between fissures and fractures of the rock matrix. However, when the precipitation increases above the average of 638 mm/year (Portoghese et al. 2013), the ACF firstly shows a steep slope representing the quick flow through karst conduits and large fractures; the baseflow process emerges only much later thanks to minor fractures and discontinuities that drain the unsaturated zone (Panagopoulos & Lambrakis 2006). A few wells (LE_NC4 and LE_2/BS), more clearly show that the baseflow is the unique water transfer process, independently from rainfall intensity and frequency; this indicates that some parts of the aquifer are much more inert.

We may assume in general that sub-horizontal flow is the prevalent component of hydrodynamic evolution at Salento aquifer. Based on the memory effect and the ACF values calculated for all wells, the hydrodynamic functioning of the Salento aquifer resembles more that of a porous aquifer rather than a karst aquifer, as proposed by Angelini (1997) for two karst springs of Central Italy.

Evidence of those aspects highlighted by the ACF is provided by the Spectral Density Function. The spectral density of GWLs strengthens the conceptual model of functioning of the Salento aquifer: it weakens the high frequencies of the precipitation input acting as a low-pass filter (Molénat et al. 1999; Imagawa et al. 2013). The large storage capacity of the Salento aquifer is further supported by the regulation times (Table 4). They range from 50 to 100 days, highlighting the noticeable regulation capacity of the system.

Overall, the CCFs decrease very slowly, demonstrating the inertial behaviour of the Salento aquifer. A comparison with the results obtained by Delbart et al. (2016), indicates that, concerning the water transfer, the Salento aquifer is a slow response aquifer, where the rock matrix with fractures and fissures appear to play the controlling role.

However, at the local scale, even within the framework of this general behaviour, each well shows some differences, as expected because of the complex features of the coastal karst aquifer.

The different response times of the considered wells are, as already highlighted, a function of the amount of precipitation but also of the presence or absence of surface karst forms, which control the link between surface and groundwater. The wells showing the fastest response (LE_12/IIIS, LE_2/BS, LE_LS21LE, LE_PS24LE, and LE_P1TAU) are in the Salento recharge area, the large central endorheic area, characterized by numerous dolines, a more incised drainage network and overall a system of major faults behave as conduits for flow parallel to the fault planes. In fact, the above-mentioned wells are located nearby dolines showing the direct link between superficial drainage and groundwater recharge (Mastronuzzi & Simone 2015). Thus, the presence of these karst forms has an influence on the rapid aquifer response time. To further validate these results, investigations related to the hydraulic role of the major faults are required. The CAF analyses indicate a small influence of the quick flow on the Salento aquifer for each of the 11 analyzed wells, highlighting the importance of baseflow. Nevertheless, during the wettest hydrological year, although with very low cross amplitude ordinates, the CAF systematically reproduces the peaks of the associated precipitation spectral density, which suggests quick flow may be contributing significantly to Salento's hydrodynamic evolution under such hydrological conditions. On the contrary, for the driest year, the CAF tends to zero at frequencies much less than 0.1 and does not clearly show peaks as in the wettest hydrological year.

The PHF results indicate a strong influence of the hydrological year on them. The absence of a linear trend depends on the time series considered and not on the analyzed well. Moreover, for the medium–high frequencies, the PHFs show that the input signal is much attenuated, distorted, and inconsistent.

In general, COF results underline and confirm the role carried out by the rock matrix in damping the input signal: this damping causes for the Salento karst aquifer poor correlation and coherency between precipitation and GWLs. These results strengthen the idea that the water transfer at the regional scale occurs by redistributing water through sub-horizontal pathways.

Results of GAF confirm the role of inertial filter of the Salento aquifer, where the matrix plays a major role in slowly releasing the infiltrated water (Benavente & Pulido-Bosch 1985) and attenuating the input signal at low frequencies, independently from the precipitation course of every hydrological cycle. GAFs for each well and hydrological year, indeed, show values never higher than 0.1.

The results of the ACF long-term analyses reveal that in most cases the response of GWLs in terms of memory effect is longer than 100 days, with few exceptions as shown in Figure 12(a).

The common behaviour of the long-term ACFs outlines the important role of the rock matrix: the slow release of the infiltration water provides a large buffering capacity that acts as a natural means to effectively preserve the aquifer system's reserves. In fact, the Salento aquifer exhibits great storage capacity in time. As to the long-term CCF, the absence of statistically significant correlation between precipitation and groundwater levels related to LE_SG3 and LE_1/LR, and the low cross-correlation detected for the remaining wells, suggest that the Salento aquifer, contrary to most karst aquifers where exploitation occurs by the uptake of springs, may be influenced by additional factors, as the heavy exploitation of groundwater by wells.

There are no data describing the degree of karstification of the Salento aquifer, which is why the time series analysis was chosen. However, the only data coming from the stratigraphy of the monitoring wells often reveal low permeable strata made up of low permeability carbonate levels as shown in Table 3. This generally confirms the result obtained and provides a new key to understanding the hydrodynamic behaviour of this peculiar karst aquifer.

The coastal state of the aquifer poses other questions. In continental aquifers, the correlation between precipitation and GWL is often significant (Zhang et al. 2013; Delbart et al. 2014; Mayaud et al. 2014; Pavlić & Parlov 2019). In the Salento coastal karst aquifer, the flow system does not have a permanent lower limit (as a geological impermeable limit) but, for its entire extent, an unsteady hydraulic boundary between freshwater and saltwater; this, in principle, implies that the response of GWL to precipitation is much more complex. In a recharge area of this coastal aquifer, when infiltration waters reach the saturated zone, the water level increase disturbs the freshwater–saltwater equilibrium and the buoyancy effects locally change the transition zone elevation, while groundwater spatially redistributes over time on the regional scale. The overall buoyancy effects may control the extent of GWL and its oscillations, being likely able for the same reason to hide the effects of overexploitation. Another concern in coastal aquifers is the effects of GWL fluctuations on water quality. Decreases in levels (seasonal or related to exploitation) cause worsening of the density distribution along the vertical of the aquifer (upwards propagation), which, due to the high frequency of fluctuations, is not easily restored during level rise phases (Fidelibus & Pulido-Bosch 2019). Thus, even when over exploitation conditions imposed on the aquifer ceases and water balance is restored, it is not certain that the quality of the resource has not been, and remains, compromised.

Even under these limitations, the results emphasize and confirm the complexity of the permeability structure of the Salento aquifer and its high anisotropy. Also, they indicate a poor degree of karstification, and a not well-developed hierarchical karst network characterized by weak water drainage capacity. Such features establish dominance of horizontal flow conditions and guarantee large storage capacity, as also suggested by Panagopoulos & Lambrakis (2006).

In conclusion, the recognized hydrodynamic behaviour based on the analyzed monitoring wells of the Salento aquifer is unimodal with a dominant water transfer process due to baseflow. Only occasionally this behaviour is bimodal with the concurrent presence of quick flow and baseflow. The recognition of the quick flow may be the proof of the existence of an epikarst. According to the conceptual model, during precipitation water infiltrates through the karst surface reaching the epikarst; there, because of its high storage capacity, infiltrating water can be temporarily retained during its migration to the saturated zone (Mangin 1975; Bailly-Comte et al. 2008). Water stored in the epikarst during periods of low recharge as in the summer period is released under the pressure of autumn–winter precipitation: epikarst starts behaving as a perched aquifer with a horizontal flow towards conduits, sinkholes and master joints. Its existence and role was proven in the neighbouring Murgia karst coastal aquifer studying the water transfer processes through natural tracers in its principal recharge area, also characterized by major discontinuities and endorheic basins (Fidelibus et al. 2017).

Thus, in an area where the census lists 280 swallow-holes, 1,172 dolines with a variable density ranging from 0 to 7 km², 223 caves, numerous major discontinuities and hundreds of endorheic basins, quick flow also points out the existence of epikarst, the role of which in the transfer process of the infiltration water only becomes apparent during the time series analysis when the system is subjected to high precipitation stress. The role of major faults in the water transfer process could be connected to the evidence of the quick flow, because they normally enhance permeability, resulting in increased recharge (Wendt et al. 2020). However, considering their role in groundwater flow, faults may behave as barrier to flow when filled with clay materials and/or red soil deposits, or by secondary cementation, or as effective conduit when dissolution along the fault or fracture planes is dominant. Recent studies (Fidelibus & Pulido-Bosch 2019) allowed hypothesizing that, in the Salento aquifer, many of the faults crossing the overburden and corresponding to the buried karst behave simultaneously as a barrier for horizontal flow transversely to the faults themselves, and as conduits for flow parallel to the fault planes. Other major discontinuities prove, instead, to have very low conductive capability.

The research provided an interpretation of the hydrodynamic behaviour of Salento aquifer, as deduced by the identified relationships through the application of short time series analyses, between precipitation from 7 gauging stations and GWL originating from 11 monitoring wells spread in the Salento coastal karst aquifer. This study was carried out considering different hydrological years in addressing the influence of precipitation frequency on the aquifer response.

The results evidence the capability of these methods to describe the hydrodynamic characteristics even of a coastal karst aquifer, far from typical karst aquifers considered in previous studies. Even though the analyses considered the variable response time of each individual well subject to the local conditions, as a whole the significant inertia of the Salento karst system in terms of transmissivity and significant storage capacity was deduced.

The analyses show the role of the rainfall intensity and frequency on the hydrological system: differences in the precipitation pattern trigger hydrodynamic mechanisms at the epikarst and unsaturated zones.

The considered period does not allow extending the analysis on the effects of a single precipitation event and its influence on the aquifer recharge process. The response of Salento aquifer to single rainfall events and their size is however worthy of further study in the light of the occasionally observed quick flow and the potential associated groundwater vulnerability (Fidelibus et al. 2017; Delle Rose & Martano 2018).

Local conditions (because of geology, surface morphology, permeability structure of the unsaturated zone, and precipitation pattern) are key factors in shaping the GWL evolution. However, when recharge water reaches the saturated zone, its move is controlled by the prevalent saturated permeability structure and is likely to migrate over large areas by (often preferential) horizontal flow. We could parallel the aquifer behaviour to that of a water distribution network supplied by several reservoirs connected to different nodes. These reservoirs supply the system with water heads that depend on space and time with a consequent time-dependent redistribution of water pressure. In this context, the water table of the Salento aquifer would be continuously subject to reshaping due to variable increase or decrease of hydraulic heads, which depend on the spatial distribution of precipitation and features of recharge areas.

Even with these limits, we rely on results that reveal a slow response to recharge of the Salento aquifer, a high storage capacity, and an annual cycle achieved anyway within a year in each considered case.

The results of this study describe the Salento aquifer as a slow system in terms of response to meteoric events but, at the same time, characterized by a relevant memory effect. The slow release of infiltration water indicates a large storage capacity for the Salento aquifer, which guarantees, despite ongoing climate change, the water reserve for drinking and irrigation uses for this peninsula characterized by the total absence of surface water resources. This aspect is interesting and favourable from the point of view of the water balance but does not exclude concerns about the groundwater vulnerability to seawater intrusion that is strictly linked to GWLs variations. Evidently, coastal aquifers, as complex systems, comprise many strongly interdependent variables, that can result in a non-linear behaviour and exhibit feedback loops: because they are dynamic systems, they change over time, and present-day states depend on previous states. This last characteristic suggests attention should be focused on the ‘evolution’ of the system, which, in turn, may only be studied through analysis of time series from carefully cited monitoring points that record data at appropriate frequencies to capture the coastal aquifers’ responses. To this end, continuous recording of precipitation and GWLs is of great value.

This research is conducted in the context of MEDSAL Project © (www.medsal.net), which is part of the PRIMA Programme supported by the European Union's Horizon 2020 Research and Innovation Programme and funded by the national funding agencies of GSRT (grant number 2018-7), BMBF, RPF, MIUR (grant number 1421), MHESR (grant number 2018-12), and TÜBİTAK (grant number 118Y366).

All relevant data are included in the paper or its Supplementary Information.