Factors Affecting Water Drainage Long-Time Series in the Salinized Low-Lying Coastal Area of Ravenna (Italy)

Abstract

The low-lying coastal area of Ravenna (North-eastern Italy), like the majority of delta and coastal zones around the world, is affected by groundwater salinization due to natural processes (such as low topography, natural land subsidence, seawater encroachment along estuaries, etc.) and anthropogenic activities (i.e., increased anthropogenic subsidence rate, sea level rise, geofluids extraction, and drainage). Among all factors causing aquifer salinization, water drainage plays an important role in lowering the hydraulic head and favouring saltwater seepage in the Ravenna coastal aquifer. A network of drainage canals and water pumping stations first allowed for the reclamation of the low-lying territory and today are fundamental to keep land and infrastructures dry and maintain effective soil depth for agriculture practices. The aim of this work is to identify and assess factors affecting water drainage long-time series (1971–2017) of the most important mechanical drainage basin in this low-lying coastal area. Statistical analyses of drainage, climate, and land use change datasets help constrain the relative weight of each single factor potentially causing an increase of water drainage through time. The results show that, among these factors, subsidence rates and seepage processes are the most significant. The data trends also indicate that the climate, especially in terms of precipitation amount and extreme events, played no important role during the studied time interval. The process of infiltration soil capacity loss due to urbanization and consequent soil sealing probably has a small secondary effect. Moreover, an increase in pumping through time will exacerbate aquifer salinization and compromise freshwater availability in the coastal area.

Highlights:

• We assess factors that may cause drainage increase in a low coastal plain;
• Land subsidence related seepage causes most of drainage increase;
• Climate and land use change have minor effects on increase of drainage;
• Increase in drainage exacerbates aquifer salinization in low-lying coastal areas.

1. Introduction

Coastal areas are potentially vulnerable to climate change, land use change, land subsidence, sea level rise, groundwater and soil salinization. Salinization processes have occurred to some degree in many of the coastal aquifers of Italy [1,2,3,4,5,6,7]. Such processes can seriously affect coastal populations and their activities [8,9,10]. Coastal vulnerability to water and soil salinization is especially high in low coastal plains, where the coastal basins are at or below sea level, the so-called polder basins [11]. Polder basins need to be constantly drained to keep the coastal land dry; this is accomplished by pumping stations that supply the mechanical energy for drainage water that gravity cannot provide in a basin below sea level. Coastal land drainage in coastal basins and its relationship to an increase in groundwater salinization has been recognized throughout the world [12,13,14], as well as locally in the low-lying coastal area of Ravenna (North-eastern Italy; [1,15]). Coastal areas are expected to be more vulnerable to future changes in sea level and salinization than natural systems especially where the water table is controlled by drainage systems and where the inland hydraulic head is kept unchanged [16,17]. Moreover, heavy drainage systems can be responsible for the process of autonomous salinization, by creating a vertical gradient and forcing groundwater to flow from the bottom toward the top of the aquifer. Therefore, connate saline groundwater previously preserved in the deepest portion of the aquifer can be mobilized and transported upwards [18]. While drainage and lowering of the hydraulic head has been recognized as one of the main causes of Ravenna coastal aquifer salinization [18,19], no drainage data analysis and assessment of the factors affecting the increase of water drainage through time have been performed so far.

A unique long time series dataset of drained water in a polder basin within the low coastal plain of Ravenna allowed us to study the change in drainage rates through the time from 1971 to 2017 and to highlight an unexpectedly high positive trend. The causes that affect the increase in drainage may be multiple, such as a change in the precipitation patterns or climate extremes, a change in land use that allows for less infiltration and more runoff, or land subsidence that causes an increase of the water volume to be drained from the soil [9,20,21,22,23].

The main objective of this work is to identify the processes, and their relative contribution, which are responsible for the drainage increase detected in long time series of pumping rates from the land reclamation stations of the Quinto Basin in the Ravenna coastal plain. This basin is the most important mechanically drained basin in the area and could be a proxy for many other basins of the Northern Adriatic Coast and throughout low coastal plains around the world [24,25].

Our study is also preliminary to identify mitigation strategies for integrated coastal and water management and to devise some fair compensation plans to charge those stakeholders running activities responsible for the increase in pumping rates and drainage [26].

2. Materials and Methods

2.1. Study Area

The Quinto Basin (Quinto Bacino) is a coastal basin of about 10 km² located south of the urban area of Ravenna (Italy, Figure 1). It is bordered to the South by the Bevano River, to the North by the Fiumi Uniti River, to the East by the Adriatic Sea, and to the West from part of the Ronco River. The basin has an area of 9252 hectares and includes several environments, such as a beach-dune system, coastal and historic pine forests, and a large agricultural area with several gravel pits. Once abandoned, the gravel pit lakes became artificial brackish-saline water bodies. Most of this area suffers from aquifer salinization [19,25,27,28,29].

Most of the Ravenna territory is a reclamation area and subject to mechanical drainage [30]. In total, there are 15 drainage basins: eight coastal basins require mechanical water drainage by pumping stations in order to discharge water toward the Adriatic Sea, while the other seven basins, in the most inland and topographically highest portion of the city, have a natural hydraulic gradient. Drainage channels, connected to numerous water pumping stations, allowed the reclamation of the topographically lower areas in the past and are nowadays fundamental to lower water table, guaranteeing agricultural activity and avoiding flooding of the territory. The hydrographical regime is entirely artificially controlled and managed by the Land Reclamation Consortium, which has divided the territory into several mechanical drainage basins (Figure 1). The aim of this reclamation system is to keep the water table about 1.5–2 m below ground level during the year (levels may slightly vary in winter and summer). In the coastal area the water table is on average close to sea level with a general inland-directed hydraulic gradient mainly forced by the drainage system [19]. Moreover, the drainage causes a decrease in the effective recharge of the aquifer [31] and a vertical hydraulic gradient that forces saline groundwater to flow upward from the bottom of the aquifer [1,18,31].

Another phenomenon that has always played a fundamental role in the evolution of the study area is land subsidence. The natural component of land subsidence is due to deep tectonic causes and compaction of the Quaternary sediments. The anthropogenic activities (urbanization and above all geofluids extraction) during the past century have significantly increased the rates of ground sinking. Geofluids extraction has occurred from deep reservoir and aquifers well below the surface aquifer that is in contact with the network of drainage channels. The surface aquifer, therefore, has not undergone any significant compaction in the last 50 years, so that its hydraulic properties have not changed. More details on the subsidence rates of the last 100 years for the study area can be found in Teatini et al. [32], Carminati and Martinelli [33], Gambolati et al. [34] and Antonellini et al. [35]. Among several impacts, land subsidence has led to the rise of water table in the shallow unconfined aquifer that has consequently required an increase in mechanical pumping drainage rates across Ravenna to avoid flooding of large areas.

2.1.1. Hydrogeological Setting

The study area is in the South-Eastern portion of the Po River plain. The sedimentary dynamics of this area are mainly driven by climatic and eustatic fluctuations and the geological evolution of the region was governed by the alternation of continental and marine sedimentation as well as erosive and depositional phases [36,37]. Sediments of the upper shallow coastal aquifer were deposited during the Holocene sea-level fluctuations [38,39] when a barrier-lagoon-estuary system was moving eastwards with the receding sea resulting in an alternation of highs and lows in the topography, which correspond to different dune belt (palaeodunes) and coastlines (Figure 2). As a result, the shallow coastal aquifer consists of a wedge-shaped dune and beach sand body pinching out in westerly and northern directions, which is sealed at the bottom and top by clays and peat formed in lagoons, marshes, and alluvial plains (Figure 3). There are two main units: a relatively thick medium-grained sand unit at the top of the aquifer (about 10 m in thickness) and a lower fine-grained sand unit of small thickness (from −20 to −25 m. a.s.l.). These two sand bodies are separated by an alternation of clay-silt and sand-silt layers in the coastal area and merge landward where the prodelta sediments disappear. The beach gravel deposits are found at the inland tip of this wedge; gravel deposits were excavated for commercial use (gravel pits). The quantity of sand decreases in the western part, towards the agricultural area, where the aquifer manly consists of recent alluvial plain deposits and back-barrier facies association, which includes fine-grained sediments (clay, silty clay with thin sand layers intercalated) with abundant organic matter, typical of lagoon and marsh environments [40]. Lastly, the Flandrian continental silty-clay basement is at a depth varying from −15–−20 m., in the west, to −25–−30 m. a.s.l., at the present shoreline [40,41,42].

Given the heterogeneity of the aquifer [40], hydraulic conductivity (K) values vary from those typical of medium sand grain size in the shallow aquifer (10⁻³ < K < 10⁻¹ cm/s), to fine sediments of prodelta (10⁻⁶ < K < 10⁻⁴ cm/s), and finally to those typical for the deep fine sands of the aquifer (10⁻⁵ < K <10⁻³ cm/s) [15,43].

2.2. Data Acquisition and Elaboration

The objectives of data acquisition and elaboration were to identify trends in the change of drainage through time and to correlate them with changes in climate time series, land use, land subsidence rates, and the hydraulic of groundwater flow in the study area.

2.2.1. Drainage Data

The daily and annual drainage data and the number of drainage days in a year from 1971 to 2017 at the Quinto Basin pumping station were obtained from the Land Reclamation Consortium and used to assess data trends.

Correlation coefficients (r) between daily pumping and daily precipitation were computed at different lag times (t = 0, 1, 2… 6 days) for each year and over two different time intervals (1971–1981 and 2006–2016) to assess an eventual change in the efficiency of the drainage system. The two different time intervals were chosen as representative of two different urbanization status: minimum (1971–1981) and maximum urbanization (2006–2016). The assumption is that an increase in urbanization in recent years, in fact, would give a better correlation at smaller lag times with respect to the original land use situation, because of the increase in impervious surfaces. Also, an eventual change in the distribution of the rainfall over time (increase in intensity of rainfall phenomena) could cause a change in the lag time for best correlation.

We performed time series analysis of the rainfall, average temperature, potential evapotranspiration (Thornthwaite and Mather method [44]), real evapotranspiration (Turc method [45]), and pumping data (annual pumping, drainage days in a year, max daily pumping per year) by applying Kalman filtering to obtain a space state model, which eliminates the irregular component in the time series (package KFAS in R; [46,47]). Kalman filtering belongs to the exponential family of state space models [47]. State space models allow modelling of several types of time series and other data (structural time series, autoregressive integrated moving average (ARIMA) models, simple regression, generalized linear mixed models, and cubic spline smoothing). The type of Kalman filtering used in our paper is a simple linear Gaussian state space model that is analytically tractable and often used in hydrology and many other fields of science [47]. This filter analytical formulation is given by:

y(t) = Z(t)α(t) + ε(t) (observation equation)

(1)

α(t + 1) = T(t)α(t) + R(t)η(t) (transition equation)

(2)

where the state disturbances ε(t) ~ N(0, H(t)), η(t) ~ N(0, Q(t)) and α(1) ~ N(a(1), P(1)) are independent of each other; y(t) is the observation matrix; Z an array corresponding to the system matrix of observation equations; H an array corresponding to the covariance matrix of observational disturbances Δ; T an array corresponding to the first system matrix of state equations; R an array corresponding to the second system matrix of state equation; Q an array corresponding to the covariance matrix of state disturbances η; a1 a matrix containing the expected values of the initial states; P1 a matrix containing the covariance matrix of the non-diffuse part of the initial state vector; P1 inf a matrix containing the covariance matrix of the diffuse part of the initial state vector.

All system and covariance matrices Z, H, T, R, and Q can be time-varying, and partially or totally missing observations y(t) are allowed. Covariance matrices H and Q must be positive. Please see the original Koopman and Durbin [46,47,48] articles describing the methodology and the Helske implementation in the R package description [49] for a detailed explanation of the Kalman filter application also for the case in which the series in the model is defined as non-Gaussian.

2.2.2. Land Use Change

Land use change and its eventual effect on drainage data was evaluated using land use maps of 1976 and 2014 elaborated by the Emilia-Romagna Region, Geological Office. The classification of land use cover is coherent with the directives of the Corine Land Cover [50], and it is based on five macro-classes and specifically: Class 1—Urban areas, Class 2—Agricultural land, Class 3—Forests and bushland, Class 4—Wetlands, and Class 5—Surface water.

Land use change was assessed in ArcGIS^® environment. The area of the different land use classes was computed for the two reference years 1976 and 2014 and compared by using the Patch Analyst, and specifically the Patch Grid—Intersect (Combine) Grid tool (http://www.cnfer.on.ca/SEP/), which facilitates the spatial analysis of landscape patches and modelling of associated attributed. Using this tool, it was possible to assess changes in land use cover from 1976 to 2014 with three evolutionary classes: (i) changes from classes 2, 3, 4, and 5 to class 1: urbanization—loss of infiltration capacity; (ii) changes including all classes apart from class 1: other variation; and (iii) no change of land cover class: unchanged.

The total change (%) for each land cover class (Δ_tot) from 1976 to 2014 was calculated as follows [51]:

$${\mathsf{\Delta }}_{tot}=\frac{Area Class n \left(2014\right)-Area Class n \left(1976\right)}{Area Class n \left(1976\right)}100$$

(3)

where n = 1, 2… 5 is the class number.

The artificialization index (I_art), which expresses the intensity of urbanization, was computed as:

$$I_{art}=\frac{Area Class 1 \left(2014\right)-Area Class 1 \left(1974\right)}{Total Basin Area}100$$

(4)

All areas are expressed in hectares.

2.2.3. Climate Data and Climatic Extremes

Climate data (precipitation and temperature) were analysed to assess if a change in precipitation pattern or intensity could have affected the pumping data time series. The meteorological daily data (precipitations—min, max, and average temperatures) from 1971 to 1975 were obtained from the Hydrological Annals published on the ARPAE website (www.arpae.it), while data from 1976 to 2017 from the database DEXT3R (http://dexter-smr.arpa.emr.it) for 3 stations in Ravenna. The monthly potential evapotranspiration (ETP) was computed with the Thornthwaite equation [44], while the annual actual evapotranspiration (ETR) was computed with the Turc equation ([45,52]; Supplementary Materials, Table S1).

A change in the intensity and pattern of rainfall related to climate change in the Ravenna area was evaluated using the indexes for extreme climate events published by ISPRA [53] for the Italian territory. The indexes for extreme rainfall events that we used are the following:

• R10—Number of days in a year with intense rainfall ≥10 mm.
• R20—Number of days in a year with very intense rainfall ≥20 mm.
• Wet days WD—Number of days in a year with rainfall ≥1 mm.
• SDII—Daily precipitation intensity obtained by dividing the cumulative annual precipitation by WD in the year.

The trends of the rainfall time series and trends of the indexes for extreme rainfall events were obtained from Kalman filtering [46] and compared against the drainage dataset.

2.2.4. Subsidence Data

The vertical land subsidence rates were obtained from different sources. The data of Teatini et al. [32], acquired with traditional leveling techniques (DGPS), were used for the period 1972–1992 after importing the original maps and digitizing them in a GIS environment for data elaboration. The average land subsidence rates were obtained for the periods 1972–1977, 1977–1982, 1982–1986, and 1986–1992 by averaging the data over the basin area.

The land subsidence rates for the period 1992–2016 were obtained from the published maps of ARPAe Emilia-Romagna (www.arpae.it/index.asp?idlivello=1414). The maps were elaborated based on GPS and satellite radar interferometry data (PSInSAR) [54] and are available for the intervals 1992–2002, 2002–2006, 2006–2011, 2011–2016. Also, in this case, an average land subsidence rate for each sub-period was obtained by elaboration in ArcGIS©.

Since Carbognin et al. [55] state that the area was relatively stable in recent years, the land subsidence rate for 2017 was considered the same as that of the period 2011–2016. The average land subsidence rates for the different periods were used as input for the model of vertical seepage in (Section 2.2.5) that can be related to the increase of pumping (ΔPu expressed in mm) through time.

2.2.5. Seepage Calculation

The vertical head gradient between the bottom of the aquifer and the water table elevation, which is below sea level, determines the seepage flux in a polder basin [56,57]. The heads at the bottom of the aquifer are influenced by the regional system (the surrounding heads including sea level), so that we can neglect the horizontal head gradient from the sea towards the basin in computing the seepage (Figure 4). Land subsidence causes an increase in vertical head gradient (Δh) between sea level and inland water table. Being a coastal unconfined aquifer, the difference in head increases whereas the thickness of the aquifer b remains constant (Figure 4). As a result of subsidence, therefore, there is an extra amount of seepage water entering the Quinto Basin that can affect the amount of water that has to be drained by the water pumping station. More than 20% of the Quinto Basin topography is at or below the sea level and the land reclamation pumping machines keep the water table at a depth of −1.5 to −2.0 m (a.s.l.) for allowing agricultural activities and avoiding low areas flooding. In these conditions, there is a vertical hydraulic gradient, which continuously increases as the drainage rate keeps pace with the land subsidence rate [1].

In order to obtain a first order estimate of the extra vertical seepage caused by the land-subsidence-related increase of vertical head difference, we use the following model based on Darcy law:

$$\mathsf{\Delta }{Q}_v=A{\int }_0^{t}A\left(t\right)q_{v}dt=K_v{\int }_0^{t}A\left(t\right)\frac{h\left(t\right)}{b}dt$$

(5)

where ΔQ_v is the increase in total vertical seepage from time 0 to time t (we set 1971 as t = 0 and the final t is 2017); q_v is the vertical seepage per unit surface area (mm year⁻¹ m⁻²); K_v is the hydraulic conductivity in the vertical direction (m/day); h(t) is the hydraulic head (m) as a function of time due to land subsidence in the period 1971–2017; A(t) is the area (m²) of the Quinto Basin affected by vertical seepage, where the topography and water table are below sea level and that is also changing through time as a result of land subsidence (~2000 ha); b is the thickness of sand units (m) of the unconfined surface aquifer, which we assume constant as land subsidence is the result of deeper processes happening in the deep portion of the aquifer [35]. We only consider the vertical seepage occurs in the most conductive sandy units (thickness ranging from 10 m, close to the coast, to 15 m in the most inland part of the basin, [27]). The assumption of uniform upward seepage through the aquifer is an approximation (see conceptual sketch in Figure 1, which is based on the finite elements numerical model) and indeed, a few studies [56,57,58,59] have shown that groundwater seepage could occur through preferential pathways as permeable, sandy units or paleochannel belts [56] where groundwater discharges vertically at high velocities.

The vertical hydraulic conductivity K_v is calculated by using arithmetic, geometric and harmonic means [60] of horizontal hydraulic conductivity values of sand units and prodelta as reported in Section 2.1.1. [15,43] and considering an anisotropy ratio, vertical versus horizontal hydraulic conductivity, K_v/K_x, of 1/3 [15]. The values of K_x considered are 1 m/day and 0.1 m/day for the shallow medium sand unit (10 m thick) and the fine sand unit (5 m thick), respectively, and 0.01 m/day for the prodelta unit (10 m thick) [15,27].

Finally, vertical seepage values are calculated for the ranges of sand unit thickness and K_v values that characterized our study area and then compared to the estimate from water budget considerations.

2.2.6. Water Budget for the Quinto Basin

A change in the components of the water budget between 1971 and 2017 can be done assessing the variations of the individual components. The water budget equation for the beginning of the examined period can be written as

$$P+Q_v=Pu+ETR$$

(6)

while its variation through time and at the end of the examined period as

$$\mathsf{\Delta }P+\mathsf{\Delta }{Q}_v=\mathsf{\Delta }Pu+\mathsf{\Delta }ETR$$

(7)

where P is the precipitation (mm/year), Q_v is the vertical seepage (mm/year) that is always considered in a polder basin, Pu is the land reclamation pumping (mm/year), and ETR is the real evapotranspiration (mm/year). The Δ quantities refer to the variations over the time period considered.

Using the model trends of the available time series, we can estimate ΔP, ΔPu, and ΔETR, so that we can evaluate ΔQ_v. Once obtained ΔQ_v from water budget considerations, we can compare it with the value of ΔQ_v obtained from the model of Equation (5).

3. Results

3.1. Time Series Analysis

The results of time series analysis on drainage (pumping), precipitation, potential evapotranspiration, real evapotranspiration over the period 1971–2017, and indexes of extreme rainfall in the Quinto Basin are shown in Figure 5, Figure 6 and Figure 7, as well as Figures S1 and S2 in Supplementary Materials. In all figures, we express the total equivalent amount of water in mm/year by dividing the total volume of water by the area of the basin.

Figure 5 shows the space state models (red lines) obtained by applying Kalman filtering to the pumping dataset. The annual drainage (Figure 5a) reaches its maximum value (469.7 mm) in 2010, and the minimum (57.2 mm) in 1988. The annual drainage tends to increase significantly over the period 1988–2017, with a cumulative total increase (ΔPu) of 150 mm. The number of operation days at the pumping station (Figure 5b) indicates a minimum of 113 days drainage days in 1989 and a maximum of 365 days in 2016. The number of drainage days per year increases by 141.4% in 46 years, indicating a positive trend. In 2015, the maximum (365 days) has been reached.

Figure 5c indicates that the maximum daily pumping per year has not increased significantly as shown by the trend of the space state model (red line), meaning that it has been compensated for by the increase in drainage days.

The analysis of the rainfall dataset indicates an average precipitation of 640 mm as also confirmed by previous studies ([24,31]), with the 1999 as the wettest year (1057.7 mm) and 1983 as the driest one (300.2 mm) (Supplementary Materials, Table S1). Figure 6a shows the yearly fluctuation in rainfall (black line) and the rainfall trend of the space state model obtained from Kalman filtering (red line). The rainfall trend is sinusoidal with a period of about 22 years; the two lows of the cycle are at 590 mm in 1985 and at 650 mm in 2008. The increase in rainfall trend (ΔP) over the period considered is in the order of 30–50 mm. Figure 6b shows the yearly ratio between cumulative pumping and rainfall, which helps filtering out the effects of different yearly precipitation amounts on pumping. The pumping-rainfall ratio has a minimum 0.2 in 1990, it steadily increases to the maximum of 0.42 in 2011, and then it stabilizes around 0.4.

The values and trends of climate extremes indexes are shown in Figure 7. Figure 7a shows the number of days in a year with intense (R10; p ≥ 10 mm/day) and Figure 7b with very intense rainfall (R20; p ≥ 20 mm/day). The R10 trend (Figure 7a) is sinusoidal with a 25 years long period. The period has a low of 14 days at the beginning of the cycle in 1985 and a second low of 16 days in 2010 at the end of the cycle. The R20 indicator trend (Figure 7b) is flat showing no increase in the extreme precipitation days during the period considered. The year with the greatest R20 is 1994 with 14 days of intense rain. In 1983, on the other hand, there are no R20 days. Figure 7c shows the wet days (WD) trend that is rather smooth with an average around 70 days. We observe a relative minimum around 68 days in the mid-eighties of the previous century. Figure 7d reports the daily precipitation intensity (SDII) and its trend. We observe a flat line for the trend showing no increase in the precipitation intensity during the period considered (1971–2017). The index has a maximum of 12.6 mm/day in 1999 and its minimum of 6.1 mm/day in 1983.

The potential evapotranspiration trend in the period 1971–2017 has increased from 740 mm/year to 850 mm/year (Figure S1 in Supplementary Materials). The real evapotranspiration trend has increased from 480 mm/year to 530 mm/year in the same period (Figure S2 in Supplementary Materials) with a minimum of 430 mm/year in 1984 and a maximum of 530 mm/year in 2015.

3.2. Land Use Change

Figure 8 shows land use change between 1976 and 2014, categorized according to three classes: urbanization, unchanged, and other variations. A small portion of the basin area (12.5%) has not undergone any variation over the years. Unchanged areas are those occupied by artificial water basins, beaches, coastal dunes, and pine forests protected by environmental restrictions, as well as the infrastructural settlements present since the 1970s. The urbanization process affects 13% of the entire basin area, and it includes land located near the road network. The Table in Figure 8 shows the percentages of total change for the five Corine Land Cover macro-classes along with the artificialization index. The two most significant variations refer to class 1, indicating a gradual urbanization process, and class 5, which includes all land reclamation canals, irrigation canals, rivers, and artificial basins that have expanded over time in order to manage water in the basin. The negative change of class 2 indicates that the arable crops, orchards, vineyards, and other cultivation systems have decreased over time and the agricultural land has been designated to other uses, above all urban.

With the aim to establish the effective influence of the urbanization on drainage in the basin, the correlation between rainfall and the equivalent amount of drained water at different lag time (in days) were calculated (Supplementary Materials, Table S2). Correlation values at lag time t = 0 (lag0 = same day) have a 0.3 average and maximum and minimum of 0.6 in 1973 and 0.1 in 1985, respectively. In order to include or exclude the influence of urbanization on drainage series, the datasets is split in two periods, one of minimum urbanization (1971–1981) and the second of maximum urbanization (2006–2016) and correlations are calculated to indirectly estimate how long rainfall takes to exit the system via the drainage system. The results are reported in

Table 1. Correlation coefficients between rainfall and drainage at different lag time (in days) for the period of minimum (1971–1981) and maximum urbanization (2006–2016).

Time Interval	r (lag0)	r (lag1)	r (lag2)	r (lag3)	r (lag4)	r (lag5)	r (lag6)
1971–1981	0.25	0.33	0.24	0.15	0.12	0.09	0.07
2006–2016	0.27	0.36	0.26	0.19	0.13	0.10	0.08

. The maximum correlation coefficient 0.33 in the period 1971–1981 is at a lag time 1 and for the period 2006–2016 it is 0.36 also at lag time 1. The values do not show any increase or decrease in the correlation values over time (see Supplementary Materials, Table S2).

3.3. Water Seepage

The average subsidence rates over the entire basin for different time intervals are listed in

Table 2. Average subsidence rate (mm/year) over the entire basin extrapolated from Teatini et al. [32] and Arpae dataset.

Time Interval	Subsidence Rate (mm/year)
1972–1977	45
1977–1982	5.5
1982–1986	5.0
1986–1992	6.0
1992–2002	5.0
2002–2006	5.8
2006–2011	5.7
2011–2017	4.4
	Total subsidence (mm)
1971–2017	482.2

. The increase of the seepage flux with on-going subsidence leads to the observed increasing trend of the drainage flux. The amount of water entering the system due to the vertical component of the seepage process, triggered by subsidence, is calculated by Equation (5) and using the land subsidence rates reported in Table 2 that are necessary to obtain the variation in vertical head difference through time. The amount of vertical seepage calculated for different values of K_v and sand unit thickness is listed in

Table 3. Total vertical seepage q_v (mm) calculated over the entire basin area for the considered time interval 1971–2017 considering different sand unit thickness (10–15 m, [15,27,43]) and vertical hydraulic conductivity (0.14–0.01 m/day) ranges. Note for the calculations the following parameters are used: the total area of Quinto Basin is equal to 9252 hectares; area affected by vertical seepage where the topography and water table are below sea level is equal to 2000 hectares; total Δh change through time as a result of subsidence is equal to = 0.53 m. K_v is assumed 1/3 of K_x (refer to paragraph 2.2.5. for K_x values).

qv (mm)		Kv (m/day)
		Arithmetic Mean	Geometric Mean	Harmonic Mean
		0.14	0.03	0.01
Sand unit thickness (m)	10	588	139	33
	11	534	126	30
	12	490	116	27
	14	420	99	23
	15	392	92	22
Mean qv		485	114	27

.

The results indicate that total equivalent amount of water due to vertical seepage (q_v) ranges from 485 to 27 mm considering maximum (arithmetic mean) and minimum (harmonic) K_v of the coastal aquifer, respectively. The geometric mean of K_v returns total seepage amount ranging from 139 to 92 mm that are comparable to the water budget calculation, with maximum value in the 1970s, when the highest subsidence rate were recorded (Table 2), and minimum value since the 2012, the period characterized by the lowest subsidence rate.

3.4. Water Budget

Using the water budget Equations (6) and (7), we estimated ΔQ_v, the increase of vertical seepage during three sub-periods where the trends of the water budget components have distinct trends (increasing or decreasing); the sub-periods considered are 1971–1986, 1986–2000, and 2000–2017. The cumulative ΔQ_v for the whole period considered (1971–2017) can be obtained by adding up the sub-period fractions (Table S3 in Supplementary Materials). The ΔQ_v = 125 mm/year estimate for the whole period is obtained from the water budget Equation (7) using ΔP = 40 mm/year, ΔETR = 40 mm/year, and ΔPu = 130 mm/year.

Note that in 1971, by using the trend components of the water budget obtained from Kalman filtering, we estimate P = 645 mm/year, ETR = 485 mm/year, and Pu = 160 mm/year, so that applying Equation (6), we obtain Q_v = 0 mm/year.

4. Discussion

The drainage data time series show the actual workload that the water pumping station constantly faces throughout the year. The ascending trend from 1971 to 2017 (Figure 5a) demonstrates the increasing difficulty of the pumping station to efficiently and effectively contrast the processes that convey water into the basin.

The drainage rates in the Quinto Basin vary widely over the 46 years examined (Figure 5a) and, specifically, they tend to increase sharply with a cumulative total increase of 150 mm from 1988 to 2017. This information confirms that the water pumping station has been working more to cope with a larger amount of water to be drained. If we consider the number of drainage days in a year (Figure 5b), it is apparent how drainage is not concentrated anymore in a few days, but, today, it is distributed throughout the year. The maximum daily pumping per year (Figure 5c) has only increased slightly (from 6.2 to 6.6 mm/day) given that it is limited by the total lifting capability of the pumping station. The increase in annual drainage, therefore, is accommodated by draining throughout the year (with an average of 306 days per year in the last 17 years). This is evidence that the system is under more stress compared to the past and that there is more water to be lifted.

It is well known that a heavy drainage system can be responsible for saltwater intrusion in coastal regions [17,61,62,63]. In our case, low hydraulic heads in the coastal aquifers due to drainage cause an upward seepage of saltwater from the bottom of the aquifer as demonstrated from our calculations in Table 3. As suggested by our modelling of Equation (5), the upward seepage seems to origin at the bottom of the sandy layer, which has a total thickness ranging from 10 m, on the eastern part of the study area, to about 15 m inland, where the prodelta has a limited thickness or is absent ([27,40,41,42,43], Figure 3). The ongoing subsidence through time and the consequent increase in drainage have created an increased hydraulic vertical gradient resulting in about 139–92 mm of total equivalent amount of water due to vertical seepage.

It is possible, however, that the increase in vertical gradient induced by land subsidence causes the connate saline groundwater, previously preserved in the deepest portion of the aquifer in stable density stratification, within low permeable layers or in hydraulic traps, to be mobilized and slowly transported upward [18]. As a result of this, the shallow coastal aquifer of Ravenna is completely compromised by salinity, with shallow and thin freshwater lenses (2–3 m) below the coastal dune belt and palaeodunes [7,15,19,27] and hypersaline water in the deepest portion of the aquifer [18].

It is therefore necessary to determine which phenomena cumulatively allowed the entry of the extra water that drainage data indicate. Considering the factors affecting the increase in pumping, climate datasets, and related climate extreme indexes do not highlight any clear trend (Figure 7). The ratio between drainage and rainfall (Supplementary Materials, Table S1, Figure 6b) has increased greatly in recent years even if we consider years of equal rainfall. According to the calculations performed over the whole dataset, the ratio is on average equal to 0.3. The trend indicates an increase of 67% since 1971, with an absolute maximum value of 0.8 in 2011 and a minimum value of 0.08 in 1975. Taking, for example, two years (1972 and 1996) with similar annual precipitation (about 950 mm/year), we observe a 244.3 mm/year annual drainage in 1972, which is much lower than the 433.7 mm/year drained in 1996. Obviously, this type of variation strictly depends also on the annual distribution of precipitation. We observe, however, that in the last twenty years, despite not having increased annual rainfall, the drainage rates are constantly increasing (Figure 5 and Figure 6). This suggests that there must be other factors affecting the amount of water entering the basin.

The same conclusion can be reached if runoff is indirectly estimated by subtracting the annual real evapotranspiration (ETR) from the annual precipitation (P) (Supplementary Materials, Table S1). Theoretically, the drainage infrastructure should completely drain the runoff water after the rainfall events. The average runoff estimated over 46 years is 153 mm, varying from a minimum of 10 mm in 1983 to a maximum of 423 mm in 1999. The annual drainage differs from runoff (Figure 5 and Table S1 in Supplementary Materials). In recent years, the annual drainage is higher than estimated runoff; the trend is visible starting from 2005 when, for 255 mm of runoff, 373 mm of equivalent water are drained. The large mismatch between the two datasets is also observed in 2011, when 310 mm of equivalent water are drained against 27 mm of runoff.

Since drainage is intrinsically linked to precipitation events and their distribution, the analysis of climate data time series must be considered. In general, all climate indexes (Figure 7) that hypothetically could highlight climate extremes have not shown clear trends. Climate extremes, therefore, cannot explain the increase in pumping rates.

These results are also in agreement with the ISPRA report [53], which highlights how the analysis of variations and trends of precipitation extremes in Italy do not clearly indicate any increase or decrease, demonstrating high spatial variability among stations and no clear climate signal at national level.

The Quinto Basin, like all areas surrounding Ravenna, has undergone an important urbanization during the years of the post “economic boom” (1970–1980s) [64,65] (Figure 8). The major factors linked to urbanization are the extension of the Classe coastal area, the construction of a regional highway, gravel pits exploitation, a gradual enlargement of farms and commercial infrastructures, etc. The analysis of land use change shows an index of artificialization (I_art) in line with values observed for the main cities of the Emilia-Romagna region [51,66]. Urbanization and shift to impervious surfaces (e.g., buildings, paved roads) causes an increase in runoff by reducing precipitation infiltration and shortening flow routing time [67,68]. The rainwater that cannot infiltrate into the soil flows at the surface and is collected quickly by the drainage system. The Quinto Basin water pumping station is in one of the most depressed areas of the entire study area [31] and it is the endpoint of a dense network of drains and ditches. Most of the runoff water, therefore, is forced to flow towards these drains and, finally, to the station pumps. In a scarcely urbanized area, water would exit the system with a certain delay after the rainy event, respecting the time required for the processes of infiltration, percolation, storage, and water retention [69], depending on the characteristics of the vadose zone above the aquifer and aquifer parameters. In this case, the largest correlation values between precipitation and drained water would be associated with high lag times (i.e., lag3, lag4, lag5, lag6). On the contrary, in a strongly urbanized system with large impervious areas, the water does not easily infiltrate and tends to flow on the surface as runoff, being collected from the hydrographic and drainage system almost instantly [70]. In these specific situations, the correlation values between the amount of rainfall (mm) and the equivalent drained water (mm) would be higher at the lowest lag times (i.e., lag0, lag1, lag2). Table 1 and Table S2 (Supplementary Materials) show the correlation values between rainfall and drainage at different lag times. If the artificialization in the basin had increased in such a way as to play an important role in the amount of water flowing to the pumping station, then the correlation values should be maximum at lag time n + 0 (lag0); specifically, this trend should be expected in the years of maximum artificialization, from 2006 to 2016 (Table 1). The correlations for the different lag times in Table 1 and Table S2 (in Supplementary Materials), however, do not follow a well-defined trend and the highest values, in most cases, do not correspond to the lag0. In general, the response of pumping with respect to rainfall events does not seem to follow a trend related to artificialization, probably because urbanization is not uniform and does not trigger an unequivocal response. Urbanization, therefore, does not seem to be the most important factor for the increase in drainage. The lack of correlation between drainage and urbanization is also due to the drainage management operated by the Land Reclamation Consortium, which proceeds, if heavy rainfall events are forecasted, to unload the channels and increase the drainage rate in advance in order to manage any flooding emergency.

The drainage data show an increasing trend, which means that with time the yearly drainage flux (amount) increases. The input from the nearby reclamation basins into which the territory of Ravenna is divided can be excluded, because each basin refers to a water pump station, which lifts and discharges the drained water to the sea by a single drainage channel; from this point the water is not used anymore for other purposes (i.e., irrigation, etc.; [30]). Between these basins, there is not water exchange and flux, rather they act as independent polder basins.

The geology can also exclude the input from upstream aquifers since the shallow coastal aquifer is confined at the bottom by the Flandrian continental clay basement, and in the West part by fine sediment of back-barrier/marsh and lagoon deposits (Figure 3) and recharge only occurs by rainfall infiltration.

It should be pointed out that the shallow coastal aquifer is used neither for drinking purposes nor for irrigation and industrial activities. For drinking needs, Ravenna is served by surface water from the main rivers and partially by the Ridracoli water reservoir. Irrigation water is also derived from rivers and distributed through the dense network of channels managed by the Land Reclamation Consortium. Moreover, the increase in drainage flux cannot be explained by the abrupt decrease is groundwater abstraction that took place thereafter the construction of a new industrial aqueduct (in the late 1970s) and the activation (at the end of the 1980s) of the public aqueduct which produced a generalized curtail of the groundwater withdrawals, from the deepest confined aquifers [32,33,34], but that did not affect the shallow aquifer.

An explanation of this increase in drainage flux, could be the extra vertical seepage from the bottom of the aquifer towards the surface, which is caused by the increase in vertical hydraulic gradient driven by land subsidence. As a result, the amount of vertical seepage keeps increasing with the on-going land subsidence since the vertical head gradient in the polder increases every year due to the land subsidence, resulting in a seepage flux increase linear to the head difference increase. The increase of the seepage flux with on-going subsidence, therefore, leads to the observed increasing trend of the drainage flux.

Based on water budget considerations (see Table S3 in Supplementary Materials), we have shown that starting in 1971 with a negligible Q_v (~ 0 mm/year), the ΔQ_v in the period 1971–2017 has been 125 mm/year. This figure is also confirmed independently by our calculations based on Darcy law (Equation (5)), which highlight 139–92 mm of vertical seepage for the same period considering mean vertical hydraulic conductivity (K_v) of 0.03 m/day. The seepage values are very sensitive to K_v values, so it is important to have an accurate estimate of the hydrogeological parameters of the studied aquifer. There are several studies that discuss the best assumption to calculate the effective hydraulic conductivity in the vertical direction based on the aquifer characteristics [71,72,73,74]. Freeze [75] showed that, in the case where a heterogeneous aquifer system is composed of alternating units of different lithologies (for example sand and clay beds), the effective hydraulic conductivity parallel to the layering of sediments (usually the horizontal directions) is commonly equal to the arithmetic mean of the hydraulic conductivity of the individual beds weighted by the thickness of each [75]. The effective vertical hydraulic conductivity is calculated as being equal to the harmonic mean of the hydraulic conductivity of the individual beds. Madden [76] modelled two-dimensional heterogeneous media and concluded that the geometric mean can be used to predict conductivity of heterogeneous mixtures of materials, and can be used, under certain flow conditions, as representative of the overall conductivity of the medium. In our case we used the three types of mean and the geometric averaging that seem to work very well if compared against water budget calculation.

In the coastal aquifer of Ravenna and Ferrara, Antonellini et al. [1] and Giambastiani et al. [15,18] highlighted the presence of vertical hydraulic gradients due to the drainage system, which control the hydrodynamics of these areas, forcing groundwater to flow from the bottom toward the top of the aquifer. In some cases, this upward flow is contrasted by small lenses of silt and clay that confine the aquifer at the top and prevent the complete salinization of the shallow aquifer.

The ETP has been increasing steadily since 1971 due to a constant increase in average temperature (Figure S1 in Supplementary Materials). The ETR, on the other hand resents of both temperature and precipitation trends but, overall, it has also been increasing (Figure S2 in Supplementary Materials). We assume that the vertical seepage is intercepted by the drainage network that must guarantee agricultural activities and keep the water table below the crop root zone, so that we do not include it in the computation of ETR. Despite a possible effect of the variations in ETP and ETR and the fluctuations of P and ETR on the computation of the vertical seepage, today Q_v adds about 125 mm/year of water to the input flux in the Quinto Basin with respect to 1971, which is the same order of magnitude to what obtained in the vertical flux computation by considering the increase in hydraulic head difference driven by land subsidence (~ 0.5 m). The increase of the precipitation/pumping ratio (Pu/P) over time (Figure 3b) reflects the fact that the pumps need to lift not only the surface runoff caused by precipitation, but also the extra vertical seepage input. This extra input consists of saltwater, which contributes to the increasing trend of salinization observed in the area [1,15,18,19,77]. The increase in pumping and vertical seepage is most apparent starting from 1986. In the period 1971–1986, in fact, there was a strong decrease in yearly rainfall and ETR, which probably masked the effect of the vertical seepage, which seems to be delayed with respect the land subsidence rates. The highest subsidence rate recorded in the Ravenna territory, in fact, were observed during the 70s of the past century (maximum values of 110 mm/y; [32]). Those were the years of the ‘economic boom’, characterized by a strong extraction of geofluids from the subsurface (water and methane gas). The subsidence process today is still important and diffuse but has certainly less intensity than in the past. Subsidence rates in the last twenty years, in fact, have stabilized [55,78,79]. Today, the average rate of total land subsidence in Ravenna is about 4–5 mm/year with the highest values recorded along the coast of the Quinto Basin (values ranging from 17 to 10 mm/year) [35].

5. Conclusions

The objective of this work was to identify and quantify the relative contributions of the main factors (climate extremes of precipitation, land use change, land subsidence, and seepage processes) affecting water drainage long time series (1971–2017) in a polder basin (Quinto Basin), within the low lying coastal area of Ravenna (Italy). This basin is also representative for many other polder basins of the Northern Adriatic coast and for other basins in low coastal plains worldwide (The Netherlands, Bangladesh, China, etc.) affected by saltwater intrusion problem.

The analysis of climate extremes does not follow well-defined trends and cannot justify the increasing drainage trend. Rainfall within the Quinto Basin has not changed significantly in the last fifty years in terms of intensity and extreme events.

The analysis of land use change indicates that the loss of infiltrating soil capacity, as a result of urbanization and soil sealing, took place at such an intensity that only slightly, if anything, affected the drainage rates and without following a well-defined trend.

The amount of extra water entering the basin is due to vertical seepage from the bottom of the upper sandy portion of the surface aquifer, which is caused by an increase in hydraulic head difference driven by the land subsidence rate. This extra vertical seepage justifies the increase in drainage recorded between 1971 and 2017. The increase in drainage flux is much lesser influenced by the process of urbanization.

The drainage in the Quinto Basin is currently less efficient in contrasting water inflow occurring in the study area compared to the past. Nowadays it is necessary to drain more continuously—almost every day—than in the past to compensate for the effects of the former and current land subsidence.

Since, among all factors causing costal aquifer salinization, water drainage plays an important role in lowering the hydraulic head and favouring saltwater up-coning (seepage), the area will experience highly negative effects on coastal groundwater resources and freshwater availability as long as land subsidence will continue in the future.

This study also shows a way to quantify the relative contribution of different processes to the damage caused by land subsidence to a coastal area. The damage, in fact, has some economic costs and these include operation costs (electrical power, infrastructure upgrading, etc.), hydrologic risk increase (flooding of low-lying areas), salinization of ground and surface water, as well as the damage to ecosystems that are threatened by the changed hydrologic conditions. Further knowledge of these contributions could help to better devise fair compensation plans to be charged to those activities responsible for the damage.