Combined Discharge and Thermo-Salinity Measurements for the Characterization of a Karst Spring System in Southern Italy

The hydrological monitoring of springs is an auxiliary and indispensable tool that goes alongside investigations in wells to reconstruct a conceptual phenomenological model of an aquifer–groundwater system and its interactions with surface waters. There are manifold ways to carry out this monitoring, but the choice of which way is significant for a correct qualitative and quantitative knowledge of spring systems. The present work focuses on the characterization of the thermo-saline and flow regimes of the Tara spring system along the northern coast of Taranto (southern Italy), where a karst groundwater basin is the major source of the Tara River and the surrounding coastal wetland. A series of measurements was carried out on the spring system to support a technical feasibility study on the possible use of the brackish water of this river to feed a future desalination plant. To estimate the flow rate, a comparison was made between different flow measurement methods in a derivation channel. Through an analysis of the available dataset, the response of the aquifer to the autumn–winter recharge, for which updated hydrologic measurements were not available, is highlighted.


Introduction
As in many other Mediterranean regions, the economy of southern Italy is mainly based on irrigated agriculture. As a result, water resources have seen a serious qualitative and quantitative depletion caused by overexploitation. Due to the lack of significant surface watercourses, there has been a greater dependence on groundwater, and a number of irrigation wells have been drilled to provide the necessary water. In some groundwater basins of southern Italy, irrigation channels draining groundwater flows are an alternative source, providing high-quality water to croplands and regular discharge regimes.
In more recent years, policies related to freshwater conservation are likely to become obsolete or no longer responsive to new and changing environmental and socioeconomic conditions [1]. It has been observed that climate change, combined with anthropogenic pressure and infrastructure, is consistently associated with changes in several components of the water systems [2], particularly concerning the Mediterranean region [3].
Springs, in particular, are very vulnerable to long-term climate alterations, land surface activities, intensive agriculture, growing water demands, and landscape alterations [4,5]. Therefore, there is a need to update studies on large groundwater spring sheds and investigate the impact of climate change on the potential groundwater recharge [6][7][8]. The monitoring of springs involves the study of its physical, chemical, and biological phenomena, using appropriate equipment with the aid of analytical techniques. Moreover, spring monitoring is an indispensable instrument (after the investigation of wells) in reconstructing reliable conceptual models of groundwater-aquifer systems [9]. A qualitative and quantitative description of groundwater spring processes can only be achieved through a thorough understanding of relationships between the geological, hydrogeological, and chemical-physical elements that characterize the aquifer system, as well as their evolution over time. In this regard, special attention has recently been addressed to the principal karst aquifer systems in central and southern Italy due to observed climate-induced alterations in the drinking water supply (supplying about 12 million people) [10].
Among the existing investigation techniques, water tracing methods have been carried out over the years to define catchment boundaries, estimate groundwater flow velocities, determine areas of recharge, and identify sources of pollution of spring water [9]. Luhmann et al. [11] presented an overview of the different methods used for the study of spring behavior within a karst aquifer as a means of characterizing a groundwater basin from spring characteristics. They outlined the hydraulic [12], chemical [13], and isotopic [14] responses from spring monitoring to provide information about the aquifer. Additionally, a combination of temperatures and heads was jointly suggested in order to provide unique information on springs. Indeed, Luhmann et al. [11] asserted that karst spring temperature is the most crucial source of information for determining aquifer characteristics. Alexander et al. [6] examined a further methodology for a better understanding of the overall flow system in southeastern Minnesota, using a combination of dye tracing data, the water well driller's records, and downhole gamma logging. Furthermore, electrical tracing methods based on geophysical techniques are also highly promising in investigating surface water-groundwater interactions, particularly on hillslopes [15]. Both temperature and the electrical conductivity of spring water have been observed to be closely related to spring discharge [16,17]. For example, Birk et al. [18] showed that localized spring discharge in a karst aquifer could be characterized by simultaneously analyzing the electrical conductivity and temperature of spring water.
Broadly speaking, chemical-physical parameters, which are recognized as typical "natural tracers", are closely related to the lithological characteristics of a rock mass, to the type of water circulation in an aquifer, and to the type of recharge, which can be autogenic or allogeneic [19][20][21]. In this type of monitoring, the absolute values are not as important as the variations that these tracers may undergo over time as a result of the increase in flow rate. In particular, the temperature of spring water is an easy parameter to measure continuously. Furthermore, the interactions between rock and water are governed by relatively simple laws such as heat exchange, which is linked to water-rock contact and residence times. Heat exchanges between rock (characterized by low thermal conductivity) and water (high thermal conductivity) affect the temperature at the source. The specific electrical conductivity is relatively more difficult to measure continuously and is linked to the total mineralization of the water. It is a physical parameter that depends on the content of the main dissolved ions (essentially bicarbonates, calcium, magnesium, sodium, potassium, chloride, and sulphate) and therefore indirectly reflects the chemical load of the groundwater. It is therefore the concentrations of these ions and their variations over time that determine the trend in the electrical conductivity at the source and therefore the degree of mixing between waters of different origins.
The present research sets the stage for a technical and environmental feasibility study about the construction of a desalination plant by the local water utility using waters emerging from springs. A three-year monitoring program was implemented to determine the quantity (i.e., the flow rate) and quality (i.e., the temperature and conductivity) of the spring system, called Tara, for which updated hydrological measurements were not available. The results were then used to estimate the feasibility of the project, as well as to provide important information on the status of the aquifer. While the specific results of this study are of local interest, the approach of using a combination of (i) established methods to estimate the channel flow rate and the residual capacity of the karst spring system and (ii) physical Sustainability 2020, 12, 3311 3 of 21 parameters to determine the qualitative characterization of the spring system may also be considered relevant. This applies to the entire scientific community dealing with the evaluation, development, and management of groundwater resources, particularly in less developed regions, where hydrological monitoring is often discontinuous or unavailable.

Study Area
From a geological point of view, the study area is located in the SW sector of the Murge relief, which corresponds to the southern Apennines foreland (Figure 1). From the bottom to the top, the following geological units can be distinguished [22]: Altamura limestone (Cretaceous), Gravina calcarenite (Upper Pliocene-Lower Pleistocene), Subappennine clays (Lower Pleistocene), terraced marine deposits (Middle-Upper Pleistocene), and alluvial and coastal deposits (Holocene). Several normal faults dislocate the Cretaceous bedrock, creating a horst and graben setting oriented NW-SE. As a consequence, a Plio-Quaternary deposit outcrop on the lower structural sectors of the area conceals faulted Cretaceous limestone, which, conversely, is found only on the highest structural sectors (Figure 1a). This geological setting also includes two hydrogeological complexes indicative of specific aquifer structures: a deep limestone aquifer (the main aquifer) characterized by secondary permeability caused by fracturing and karstic processes and a shallow porous aquifer corresponding to a sand-conglomerate-calcarenite complex of terraced marine deposits and alluvial-coastal successions [23].
Sustainability 2020, 12, 3311 3 of 21 development, and management of groundwater resources, particularly in less developed regions, where hydrological monitoring is often discontinuous or unavailable.

Study Area
From a geological point of view, the study area is located in the SW sector of the Murge relief, which corresponds to the southern Apennines foreland (Figure 1). From the bottom to the top, the following geological units can be distinguished [22]: Altamura limestone (Cretaceous), Gravina calcarenite (Upper Pliocene-Lower Pleistocene), Subappennine clays (Lower Pleistocene), terraced marine deposits (Middle-Upper Pleistocene), and alluvial and coastal deposits (Holocene). Several normal faults dislocate the Cretaceous bedrock, creating a horst and graben setting oriented NW-SE. As a consequence, a Plio-Quaternary deposit outcrop on the lower structural sectors of the area conceals faulted Cretaceous limestone, which, conversely, is found only on the highest structural sectors (Figure 1a). This geological setting also includes two hydrogeological complexes indicative of specific aquifer structures: a deep limestone aquifer (the main aquifer) characterized by secondary permeability caused by fracturing and karstic processes and a shallow porous aquifer corresponding to a sand-conglomerate-calcarenite complex of terraced marine deposits and alluvial-coastal successions [23]. The Tara spring system represents the most important karst spring in the Apulia region fed by groundwater circulating within the deep limestone aquifer, which originates from the Murge relief.
The presence of the impermeable covering deposits (Subappennine Clays), going toward the sea, makes the aquifer confined and locally even artesian, as confirmed by several flowing artesian wells in this area. Moreover, the deepening of the clay deposits hinders the seaward water flow and thus forces it to flow upwards along higher permeability paths controlled by structural discontinuities in the bedrock or by permeability changes in the underlying deposits. Consequently, a large water-rising zone results, located in the NW sector of Taranto (Figure 1b). It consists of a main spring area about 2 km from the sea that gives rise to the Tara River and to several secondary springs (about 20) spread over a larger area. All springs form a single spring system with moderately saline water (~3 g l −1 ). Studies dealing with the chemical characteristics of brackish groundwater sampled in different coastal springs of the Apulian karstic aquifers have revealed that the salinity of the Tara spring system is linked to aged saltwater, which is very different from modern seawater [24,25]. The Tara spring system generates regular outflows with peaks of about 4 m 3 s −1 , according to measurements made in the first half of the last century. The remarkable outflows of the Tara spring system are fed by a hydrogeological basin that is much more extended than its hydrographic one The Tara spring system represents the most important karst spring in the Apulia region fed by groundwater circulating within the deep limestone aquifer, which originates from the Murge relief.
The presence of the impermeable covering deposits (Subappennine Clays), going toward the sea, makes the aquifer confined and locally even artesian, as confirmed by several flowing artesian wells in this area. Moreover, the deepening of the clay deposits hinders the seaward water flow and thus forces it to flow upwards along higher permeability paths controlled by structural discontinuities in the bedrock or by permeability changes in the underlying deposits. Consequently, a large water-rising zone results, located in the NW sector of Taranto (Figure 1b). It consists of a main spring area about 2 km from the sea that gives rise to the Tara River and to several secondary springs (about 20) spread over a larger area. All springs form a single spring system with moderately saline water (~3 g l −1 ). Studies dealing with the chemical characteristics of brackish groundwater sampled in different coastal springs of the Apulian karstic aquifers have revealed that the salinity of the Tara spring system is linked to aged saltwater, which is very different from modern seawater [24,25]. The Tara spring system generates regular outflows with peaks of about 4 m 3 s −1 , according to measurements made in the first half of the last century. The remarkable outflows of the Tara spring system are fed by a hydrogeological basin that is much more extended than its hydrographic one ( Figure 1a). Moreover, considering the limestone outcroppings, it is possible to identify a feeding area close to the hydrographic catchment of the Tara River and a farther one connected to the inland limestone outcroppings [23].
The Tara River, which flows into the Gulf of Taranto after a short stretch of about 3.5 km (Figure 1), is almost exclusively fed by deep groundwater input, as evidenced by the high electrical conductivity values of its water and an almost constant temperature of 18-20 • C throughout the year. In the early 1950s, a capture plant for one of the secondary springs was built about 1.2 km from the source of the Tara River by constructing a pumping station consisting of an intake chamber carved out of limestone and equipped with a pumping system. In addition, an artificial channel 3 m wide, 3.5 m deep, and 1000 m long was built as a derivation channel connecting the Tara River to the pumping station. The side walls of this channel, which have a concrete lining, prevent any lateral inflow and outflow, thus allowing the channel to drain any groundwater outflow from the channel bottom. Indeed, the function of the channel is to intercept additional sources and collect them toward the pumping station. The latter consists of eight different electric pumps with a total power of 750 kW that are able to operate at a total flow rate in the range of 0.40-3.70 m 3 s −1 . A movable dam was built on the natural stream downstream of the confluence with the derivation channel to enable a more efficient and flexible operation of the pumping station. In this way, water can flow in either direction through the derivation channel, i.e., from the intake chamber to the confluence with the river and vice versa, depending on the hydraulic gradient between the water heads in the river and the intake chamber.
From its construction until the early 1990s, the pumping station was operated at around full capacity, mainly supplying water for industrial needs linked to steel production and subordinately for agricultural use. However, at the present time, only one or two electric pumps are simultaneously active, so that water is withdrawn almost completely from the spring, which is intercepted just below the pumping station. Nowadays, the extracted flow rate ranges between 0.45 m 3 s −1 and 0.65 m 3 s −1 due to a drastic reduction in blast furnace activity in the Taranto industrial zone. Under such an exploitation regime, the flow direction in the derivation channel occurs from the pumping station to the natural stream. Consequently, a renewed interest in the hydrogeochemical characterization of the Tara spring system is justified both for the purposes of wetland conservation as well as for its potential exploitation as feed water for a water desalination plant.

Methods
The survey method adopted for the present study aimed to (i) detect the main springs intercepted along the derivation channel by means of thermoconductivity measurements; (ii) identify the optimal measurement method as a tradeoff between simplicity of execution and reliable estimates, taking into account the characteristics of the channel (narrow and deep) and the observed flow regimes; (iii) study the variability of the flow regime during a hydrological year through an appropriate measurement program; and (iv) evaluate the medium-long-term variations in streamflow regimes, also as a result of the decreased spring withdrawal.

Measurement Techniques
Natural traces provide valuable information about how groundwater moves in the subsurface and the characteristics of karst springsheds. Electrical conductivity (EC) can be used as an index of total dissolved solids and, in some cases, as a predictor of concentrations of individual ions. EC can also be used to interpret the changing sources of stream runoff at different time scales and to provide information about the contrasting hydrological behavior of specific catchments [26]. Because groundwater commonly differs chemically from stream water, groundwater discharge zones often coincide with relatively sudden changes in water chemistry along a stream, which can be detected Sustainability 2020, 12, 3311 5 of 21 by measuring along-stream variations in EC. Inferences regarding groundwater discharge can be made more confidently by combining EC measurements with other observations, such as hydraulic gradients across the streambed, water temperature, and streamflow measurements. Furthermore, EC can be used to compute the relative contributions of two tributaries to flow below the confluence, or to separate quantitatively the contributions to streamflow from two distinct sources. As long as the limitations of EC are borne in mind, its measurement can provide useful and rapid insight into the chemical and hydrological characteristics of water systems. EC is relatively easy to measure either manually, using a handheld conductivity probe, or almost continuously using a probe connected to a data logger. EC measurements are strongly temperature-dependent, so they must be standardized to a reference temperature, typically 25 • C or 20 • C, to make the measurements comparable to each other.
In a recent review of methods used for streamflow monitoring, Dobriyal et al. [27] wrote that streamflow monitoring methods are specific to stream types, which can be classified on the basis of eight major variables. These are their width, depth, velocity, discharge, slope, roughness of bed and bank materials, sediment load, and sediment size. The different methods available to quantify and monitor surface water flow are grouped into four categories [28][29][30]: (a) direct measurement methods, (b) velocity-area methods, (c) formed constriction or constricted flow methods, and (d) noncontact measurement methods. A brief account of the suitability of each method to different terrains, with their advantages and disadvantages, is summarized in Dobriyal et al. [27]. Moreover, the selection of measurement method should be based on the volume of water to be measured, the degree of accuracy desired, whether the installation is permanent or temporary, and the required costs [31,32].
Considering that the investigated channel has regular, artificial walls, a width of 3 m, a natural stony-sandy bed, and a water depth of around 2 m, certain measurement methods were chosen and tested on multiple sections, which had been previously selected according to the thermo-salinity anomalies detected along the stream current. Based on the channel characteristics, two measurement methods belonging to the velocity-area category were adopted, namely the current meters method (CM method) and the float method (F method). Though they are very different in terms of complexity of operation and the accuracy of the measured discharge, the two methods were tested and compared, and we highlight the pros and cons in the given field conditions.

CM Method
The CM method is considered to be the most accurate of the velocity-area methods when appropriate procedures are adopted both for the measurement devices and field acquisition and processing. This method consists of determining a flow velocity in a cross-section of a stream or channel by means of a current meter and computing the discharge using known geometric relationships [33]. We adopted the midsection method for computing discharge, assuming that the mean velocity in each vertical represented the mean velocity in a partial section (segment). The mean velocity in each vertical was achieved from velocity observations at several points in that vertical and by using a known relation between those velocities and the mean in the vertical. As suggested by most common methods, velocity measurements were taken at 60% of the water depth, as we assumed this velocity to be representative of the average velocity along the given vertical. The sum of the discharges for all the partial sections was assumed as the total discharge of the stream section. Moreover, the vertical-velocity curve method was adopted in some cases where detailed current measurements were required to investigate the effect of submerged weeds on the velocity profiles and on the two-point mean velocity estimated values. In these cases, a series of velocity measurements was carried out at each of the verticals, i.e., at 0.1-m depth increases between 0.1 and 0.9 m of depth. The results obtained using the vertical-velocity curve method were compared to the two-point method for a mean velocity measurement.
A drawback of indirect methods is the presence of flexible (elastic) vegetation in the channel, which can lead to errors in velocity measures, especially when using reels. In fact, the submerged vegetation that flexes due to the dragging action exerted by water flow offers additional resistance to Sustainability 2020, 12, 3311 6 of 21 motion and therefore affects flow velocity. Moreover, vegetation roughness is not constant but depends on the flow condition (depth and velocity) as well as the vegetation condition (type and density). According to the theory of Kouwen et al. [34], in the presence of flexible vegetation, the velocity profile exhibits a logarithmic trend, described as follows: where V m is the mean flow velocity; u* is the shear velocity, defined as (gRJ) 0.5 , with g being gravitational acceleration, R being the hydraulic radius, and J being the energy gradient (i.e., the hydraulic slope, i); A is the area of the channel cross-section and A v the area of its vegetated part; C 0 is a parameter based on the vegetation density; and C 1 is dependent on the vegetation stiffness. The parameters C 0 and C 1 are tabulated (Table 1) based on the ratio between u* and the critical shear velocity (u* c ), which is introduced to describe the limiting friction velocity between the erect and prone states of vegetation. Indeed, vegetation can be regarded as prone if u* is higher than u* c . Table 1. Values of C 0 and C 1 in Equation (1).

Configuration
Criterion To estimate u* c and then the roughness deflection for flow over submerged flexible vegetation, Kouwen and Unny [35] introduced a stiffness parameter, MEI, which takes into account the density M and the rigidity EI (elasticity, E, and the moment of inertia, I) of the vegetation and is defined as follows: where γ is the specific weight of water, d the average depth of the water, h d the deflected roughness height of the vegetation, h v the height of vegetation, and h v /h d the inflection degree of vegetation. Therefore, u* c is defined by two different expressions, depending on the deformation behaviors of the vegetation: u* c = 0.028 + 6.33(MEI) 2 , where the first equation is valid in the case of flexible vegetation, and the second one is valid for vegetation that breaks when flattened. Based on the value assumed by the ratio u*/u* c , the values of parameters C 0 and C 1 can be chosen (Table 1) and the velocity calculated. The flow rate is then calculated using the Chezy formula: where χ is the roughness coefficient assessed considering Manning's formula, in which n nominally quantifies channel roughness or resistance to flow. In the case of vegetated river beds, n can be calculated as follows: Sustainability 2020, 12, 3311 7 of 21

F Method
Floats have somewhat limited use for stream gauging, but they prove useful when the velocity is too low to obtain reliable measurements with a current meter [33].
Similarly to the CM method, the F method is based on the velocity-area principle; therefore, the section geometry must be defined by measuring the depth and width of n subsections. A single float near the middle of the channel is used to determine surface velocity (V s ) as the averaged value of several measurements [36]. The technical literature reports that a surface float travels with a velocity about 1.2 times the mean velocity of the water column beneath it [37]. The velocity of any float (e.g., a wooden stick, a tube with a weight at its lower end), whether on the surface or submerged, is likely affected by wind, though with different impacts on measurement errors. The V s data collected from the measurements were processed by applying three different corrective methods: where d is the average depth of water in the upstream and downstream sections, and l is the sinking of the float below the free surface; • The 1963 Roche study [39], which indicated that the appropriate F for this case study was equal to 0.85 for all of the sections; and • An empirical expression taken from the Mysore Research Institute in India [40]: Finally, the transit flow rate (Q) between the upstream and downstream sections (selected for the float path) is where A m is the average area between the upstream and the downstream sections. Discharge measurements using the F method under favorable conditions may be accurate to within ±10%, but if a poor reach is selected and not enough float runs are made, the results can be as much as 25% in error [33].
In this study, the F method was selected because of its easy, efficient, and inexpensive use, while the more difficult, time-consuming CM method was adopted for its much higher accuracy. In this work, the accuracy of the two methods was evaluated by comparing the F method measurements to the CM method measurements, which were taken in some specific cross-sections of the investigated channel. We sought a fair compromise between measurement complexity and accuracy.

Measurement Equipment
The instrument used for the EC and temperature (T) measurements was a Delta OHM HD2106.2 equipped with a datalogger and internal memory, which provided EC values at a standard temperature of 20 • C or 25 • C. The thermoconductivity probe had a T measurement range of −50 • C/+90 • C (accuracy +/−0.2 • C, resolution 0.1 • C). The EC measurements covered various orders of magnitude of salinity ranging from 0.0 µS cm −1 up to 2000 mS cm −1 (accuracy +/−1% of full scale, resolution 1 µS cm −1 ) if a proper cell constant value was set. For the measurements of the Tara Springs, a standard temperature of 20 • C was chosen, and a cell constant K = 1 cm −1 was adopted, which allowed measurements to be made from low to relatively high EC. For the application of the CM method, which was based on punctual speed measurements at different depths in the channel, a microreel and a hydrometric reel were used. The measurements provided by microreels can be disturbed by many materials, such as debris or filamentous algae, which can extend from the channel bottom to the surface, thus altering the rotation velocity of the rotor. For the determination of the flow rate of the Tara derivation channel, two reels with a propeller axis parallel to the current were used, specifically a laboratory microreel (a Nixon Streamflo Velocity Meter 400 model) and a hydrometric reel (an SIAP ME4001 model). The first, which was designed to measure low stream speeds in open channels, is intended primarily for use in clean waterways. The measuring head, which has a cage about 1.2 cm in diameter, allows for measurements in very small spaces. On the contrary, the ME4001 SIAP hydrometric reel can be used for small and large speeds in large-and small-scale waterways, in canals, and in conduits with clear, turbid, and salty waters.
For the application of the F method, different floating objects, including partially submerged spherical and vertical stick objects, were tested to identify the best method (to compromise between accuracy and feasibility), considering the multiple measurement sections and flow monitoring campaigns of the planned investigation. Specifically, two types of floats were used: a spherical object approximately 10 cm in diameter (e.g., the diameter of an orange), so as to be almost completely submerged, and a velocity rod consisting of a plastic tube closed at its ends and weighed down with quartz sand to guarantee a sinking of 45 cm.

Conductivity and Temperature Measurements
Considering the investigated hydrogeological system, the monitoring of T and EC was mainly oriented toward identifying and characterizing the different spring sources present in the intake chamber, along the derivation channel, and in the Tara River. For this purpose, several measurement transects were selected along the channel, both upstream and downstream of the confluence with the Tara River ( Figure 2). For the application of the F method, different floating objects, including partially submerged spherical and vertical stick objects, were tested to identify the best method (to compromise between accuracy and feasibility), considering the multiple measurement sections and flow monitoring campaigns of the planned investigation. Specifically, two types of floats were used: a spherical object approximately 10 cm in diameter (e.g., the diameter of an orange), so as to be almost completely submerged, and a velocity rod consisting of a plastic tube closed at its ends and weighed down with quartz sand to guarantee a sinking of 45 cm.

Conductivity and Temperature Measurements
Considering the investigated hydrogeological system, the monitoring of T and EC was mainly oriented toward identifying and characterizing the different spring sources present in the intake chamber, along the derivation channel, and in the Tara River. For this purpose, several measurement transects were selected along the channel, both upstream and downstream of the confluence with the Tara River ( Figure 2).  Moving from the spring intake chamber to the channel, an abrupt rise in salinity could be observed. The EC values remained roughly constant in the channel, while a reduction in salinity was registered after the confluence with the Tara River and in the section Q7, which was the only section upstream of the confluence (Figure 3a).
Moreover, it was evident that EC variations over time were negligible for the same measurement section. A more evident variation could be observed for the 2019 measurements, which was an effect of abundant autumn rainfall that recharged the area of the aquifer feeding the Tara spring system.
The T variations along the measurement sections roughly reflected those of EC, i.e., a lower T could be noted for the spring intake chamber, as could a reduction in temperature downstream of the confluence with the Tara River (Figure 3b). Unlike EC, T was more variable over time, as it was more easily conditioned by the air temperature. Indeed, we observed a trend linked to the season during which the measurements were taken. Nevertheless, the T values of the Tara spring system remained high compared to nearby rivers with a shallow groundwater feeding system or with prevalent rainfall input.
The T and EC variations characterizing the measurement sections along the Tara spring system Moving from the spring intake chamber to the channel, an abrupt rise in salinity could be observed. The EC values remained roughly constant in the channel, while a reduction in salinity was registered after the confluence with the Tara River and in the section Q7, which was the only section upstream of the confluence (Figure 3a).
Moreover, it was evident that EC variations over time were negligible for the same measurement section. A more evident variation could be observed for the 2019 measurements, which was an effect of abundant autumn rainfall that recharged the area of the aquifer feeding the Tara spring system. measurements. At first, the geometry of each cross-section was surveyed, i.e., its width, bottom depth, and water level, moving from the right to left bank. At the same time, velocity measurements were taken along different vertical profiles. As happened with the qualitative measurements, the position of the flow rate measurement sections changed during the monitoring period. To distinguish the measurement sections of the 2018-2019 campaign from those of the 2016 campaign, a different label was used by adding the suffix "Bis" or "Ter" (Figure 4). During the first measurement campaign (29 January 2016), the velocity was measured in five sections (from 1 to 5 in Figure 4) using a microreel Streamflow Velocity Meter 400 at 60% (Y0.6) of the water level for each vertical investigated (Table A1, Appendix A) (as an approximation of the average velocity). Only for Section 3 was the velocity measured along the vertical center line (i.e., xi = 1.5 m) with a 0.10-m vertical step, starting from 0.05 m below the free water surface up to 1.49 m, with deeper measurements being impossible due to dense vegetation (Table A2, Appendix A). The presence of filamentous algae is one of the possible causes of the measurement error, together with the practical difficulty of operating in stationary conditions while using a hydrometric rod over 3 m in length with the propeller of a microreel fixed at one end.
The velocity data were processed using the CM method in order to obtain estimates of the flow rates in the sections under study. Calculations were made for the area of each segment in which the section was split (Ai) as well as for the average velocity (Vm) as a weighted average of the mean The T variations along the measurement sections roughly reflected those of EC, i.e., a lower T could be noted for the spring intake chamber, as could a reduction in temperature downstream of the confluence with the Tara River (Figure 3b). Unlike EC, T was more variable over time, as it was more easily conditioned by the air temperature. Indeed, we observed a trend linked to the season during which the measurements were taken. Nevertheless, the T values of the Tara spring system remained high compared to nearby rivers with a shallow groundwater feeding system or with prevalent rainfall input.
The T and EC variations characterizing the measurement sections along the Tara spring system could be linked to the presence of distinct springs fed by circuits with different depths and travel times (which means different interaction degrees with the rock and mixing with waters with different degrees of salinity). This is in agreement with previous studies that have highlighted the presence, within the springs area, of sectors characterized by different EC values [41].
On the basis of the above remarks, it was possible to identify with some confidence the presence of additional springs between Sections 1 and 2, Sections 2 and 3, and Sections 5 and 6.

Discharge Measurements Using the CM Method
On the basis of the initial indications regarding the location of additional springs along the channel, which were deduced from the T and EC monitoring campaigns, the first two flow rate measurement campaigns took place on 29 January and 22 March 2016. Considering the regular lineament of the channel (a single turn about 700 m from the intake chamber), the selection of the cross-sections took into account having sufficient channel portions with a smooth bottom and banks to guarantee minimal flow turbulence and a representative hydraulic gradient. Depending on the approximately 3-m width of the stream, between 5 and 6 vertical sections were adopted for the flow measurements. At first, the geometry of each cross-section was surveyed, i.e., its width, bottom depth, and water level, moving from the right to left bank. At the same time, velocity measurements were taken along different vertical profiles. As happened with the qualitative measurements, the position of the flow rate measurement sections changed during the monitoring period. To distinguish the measurement sections of the 2018-2019 campaign from those of the 2016 campaign, a different label was used by adding the suffix "Bis" or "Ter" (Figure 4).
During the first measurement campaign (29 January 2016), the velocity was measured in five sections (from 1 to 5 in Figure 4) using a microreel Streamflow Velocity Meter 400 at 60% (Y 0.6 ) of the water level for each vertical investigated (Table A1, Appendix A) (as an approximation of the average velocity). Only for Section 3 was the velocity measured along the vertical center line (i.e., x i = 1.5 m) with a 0.10-m vertical step, starting from 0.05 m below the free water surface up to 1.49 m, with deeper measurements being impossible due to dense vegetation (Table A2, Appendix A). The presence of filamentous algae is one of the possible causes of the measurement error, together with the practical difficulty of operating in stationary conditions while using a hydrometric rod over 3 m in length with the propeller of a microreel fixed at one end.
The velocity data were processed using the CM method in order to obtain estimates of the flow rates in the sections under study. Calculations were made for the area of each segment in which the section was split (Ai) as well as for the average velocity (V m ) as a weighted average of the mean velocities (at Y 0.6 ). The flow rate (Q) of the section was then calculated by multiplying the average velocity by the wet surface of the section. All data are summarized in Table 2. Unfortunately, for the abovementioned reasons, the flow rate values calculated using this method cannot be considered to be representative of the actual flow rate of the channel. However, it was possible to correct the calculated flow rate using the data collected from the vertical centerline of Cross-section 3, as follows: First, the vertical velocity profile, shown in Figure 5, was fitted by a fifth-degree polynomial function. Then the average velocity was obtained by considering an integral mean value theorem. Using this procedure, the average velocity of the vertical centerline of Section 3 appeared to be double (and therefore more realistic) that measured at the depth Y 0.6 . The absolute error was calculated and, assuming this constant error for all verticals and adding it to the velocity values assessed at Y0.6, new average velocities were calculated. On the basis of the data derived from this model, the flow rate in Section 3 was recalculated, reapplying the segments method and obtaining more realistic flow rate values ( Table 3).
The velocity measurements for the second campaign (22 March 2016) were performed using a traditional hydrometric reel for river measurements, i.e., the ME 4001 model. The measurements were limited to Sections 4, 5, and 6 due to the presence of dense vegetation that made the channel inaccessible. Moreover, in Sections 4 and 5, it was not possible to derive velocity profiles along the entire depth due to the presence of filamentous algae within the channel (Table 3A, Appendix A). Therefore, we first calculated the fraction of the flow rate relative to the upper part of the section, i.e., from the free surface up to 25% of the depth (Y0.25) for all sections. Subsequently, the lower fraction of the flow rate was calculated for the six sections. Assuming this value to be constant in all sections, it was possible to calculate the total flow in Sections 4, 5, and 6 as the sum of the two rates (Table 4). This choice was justified by the fairly regular geometry of the channel and the almost uniform distribution of the submerged vegetation.  The absolute error was calculated and, assuming this constant error for all verticals and adding it to the velocity values assessed at Y 0.6 , new average velocities were calculated. On the basis of the data derived from this model, the flow rate in Section 3 was recalculated, reapplying the segments method and obtaining more realistic flow rate values (Table 3). Table 3. Corrected values of velocities and flow rate for Cross-section 3, which were obtained during the first monitoring campaign (29 January 2016).

Section
Wet Area (m 2 )

Xi (m)
V at Y 0.6 (cm s −1 ) The velocity measurements for the second campaign (22 March 2016) were performed using a traditional hydrometric reel for river measurements, i.e., the ME 4001 model. The measurements were limited to Sections 4, 5, and 6 due to the presence of dense vegetation that made the channel inaccessible. Moreover, in Sections 4 and 5, it was not possible to derive velocity profiles along the entire depth due to the presence of filamentous algae within the channel (Table A3, Appendix A). Therefore, we first calculated the fraction of the flow rate relative to the upper part of the section, i.e., from the free surface up to 25% of the depth (Y 0.25 ) for all sections. Subsequently, the lower fraction of the flow rate was calculated for the six sections. Assuming this value to be constant in all sections, it was possible to calculate the total flow in Sections 4, 5, and 6 as the sum of the two rates (Table 4). This choice was justified by the fairly regular geometry of the channel and the almost uniform distribution of the submerged vegetation.  To confirm the results obtained from the second monitoring survey, the average velocity and the flow rate passing through the downstream section of the channel, i.e., Section 6, were calculated by adopting Kouwen et al. experimental theory [34], i.e., Equation (1) and Equation (5). The flow rate turned out to be equal to 1.439 m 3 s −1 , a value similar to that obtained using the segments method (Table 4).

Discharge Measurements Using the F Method
The 2016 flow rate measurement campaign was also done by using the simpler F method. The flow rate measurement points were the same as those used for the application of the CM method ( Figure 5), excluding Section 6, which was not measured during the second monitoring campaign. At each point, the F method mandated the identification of two sections (upstream and downstream section) at a distance approximately three times the width of the channel. The geometry of all sections was reconstructed in order to calculate the average area of the liquid section. The travel times of the two types of floats (from the upstream to the downstream sections) were measured to estimate the surface speed (V s ). Three measurements were taken for each measurement point, and the average value was assumed to be representative of the streamflow. Subsequently, the collected V s data were processed to calculate the average sectional velocity, applying three different correction coefficients, as described in Section 3.1.2. Table 5 shows the values of V s for the two types of floats, the average velocities in each section obtained by applying the three corrective coefficients (V m quadr , V m Roche , V mMysore ), and finally, the transit flow rates calculated from the three speeds (Q quadr , Q Roche , Q Mysore ). Among these three values, Q quadr was chosen as the reference because it gave intermediate values compared to the other methods.
Comparing the average flow rate obtained with the two types of floats, the values were not very different from each other (±6%). This tiny difference could have been related to the presence of vegetation, which, when closer to the surface, could have hindered the free movement of the partly submerged rod. Moreover, increasing values were observed along the channel (from Section 1 to Sections 4 and 5) until values above 1.50 m 3 s −1 were reached. Finally, compared to the flow rates in Sections 4 and 5 from the first and second campaigns, an increase in the flow rate of approximately 0.135 m 3 s −1 was observed. This was probably due to the abundant rainfall in the catchment area in the period between February and March 2016. Table 5. Values, for each of the five sections analyzed, of the surface velocity measured using the two types of floats (V s ), the average sectional velocity obtained using the three methods (V m ), and the relative transit flow rates (Q).

Section
Campaign

Comparison between Discharge Measurements
We compared the flow rates obtained through the different measurement methods. Starting from the more reliable results, which were obtained during the second measurement campaign, the reliability of the F method with respect to the more rigorous CM method was assessed by comparing the same measurement sections (Table 6). From this comparison, it could be deduced that, regardless of the type of float used, the F method tended to overestimate the water flow rates by about 11% compared to the CM method. However, given the simplicity and rapidity of the former, a relative error of 11% may be widely acceptable compared to a more rigorous but much more time-consuming method, especially if multiple measurements are required with systematic repetitions in different periods of the hydrological year [38]. For this reason, the subsequent flow measurement campaigns that took place during 2018 were carried out using only the F method.

Discussion
The need for timely and cost-effective hydrological campaigns could be crucial for new water development projects or to evaluate the current degree of exploitation under altered hydrogeological regimes due to climate change effects and water overexploitation trends. In this study, EC and T measurements along the derivation channel were preliminarily carried out with the aim of identifying possible additional sources of springs existing at the bottom of the channel. Afterwards, flow rate measurements were performed in the same transects to characterize the hydrological regime. The flow rates measured using the F method are shown in Table 7. Both the measured and net (or natural) flow rate were obtained while taking into account (i.e., adding) the flow rates from the active pumping station at the time of measurement (listed below). Figure 6 shows the net flow rates calculated in the sections during all of the measurement campaigns. During all of the measurement sessions, the direction of the current in the channel (starting from a distance of 54 m from the pumping station) was always from the station toward the confluence with the Tara River, confirming the fact that the present pumping rates are lower than drainage by the derivation channel and the main spring located in the intake chamber.
For instance, the first values shown in the graph (Figure 6), i.e., those closest to the intake chamber (Section 1), were the outcome of the water flowing from the intercepted spring. An increase in the flow rates, which was observed around 50 m from the intake chamber (section P1-bis), confirmed the presence of a further spring intercepted by the channel, which had already been highlighted by the EC and T measurements. Its flow rate, which was calculated as the difference in flow between this section and the previous one, was, on average, equal to 0.332 m 3 s −1 . During all of the measurement sessions, the direction of the current in the channel (starting from a distance of 54 m from the pumping station) was always from the station toward the confluence with the Tara River, confirming the fact that the present pumping rates are lower than drainage by the derivation channel and the main spring located in the intake chamber.
For instance, the first values shown in the graph (Figure 6), i.e., those closest to the intake chamber (Section 1), were the outcome of the water flowing from the intercepted spring. An increase in the flow rates, which was observed around 50 m from the intake chamber (section P1-bis), confirmed the presence of a further spring intercepted by the channel, which had already been highlighted by the EC and T measurements. Its flow rate, which was calculated as the difference in flow between this section and the previous one, was, on average, equal to 0.332 m 3 s −1 .
Proceeding toward the confluence with the Tara River, around 730 m from the pumping station (Sections 5 and 5-Bis), the flow rates tended to increase up to 1.700-1.900 m 3 s −1 , highlighting a further spring interception along the channel with a contribution of 0.900-1.000 m 3 s −1 . On the other hand, in the terminal branch of the derivation channel, a reduction in flow rate was generally observed, with Proceeding toward the confluence with the Tara River, around 730 m from the pumping station (Sections 5 and 5-Bis), the flow rates tended to increase up to 1.700-1.900 m 3 s −1 , highlighting a further spring interception along the channel with a contribution of 0.900-1.000 m 3 s −1 . On the other hand, in the terminal branch of the derivation channel, a reduction in flow rate was generally observed, with values of 1.100-1.600 m 3 s −1 , depending on the period. The observed reduction in flow rate was linked to the presence of a bypass channel with a siphon that deviated part of the channel discharge laterally. It is interesting to compare these flow measurements to those carried out in the derivation channel during the two-year period from 1980 to 1981 [42] (Table 8), a period in which the uptake from the pumping station varied between 2.400 m 3 s −1 and 3.000 m 3 s −1 . In those withdrawal conditions, the direction of the current was always from the Tara River to the intake chamber (and is therefore indicated with negative values). The net flow estimations closest to the pumping station were carried out at 66 m and 86.5 m from the intake chamber (Sections B and 1A, Table 8), with values between 1.000 m 3 s −1 and 1.200 m 3 s −1 , about 10% higher than the 2016 and 2018 campaigns. Furthermore, considering the extent of the withdrawals in 1980-1981, it can be deduced that a further 1.200-2.000 m 3 s −1 was derived partly from the additional springs intercepted by the channel and partly by the Tara River. In fact, in the 1980-1981 measurements, the flow values measured at the sections near the confluence with the Tara River (indicated as A, 4A, and 4B in Table 8) express a good approximation of the amount of water derived from the natural water course (varying between 0.800 m 3 s −1 and over 1.600 m 3 s −1 ).

Conclusions
The aim of this study was to verify the residual exploitation capacity of a complex of springs fed by a large carbonate aquifer in southern Italy. The complex has been used for more than 60 years for a variety of industrial and irrigation uses. Consequently, the study provides a comparison of the advantages and disadvantages of two streamflow measurement methods in field conditions in which the repeatability of measurements in the same cross-sections was conditioned by apparently random algal blooms.
Through the integration of different monitoring techniques, the survey method adopted allowed us to identify the main source springs that feed the artificial channel. Its current water withdrawal rate is in the range between 0.500 m 3 s −1 and 0.700 m 3 s −1 , and plans are in place to use this water to supply a future desalination plant for drinking use that will bring the total withdrawal to around 1.500 m 3 s −1 .
The thermo-salinity monitoring carried out at different times of the year, both along the derivation channel and downstream of the confluence with the natural watercourse, showed substantial differences due to the presence of other springs at different points of the system, which were fed by circuits characterized by different travel times and salinity. Considering the significant increase in withdrawal from the spring system, we focused on a characterization of the channel's hydrological regime, identifying the most suitable measurement method in terms of simplicity of execution and reliability of the measurements, as well as the distinctive characteristics of the channel (narrow and deep) and the observed flow regimes.
The first flow estimation sessions carried out with flow rate measurements using the CM method and the first outflow estimations using the subsections method posed considerable difficulties due to the presence of algae in the sections under consideration. The F method, although less precise (with an average error of 11%), was easier to implement and repeat over time, so it was chosen to monitor the flow in different sections of the derivation channel and during various periods of the hydrological year. The relative error between the CM and F measurements, including the development of a CM correction method based on the velocity profiles obtained in some vegetation-free cross-sections, was evaluated. Considering the flows taken from the pumping station, it was possible to estimate the natural flow rate regime, an essential element for assessing the sustainability of the pumping rate needed to supply the desalination plant. The current flow rate values from the area close to the confluence with the Tara River were determined on the basis of the measurements carried out over two years of monitoring: between 1.000 m 3 s −1 and 1.500 m 3 s −1 . The present flow regime is different from that observed in the early 1980s, when there was a withdrawal rate between 2.400 m 3 s −1 and 3.500 m 3 s −1 , and a conspicuous amount of water (0.800-1.600 m 3 s −1 ) was derived from the natural water course. The final observation is that this comparison, which is reported in terms of variations in natural flow, showed a slight reduction in terms of current spring flows, probably due to the general increase in groundwater withdrawal for irrigation purposes. It can be deduced that given the considerable decrease in withdrawals at the pumping station, the portion of the spring system drained by the derivation channel is characterized by a flow rate regime, which is both consistent and stable enough to satisfy the increase in withdrawals required for the desalination plant. As an overall benefit of the adopted approach, the combined application of the implemented techniques, although not new, provides a solid, but simple and economic methodology for characterizing springs contributing flow to rivers that can be applied by local water practitioners after limited training. Funding: This research was cofounded by the regional water utility Acquedotto Pugliese (AQP S.p.A.) within an investigative study commissioned for the qualitative and quantitative characterization of the Tara Springs.
Acknowledgments: This study was made possible following an agreement between the Regional Water Utility (AQP), the Irrigation Development Agency (EIPLI), and the Water Research Institute (IRSA-CNR). The authors express their gratitude to EIPLI for providing technical information on the Tara spring pumping system and on the local environment. The support provided by A. Fabiano, L. Ribatti, and P. Lavermicocca during the field measurement sessions and the data preparation is gratefully acknowledged.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.