Groundwater-Surface Water Interaction in the Nera River Basin (Central Italy): New Insights after the 2016 Seismic Sequence

: The highest part of the Nera River basin (Central Italy) hosts signiﬁcant water resources for drinking, hydroelectric, and aquaculture purposes. The river is fed by fractured large carbonate aquifers interconnected by Jurassic and Quaternary normal faults in an area characterized by high seismicity. The 30 October 2016, seismic sequence in Central Italy produced an abrupt increase in river discharge, which lasted for several months. The analysis of the recession curves well documented the processes occurring within the basal aquifer feeding the Nera River. In detail, a straight line has described the river discharge during the two years after the 2016 seismic sequence, indicating that a turbulent ﬂow characterized the emptying process of the hydrogeological system. A permeability enhancement of the aquifer feeding the Nera River—due to cleaning of fractures and the co-seismic fracturing in the recharge area—coupled with an increase in groundwater ﬂow velocity can explain this process. The most recent recession curves (2019 and 2020 periods) ﬁt very well with the pre-seismic ones, indicating that after two years from the mainshock, the recession process recovered to the same pre-earthquake conditions (laminar ﬂow). This behavior makes the hydrogeological system less vulnerable to prolonged droughts, the frequency and length of which are increasingly affecting the Apennine area of Central Italy. model describes both and 2020 recession curves as a result of the best-ﬁt analysis, indicating that the hydrogeological system feeding the river empties with a Darcian laminar ﬂow. The two curves show similar a recession coefﬁcient α E , − 2.1 × 10 − 3 day − 1 and − 2.3 × 10 − 3 day − 1 for 2019 and 2020, respectively.


Introduction
The interaction between groundwater (GW) and surface water (SW) influences river water quantity and quality. The understanding of the processes and dynamics of GW-SW interactions are fundamental for the accurate assessment, integrated management, and environmental protection of water resources [1][2][3][4][5]. Although the anthropic pressure on many river basins is increasing [6], climate change is negatively impacting the river discharge, especially in regions characterized by a reduction of snow and rainfall during the recharge periods [7].
The Nera River in Central Italy represents a hydrogeological system where SWs are mainly provided by GWs, thanks to a set of permanent linear springs, the water of which comes from large fractured and karstified basal carbonate aquifers [8][9][10][11]. Aquifers hosted in the upper part of the Nera River supply water to a multipurpose system (drinking water, hydropower energy production, and fish farming). The study area is located in a region affected by a decrease in rainfall during the recharge period, which occurs from autumn to early spring [12][13][14][15][16][17]. A recent review of rainfall trends published by Caporali et al. [18] reveals a more pronounced negative trend in winter periods in Central Italy than in Northern Italy. This general trend is coupled with the increase in length and frequency of drought periods in the last two decades [15,19,20]. Moreover, as reported by Diodato and Bellocchi [21], the number of snowy days declined in peninsular Italy from the end of the Little Ice Age (LIA) and, markedly, after the 1940s. Since snowmelt and rainfall affects the

Geological Structural and Hydrogeological Setting
The area investigated in this study is located in the mountainous region of the Sibillini Mountains (northern Apennines) in Central Italy. The stratigraphic sequence of this part of the Apennine consists of Meso-Cenozoic carbonate multilayer formations (Fms) of the Umbria-Marche succession and siliciclastic foredeep deposits of the Laga Fm [32]. Figure 1a shows the hydrogeological map of the upper part of the Nera River basin. It is based on a new geological survey, which checked and updated the geological map published by [32]; rocks are described considering the stratigraphic relationships and focusing on the hydrogeological properties, e.g., grouping them in hydrogeological complexes [33][34][35][36].
The sequence starts with Upper Triassic evaporites, including anhydrides and dolomites (evaporites Fm). These rocks are not exposed in the study area and are characterized by a general low permeability (Evaporites Complex, EC). Upward, the sequence continues with shallow-water carbonates (Calcare Massiccio Fm) with a variable thickness of 600 to 700 m. In the Jurassic structural lows, a thick prevalent pelagic sequence was deposited. Micritic limestones represented the basal deposition of structural lows with nodules of chert (Corniola Fm). Being high fractured/karstified limestones, the Calcare Massiccio and Corniola host the basal aquifer with a maximum thickness of ca. 900 m (Basal Limestones Complex, BLC). The pelagic sequence continues with Calcareous Siliceous Marly Complex (CSMC, including Rosso Ammonitico, Marne del Serrone, Calcari a Posidonia, and Calcari Diasprigni Fms), characterized by low relative permeability. The subsequent deposition of Maiolica Fm occurred on both structural highs and lows, with a variable thickness from 150 to a maximum of 400 m, respectively (Maiolica Complex, MAC). Stratified and fractured limestones characterize the hydrogeological complex, with high relative permeability. Due to high permeability, it represents the middle aquifer. The sequence continues upward with the Marne a Fucoidi Fm, a marly rock with low relative permeability (Marne a Fucoidi Complex, MFC). Marly-limestones (about 400 m thick) were then deposited above the Marne a Fucoidi Fm; it includes the Scaglia Calcarea Complex (SCC, Scaglia Rossa and Scaglia Bianca), characterized by a moderate relative permeability. Above the SCC, Scaglia Variegata and Scaglia Cinerea Fms are deposited, composed of marls and marly limestones with low permeability (Calcareous Marly Complex, CMC). The subsequent Miocene marly units (Schlier and Bisciaro) and siliciclastic Laga Fm are mainly exposed in the eastern sector of the studied area. Most of these Fms are characterized by moderate or low permeability (Terrigenous Units Complex, TUC).
Formations belonging to the hydrogeological complexes were involved in distinct tectonic phases since the Jurassic. First, extensional tectonics was associated with the thinning of the Adriatic continental margin during Late Jurassic [32], involving the Basal Limestones Complex. Since the late Miocene, the whole limestone multilayer was affected by a compressional phase that produced important shortening of this part of the Apennines with a typical fold and thrust belt [37,38]. The Sibillini Mountains thrust (MST, Figure 1a) represented the main regional compressional structure [39], marking the tectonic boundary between the Mesozoic-Paleogene limestone sequence at its hanging-wall and the Late Miocene-Early Pliocene Laga sequence (Flysch della Laga Fm, [40]) at its footwall. A more internal compressional structure is represented by the Pizzo Tre Vescovi thrust (PTV, Figure 1a).
The compressional structures were subsequently crosscut by the NNW-SSE trending normal faults due to the prevalent extensional Quaternary tectonic phase [41,42]. These faults are still active and responsible for historical seismicity and represent preferential drainage paths of groundwater in the basal aquifer, which usually hinder the transversal groundwater exchanges [35,43].    [44]).

Seismic Sequence and Co-Seismic Effects
The study area was struck by a long and intensive seismic sequence during 2016 and 2017 (Figure 1a), characterized by the occurrence of nine mainshocks with moment magnitude, Mw > 5.0 ( Table 1). The seismic sequence was due to a system of mainly NNW-SSE striking and WSW dipping normal faults, as shown by the focal mechanism of the main shocks, the distribution of the aftershocks, and significant co-seismic ruptures [46][47][48][49][50][51][52][53]. The main seismogenic normal faults activated during the seismic sequence were the Vettore Mt.-Bove Mt. fault systems (VBF) to the north and the Laga (LAF) fault systems to the south [55]. This study is focused on the northern sector affected by the VBF segment.
In particular, the southern portion of the VBF and the northern portion of the LAF (Figure 1a) were activated during the 24 August Mw 6.0 Accumoli event, which enucleated at their overlapping area, at a depth of 8 km [51,56]. The 26 October Mw 5.9 Visso event enucleated at a depth of 4 km [57] and activated the northern part of the VBF, while during the 30 October Mw 6.5 Norcia event, the VBF was activated to the south and center. The average co-seismic throw of the entire sequence is~0.3 m [52,53] and the activation of the VBF during the largest events produced important surface ruptures developed along major and minor fault segments as shown in Figure 1a [51,53]. Figure 1b shows a cross section roughly oriented WE with co-seismic displacements produced by the October 30, 2016, seismic sequence as taken from Valerio et al. [45]. Geodetic measurements, using the Global Navigation Satellite System (GNSS) technique together with the available surface deformation field obtained on the basis of the DInSAR technique, registered the co-seismic and immediately post-seismic deformation of the two major October shocks surrounding the VBF [44]. One of the GNSS stations ( Figure 1a) recorded westward horizontal displacements of 419 mm and subsidence of 707 mm (with 95% confidence errors), with a total off-fault vertical displacement between footwall and hanging-wall blocks of 736 mm in correspondence of the northern sector of the VBF.

Meteorological Data and Techniques of Analysis
The study of river discharge fed by regional aquifers requires reliable meteorological data. According to Cambi et al. [13], as for other mountain regions in Central Italy, often few data are available to define the hydrogeological scheme. Figure 2 shows all the thermopluviometric stations located within the Nera River catchment at Visso section. As for other locations in the Mediterranean region, the upper part of the Nera River catchment has had changes in the data acquisition network over time; the overall monitoring network has moved from the Servizio Idrografico e Mareografico Nazionale (SIMN) to the Civil Protection Agency of Marche Region (SIRMIP) during the beginning of the 2000s, resulting in lack of data often coupled with the relocation of some stations.
To check the consistency of rainfall data, the double-mass analysis was used [57]. This method compares cumulated rainfall data of a single station with that of other stations in the area, allowing the individuation of abrupt deviations produced by instrument malfunctions and/or relocations. As recently reported by Caloiero et al. [19], lack of spatially distributed data can be overcome due to the increased availability of meteorological sat- Rainfall monthly data are downloaded from https://disc.gsfc.nasa.gov/datasets/GPM_3IMERGM_06/summary?keywords=%22IMER As for other locations in the Mediterranean region, the upper part of the Nera River catchment has had changes in the data acquisition network over time; the overall monitoring network has moved from the Servizio Idrografico e Mareografico Nazionale (SIMN) to the Civil Protection Agency of Marche Region (SIRMIP) during the beginning of the 2000s, resulting in lack of data often coupled with the relocation of some stations. To check the consistency of rainfall data, the double-mass analysis was used [57]. This method compares cumulated rainfall data of a single station with that of other stations in the area, allowing the individuation of abrupt deviations produced by instrument malfunctions and/or relocations. As recently reported by Caloiero et al. [19], lack of spatially distributed data can be overcome due to the increased availability of meteorological satellites. The Global Precipitation Measurement (GPM) provides monthly rainfall data on an approximately 10 × 10 km grid across the globe from April 2014 to the present. The latest release of the Integrated Multi-satellitE Retrievals for Global Precipitation Measurement product (IMERG, Version 6B-Final) fused GPM estimates with those collected by the TRMM satellite (Tropical Rainfall Measuring Mission), operating during 2000-2015 (https://gpm.nasa.gov/data/imerg, accessed on 25 May 2021). The Nera River catchment is included in a cell having the following lat-long coordinates: 42.45-13.05; and 42. 55-13.15. Rainfall monthly data are downloaded from https://disc.gsfc.nasa.gov/datasets/GPM_3IMERGM_06/summary? keywords=%22IMERG%20final%22 (accessed on 25 May 2021). Due to the complex topography in mountain regions, checking the performance of satellite data is necessary to achieve analyses that are as realistic as possible [58][59][60]. To quantify the performance of GPM-IMERG data, Pearson's Correlation Coefficient (CC, Equation (1)), Nash-Sutcliffe efficiency (NSE, Equation (2)), Root Mean Square Error (RMSE, Equation (3)), and the relative bias (rBias, Equation (4)) are used; the latter reflects the systematic bias between the satellite and the ground observations [61].
where: S i = satellite precipitation estimate (mm/month) G i = rain gauge observation (mm/month) S = mean satellite precipitation estimate (mm/month) G = mean gauge precipitation observation (mm/month) Reliable monthly rainfall datasets are useful to check the frequency and length of droughts, which affect the groundwater recharge of hydrogeological systems. Worldwide, the Standardized Precipitation Index (SPI, [62]) is used as an effective drought index for detecting and characterizing meteorological droughts. According to WMO [63], for the computation at least 20-30 years of monthly rainfall data are needed. In the SPI computation the data are fitted to a gamma probability distribution and then transformed into a normal distribution. For a given data time series X i as X 1 , X 2 . . . .X n , the SPI i is defined by Equation (5).
where X is the arithmetic mean of rainfall and S x is the standard deviation. According to McKee et al. [62], wet periods occur when SPI values are higher than 1 while drought periods occur when SPI is lower than −1. SPI equal to zero implies that there is no deviation from the mean. Table 2 shows the classification for the SPI values. Monthly precipitation (P) and temperature (T) have been used to compute the monthly evapotranspiration (ETR) using the Thornthwaite and Mather [64] method. This method is still useful when no data concerning radiation, air humidity, etc.-necessary to apply more modern methods-are available, as occurred for the study area [59,65]. The computation of ETR was carried out by using the software developed byČadro [66]. Surplus monthly data (S = P-ETR) were then computed to individuate the period with no recharge (P = ETR), i.e., the river discharge is fed only by the groundwater stored in the hydrological system, which empties as the depletion phase proceeds.

Discharge Data
Daily stream levels (H) of the Nera River were collected from the SIRMIP on-line monitoring network on the I1 river section (Figure 2), which also includes its main tributary (Ussita stream, I2 stream gauge). By using the rating curve published by SIRMIP (http: //app.protezionecivile.marche.it/sol/indexjs.sol?lang=it, accessed on 25 May 2021), daily stream discharge values (Q) were computed from January of 2016. No substantial anthropic modification occurred in the catchment in the last decades that could have influenced SW-GW interactions through the river system. Apart from the construction of a fluent hydroelectric plant close to the P1 station ( Figure 2), a drinking water intake was built near Castelsantangelo sul Nera village in the 1980s (San Chiodo spring, Figure 1). As the water is piped out of the Nera River catchment, the amount of drinking water withdrawn was added to the river discharge data; it corresponds to 0.150 m 3 /s from January of 2016 to February of 2018 and 0.200 m 3 /s from March of 2018 to today [67]. As reported by Di Matteo et al. [24], in addition to recent monitoring data, some historical discharges from the SIMN monitoring network (I3 and I4 in Figure 2) are also available (1928-1943 period).

Models for River Recession Curves
The study of the discharge during periods with no recharge (recession curve) provides insightful hydrologic information about the hydrogeological properties of the system feeding the river, at least in terms of average or equivalent values. Several conceptual models for recession curves have been developed to describe the discharge during the recession period. In the hydrogeological practice the exponential formula (EF) is extensively used (Equation (6), [68,69]). Although EF was developed for porous media (Darcian laminar flow), it has also been widely applied to fractured aquifers characterized by low hydraulic gradients and velocity [70][71][72]. When these conditions are not met (turbulent flow), the straight-line equation (SL) can describe the recession process in fractured and karst aquifers (Equation (7), [73,74]). The rate of change of the discharge (RCD) is computed making the first derivative of Equations (6) and (7), obtaining Equations (8) and (9) for EF and SL, respectively, According to Rorabaugh [75], and Kovacs and Perrochet [76], α E is linked to the aquifer characteristics by Equation (10).

Rainfall-Water Surplus Analysis
The recharge of aquifers is affected by prolonged drought periods, the evaluation of number and frequency of which require reliable rainfall data. Figure 3 shows the results of the double mass analysis; data of Ussita station (PT-2) are compared with those of ENDESA (P1) and M. Prata (PT-3) stations. The cumulative rainfall of PT-2 and P1-located at similar altitude-is aligned along a straight line with no abrupt deviations. On the contrary, an abrupt change between the years 2008 and 2010 is observed by comparing the cumulated rainfall data of PT-2 with PT-3. The latter shows lower cumulative rainfall despite being located 1050 m higher than PT-1. The double-mass analysis results for the 2002-2020 period found the Ussita station to be the most reliable in the study area. The recharge of aquifers is affected by prolonged drought periods, the evaluation of number and frequency of which require reliable rainfall data. Figure 3 shows the results of the double mass analysis; data of Ussita station (PT-2) are compared with those of ENDESA (P1) and M. Prata (PT-3) stations. The cumulative rainfall of PT-2 and P1-located at similar altitude-is aligned along a straight line with no abrupt deviations. On the contrary, an abrupt change between the years 2008 and 2010 is observed by comparing the cumulated rainfall data of PT-2 with PT-3. The latter shows lower cumulative rainfall despite being located 1050 m higher than PT-1. The double-mass analysis results for the 2002-2020 period found the Ussita station to be the most reliable in the study area. As illustrated in Section 2.1, the computation of SPI requires at least 20-30 years of monthly rainfall data [62]. The Ussita station does not meet this requirement, i.e., it is not suitable for SPI analysis. Among the available stations (Figure 2), only the Bolognola station could be useful for this analysis (P3 and PT-5). Unfortunately, in 2007, the historical rain gauge (P3, SIMN) was removed, and the new station was located at a higher elevation than the historical one (PT-5, SIRMIP), i.e., the rain gauge was moved from 1070 m a.s.l. to 1370 m a.s.l. Under these conditions, the aggregation of the two-time periods (1951-2007 and 2007-2020) was not successful; therefore, the data were not used for the SPI computation. In this framework, satellite data can be a valid alternative for this analysis, especially in data-scarce regions. The results of the application of some accuracy indices show that the correlation among monthly GPM-IMERG satellite data and rainfall at the Ussita gauge reaches a CC value equal to 0.76, a NSE of 0.99, a RMSE of 4.5 mm, and a rBias of about 24%, confirming that satellite data tend to slightly overestimate rainfall data as already recorded by Navarro et al. [59] along the Apennine chain of the Italian peninsula. The results of statistical evaluations indicate that-within the limits of use of satellite data in mountain areas (e.g., effects of topography, etc.)-the performance of GPM-IMERG allows the use of this dataset to compute the SPI during the 2000-2020 period. Similar results were recently obtained by Caloiero et al. [18], who used the GPM-IMERG data to calculate SPI values for the Italian peninsula. Figure 4a   As illustrated in Section 2.1, the computation of SPI requires at least 20-30 years of monthly rainfall data [63]. The Ussita station does not meet this requirement, i.e., it is not suitable for SPI analysis. Among the available stations (Figure 2), only the Bolognola station could be useful for this analysis (P3 and PT-5). Unfortunately, in 2007, the historical rain gauge (P3, SIMN) was removed, and the new station was located at a higher elevation than the historical one (PT-5, SIRMIP), i.e., the rain gauge was moved from 1070 m a.s.l. to 1370 m a.s.l. Under these conditions, the aggregation of the two-time periods  and 2007-2020) was not successful; therefore, the data were not used for the SPI computation. In this framework, satellite data can be a valid alternative for this analysis, especially in data-scarce regions. The results of the application of some accuracy indices show that the correlation among monthly GPM-IMERG satellite data and rainfall at the Ussita gauge reaches a CC value equal to 0.76, a NSE of 0.99, a RMSE of 4.5 mm, and a rBias of about 24%, confirming that satellite data tend to slightly overestimate rainfall data as already recorded by Navarro et al. [60] along the Apennine chain of the Italian peninsula. The results of statistical evaluations indicate that-within the limits of use of satellite data in mountain areas (e.g., effects of topography, etc.)-the performance of GPM-IMERG allows the use of this dataset to compute the SPI during the 2000-2020 period. Similar results were recently obtained by Caloiero et al. [19], who used the GPM-IMERG data to calculate SPI values for the Italian peninsula. Figure 4a shows the SPI-12 computed over the last two decades; prolonged drought events are detected in 2001-2003, 2007, 2012, and 2018. Three moderate to near normal wet periods are interposed between these events. Similar findings-even if with slightly different intensities and with some difference for 2018-were obtained by Valigi et al. [25] for the Pescara di Arquata Spring recharge area, located about 20 km southeast of the upper part of the Nera River catchment.
Arquata Spring recharge area, located about 20 km southeast of the upper part of the Nera River catchment. Figure 4b shows the cumulated water surplus (S) during the recharge periods (September-April); weather data of the Ussita station are used for the computation. During the observation period, the mean S was about 550 mm/8-months with the highest values from 2013-14 to 2017-18, confirming the results of SPI analysis. On the contrary, in the 2018-19 and 2019-20 periods, S values were 100-200 mm lower than the mean, corresponding to the moderate/severe dry periods highlighted with the SPI-12 (Figure 4a).

River Hydrograph Analysis
The flow regime of the Nera River was studied considering the available daily discharge data. According to the SIRMP monitoring system, reliable data were recorded during the 2016-2021 period ( Figure 5). These data include both pre-and post-seismic phases.
The October 30, 2016, earthquake (n. 5 in Figure 5) produced the abrupt discharge rising observed during November-January of 2016 with an increase of about 4.8 m 3 /s. A small peak flow rate had already been recorded on August 24, 2016, connected with the Accumoli-Norcia earthquakes (ns. 1-2 in Figure 5). Figure 5 also shows the monthly water surplus, computed according to the Thornthwaite and Mather [64] method, and daily discharge data of the Ussita River. The discharge hydrograph of this sub-catchment overlaps the Nera River one with some differences during the late August-September periods of 2017 and 2018. In these two periods, the recession phase of the Ussita River continued while that of the Nera River stopped after receiving some recharge that did not occur in the Ussita catchment.
Considering periods with no aquifer recharge (water surplus S = 0 in Figure 5) the pre-and post-seismic recession curves were collected. After the individuation of months with S = 0, a deep analysis was carried out to separate consecutive days characterized by a decrease in discharge values, avoiding days with discharge increases in response to the previous effective rainfalls. These curves help to investigate the relationship between aquifer hydrogeological properties and groundwater discharge to the river. The results of the analysis of the 2019 and 2020 depletion phases are presented in Figure 6, moving to the discussion section for the comparison with previous studies published in the study area in the two years following the earthquake. The EF recession model describes both 2019 and 2020 recession curves as a result of the best-fit analysis, indicating that the hydrogeological system feeding the river empties with a Darcian laminar flow. The two curves show similar a recession coefficient α, −2.1 × 10 −3 day -1 and −2.3 × 10 −3 day −1 for 2019 and 2020, respectively.

River Hydrograph Analysis
The flow regime of the Nera River was studied considering the available daily discharge data. According to the SIRMP monitoring system, reliable data were recorded during the 2016-2021 period ( Figure 5). These data include both pre-and post-seismic phases.
The 30 October 2016, earthquake (n. 5 in Figure 5) produced the abrupt discharge rising observed during November-January of 2016 with an increase of about 4.8 m 3 /s. A small peak flow rate had already been recorded on 24 August 2016, connected with the Accumoli-Norcia earthquakes (ns. 1-2 in Figure 5). Figure 5 also shows the monthly water surplus, computed according to the Thornthwaite and Mather [64] method, and daily discharge data of the Ussita River. The discharge hydrograph of this sub-catchment overlaps the Nera River one with some differences during the late August-September periods of 2017 and 2018. In these two periods, the recession phase of the Ussita River continued while that of the Nera River stopped after receiving some recharge that did not occur in the Ussita catchment.  Table 1).

Discussion
Updates on the response of the Nera River discharge to the 2016 seismic sequence are discussed, taking into consideration the meteorological conditions during pre-and post-seismic phases. During the two years after the 2016 seismic sequence (2017 and 2018), the peak of river discharge after the autumn-spring recharge period was very high compared with 2019 and 2020 ( Figure 5). The high discharge value recorded in January of 2017 was mainly linked to co-seismic effects (e.g., increased aquifer permeability and pore water pressure). Moreover, as reported by Petitta et al. [21], the impressive increase  Table 1).
Considering periods with no aquifer recharge (water surplus S = 0 in Figure 5) the pre-and post-seismic recession curves were collected. After the individuation of months with S = 0, a deep analysis was carried out to separate consecutive days characterized by a decrease in discharge values, avoiding days with discharge increases in response to the previous effective rainfalls. These curves help to investigate the relationship between aquifer hydrogeological properties and groundwater discharge to the river. The results of the analysis of the 2019 and 2020 depletion phases are presented in Figure 6, moving to the discussion section for the comparison with previous studies published in the study area in the two years following the earthquake. The EF recession model describes both 2019 and 2020 recession curves as a result of the best-fit analysis, indicating that the hydrogeological system feeding the river empties with a Darcian laminar flow. The two curves show similar a recession coefficient α E , −2.1 × 10 −3 day −1 and −2.3 × 10 −3 day −1 for 2019 and 2020, respectively. Figure 5. Hydrograph of the Nera River and Ussita River with monthly surplus calculated by using meteorological data of the Ussita station. The blue numbers adjacent to the dashed lines indicate the main earthquakes with moment magnitude (Mw) higher than 5.0 that occurred during the 2016-2017 seismic sequences in Central Italy (information about the earthquakes are reported in Table 1).

Discussion
Updates on the response of the Nera River discharge to the 2016 seismic sequence are discussed, taking into consideration the meteorological conditions during pre-and post-seismic phases. During the two years after the 2016 seismic sequence (2017 and 2018), the peak of river discharge after the autumn-spring recharge period was very high compared with 2019 and 2020 ( Figure 5). The high discharge value recorded in January of 2017 was mainly linked to co-seismic effects (e.g., increased aquifer permeability and pore water pressure). Moreover, as reported by Petitta et al. [21], the impressive increase

Discussion
Updates on the response of the Nera River discharge to the 2016 seismic sequence are discussed, taking into consideration the meteorological conditions during pre-and postseismic phases. During the two years after the 2016 seismic sequence (2017 and 2018), the peak of river discharge after the autumn-spring recharge period was very high compared with 2019 and 2020 ( Figure 5). The high discharge value recorded in January of 2017 was mainly linked to co-seismic effects (e.g., increased aquifer permeability and pore water pressure). Moreover, as reported by Petitta et al. [22], the impressive increase in spring and river discharge observed in the upper Nera River may be correlated with the subsidence induced by the toe of VBF-as already recorded by GNSS station (Figure 1a)-which might have created an additional "squeezing effect" in the core of the Sibillini Mountains aquifer. Barberio et al. [77] reported some changes in hydrochemical characteristics of water of the basal aquifer feeding the San Chiodo spring (Figure 1a), such as an increase in the content of three metals (Cr, Fe, and V) and one metalloid (As) four months before or roughly in conjunction with the 2016 seismic sequence, respectively. Rosen et al. [78] reported no change in most of the major ions and in δ 18 O H2O and δ 2 H H2O of waters of the Nerea spring (Figure 1a), the recharge area of which is characterized by the basal aquifer extending towards SE including the Vettore Mountain faulting systems deeply affected by co-seismic ruptures ( Figure 1). In detail, concentrations of SO 4 showed no change after the 30 October 2016, earthquake and trace element concentrations returned to pre-earthquake concentrations by the end of November, 2016. The authors concluded that the hydrochemical dynamics suggested within-aquifer changes instead of mixing with another aquifer, geothermal fluids, or aquifer breaching. Recently, Fronzi et al. [79] presented updated hydrogeochemical results for the Castelsantangelo sul Nera areaalso fed by basal aquifer and located slightly upstream to the Nerea spring-showing some sampling points with an increase in SO 4 content and some others with no increase after the mainshock. The different results reported by [78,79] highlight the complexity of processes that occurs in a hydrogeological system highly influenced by normal faults and co-seismic fracturing systems. In other words, sampling points collected along the normal faults interacting with Triassic evaporites aquiclude below the Mt. Vettore-Pian Grande Plain show increases in SO 4 , which is not evident in other zones-also fed by basal aquifer-where normal faults and co-seismic fracturing did not play a significant role.
Mastrorillo et al. [23], by a mechanism named "Aquifer Fault Rupture", hypothesized a shift eastward of the piezometric divide, severely penalizing the aquifers in the eastern side of the Sibillini Mountains. The authors left open whether or not the modified piezometric divide would recover to its original position. Our results on the Upper Nera River catchment (western part of the Apennine chain) show that co-seismic effects on GW lasted up to at least 2018 by analyzing river discharges at the Visso gauge. The peak discharge that occurred in 2018 is preceded by a recharge period characterized by a water surplus in line with or slightly above the historical mean. The maximum discharge was 1.3 times higher than that recorded in 2016 (pre-seismic conditions); the latter was preceded by about four years of very high water surplus values, which did not produce the high spring flows recorded in 2018 (Figure 4b). After about five years from the seismic sequence, the comparison of pre-and post-seismic recession curves can contribute to understanding what processes are taking place and whether or not permanent changes have taken place. Table 3 shows the recession processes and recession coefficients computed on the available recession periods. Figure 7 shows all the available pre-and post-seismic recession curves, sorted from the highest discharge to the lowest. Table 3. Recession models that best describe the Nera River discharges (Section I1 in Figure 2) during the available recession periods.
Curves with similar recession processes and recession coefficients have been merged using the strip method [79][80][81]. This approach allows obtaining a Master Recession Curve (MRC); the Root Mean Square Error (RMSE) criterion has been used to choose the best time-shifting to join the curves. A recent study carried out by Di Matteo et al. [23] highlighted that the SL equation described the recession curves of the two years following the seismic sequence (2017 and 2018). A permeability enhancement of the aquifer feeding the Nera River-due to cleaning of fractures and the co-seismic fracturing in the recharge  (Table S1).
Curves with similar recession processes and recession coefficients have been merged using the strip method [80][81][82]. This approach allows obtaining a Master Recession Curve (MRC); the Root Mean Square Error (RMSE) criterion has been used to choose the best time-shifting to join the curves. A recent study carried out by Di Matteo et al. [24] highlighted that the SL equation described the recession curves of the two years following the seismic sequence (2017 and 2018). A permeability enhancement of the aquifer feeding the Nera River-due to cleaning of fractures and the co-seismic fracturing in the recharge area-coupled with an increase in groundwater flow velocity can explain this process, which lasted up to 2018. Fronzi et al. [26] also recorded consistent variations in the hydraulic conductivity distribution throughout carbonate aquifers related to the VBF system activation. Figure 7 also includes the most recent recession curves (2019 and 2020), which fit very well with the pre-seismic ones (Equation (11)), represented by the EF equation (Darcian flow).
Considering that the recession process recovered to the same pre-earthquake conditions-at least for two consecutive years (2019 and 2020)-and α E values are similar to the pre-seismic ones (Table 1), it can be deduced that hydrogeological changes induced by the seismic sequence seem to be almost recovered. As Equation (10) shows, co-seismic effects could have changed other parameters than aquifer hydraulic conductivity. It is interesting to point out that if the length of the one-dimensional domain (L) had increased (i.e., the recharge area had enlarged due to aquifer breaching), we should have expected much lower αE than pre-earthquake conditions. The permeability recovery after high magnitude earthquakes has been documented in a wide range of hydrogeological settings; most documented recovery times are about two years as obtained for the Nera River [27,28,83]. As Manga et al. [84] discussed, the mechanism of recovery is not clear because of the complexity and inaccessibility of the subsurface; a combination of mechanical and geochemical processes can cause the return of the enhanced permeability to approximately the pre-seismic value. Aben et al. [85] reported a list of possible recovery mechanisms of earthquake-related fracture damage including time-dependent mechanical recovery of micro and macro-cracks, sealing by mineral precipitation (e.g., calcite), resulting in fracture sealing, and fracture closure by pressure solution creep, the latter two act on longer time scales than mechanical ones. Based on the discharge at Visso gauge (I1), regardless of the mechanism that led to the recovery of the hydrogeological properties of the aquifer feeding the Nera River, it is important to emphasize again that, the emptying process seems to be restored to the same as it was before the 2016 earthquake. Since the reservoir in 2019 and 2020 empties at rates lower than 2017 and 2018, the hydrogeological system-as it was in the past-is less vulnerable to prolonged droughts such as those that are increasingly affecting the Apennine area of Central Italy.

Conclusions
The 2016 seismic sequence deeply affected the hydrogeological system of the upper part of the Nera River catchment-the waters of which are of strategic importance to the territory's economy. Through analysis of the geological-structural and hydrogeological data about five years after the seismic sequence, the following conclusions can be drawn: − The peak river flow recorded in 2017 and 2018 did not occur in 2019 and 2020-the latter was preceded by a drought period well detected by the SPI, which affected the recharge of groundwater, and consequently the river discharge; − The rapid depletion observed during the 2017 and 2018 recession periods (straight line equation and turbulent flow) were not detected in 2019 and 2020, where both recession process and recession coefficients seem to be restored to the same as prior to the 2016 earthquake (exponential equation and Darcian flow); − Co-seismic effects on the hydrogeological system (e.g., increased aquifer permeability and pore water pressure) appear to have recovered after two years from the 2016 seismic sequence, as documented in other systems around the world affected by strong earthquakes.
Based on the discharge data representing the whole hydrogeological system feeding the Nera River, the findings reported here contribute to understanding processes and dynamics of fractured carbonate aquifers located in geologically and climatically complex regions, which are useful for water management. The results need to be further refined by continuing the monitoring and also with consideration of the response of hydrogeological systems located east of the Apennine chain.