The Compatibility of Geothermal Power Plants with Groundwater Dependent Ecosystems: The Case of the Cesine Wetland (Southern Italy)

Abstract

The Cesine Wetland, located along the Adriatic coast, was recognized as a Wetland of International Interest and a National Natural Park. Managed by the “World Wide Fund for nature” (WWF), it is considered a groundwater dependent ecosystem which is affected by seawater intrusion. The site was selected to test the environmental compatibility of a low-enthalpy geothermal power plant (closed loop) operating in the aquifer saturated portion with purpose to improving the visitor centre. For this purpose, the long-lasting thermal impact on groundwater was assessed using a multi-methodological approach. The complex aquifer system was carefully studied with geological, hydrogeological and geochemical surveys, including chemical and isotopic laboratory analyses of surface water, groundwater and seawater. The isotopes δ¹⁸O, δD, δ¹¹B, and ³H were useful to clarify the recharge contribution, the water mixing and the water age. All information was used to improve the conceptualization of the water system, including aquifers and the boundary conditions for a density driven numerical groundwater model. The purpose was to forecast anthropogenic thermal groundwater variations up to 10 years of plant working before the plant realization and to validate the solution after some working years. All results show the environmental compatibility notwithstanding the peculiar ecological environment.

1. Introduction

Wetlands are important ecosystems that cover roughly 6% of earth surface [1,2] but show a higher relevance in the hydrological cycle [2] and much higher relevance in terms of biodiversity, being a benchmark of the global environment [3]. The coastal wetlands constitute a critical and vulnerable interface between land and sea environments, which are widely the most productive and ecologically valuable regions of the world [4,5].

The coastal Cesine Wetland, located in the Apulia Region, along the Adriatic coast, represents one of the most famous Italian Groundwater Depending Ecosystem (GDE) (Figure 1).

It was recognized Wetland of International Interest and National Natural Parks. The national authorities have entrusted the management at the “World Wildlife Fund” (WWF), creating a visitor centre in a historical farm.

Due to the visits and staff management needs, technical plants of visitable natural parks could consume relevant amounts of energy, with potential negative impacts on the surrounding environment.

Geothermal resources, classified into low (<90 °C, [6]), intermediate (90–150 °C) and high enthalpy (>150 °C, [6]) by their reservoir temperatures, can provide electric energy generation and/or heat exchange through a variety of applications.

Man has used the thermic properties of the rocks and the thermal resources of the subsoil since ancient times. In more recent times, man learned to use the high temperatures and the hot fluids of some earth surface zones for the production of electrical energy, such as in the case of the Larderello zone in Italy [7], the Geysers area and Salton Sea in United States [8,9].

After the 1973 oil crisis, low enthalpy geothermal energy began to spread in Europe, especially in France and Sweden where the first geothermal probes were installed in 1980 [10,11].

The widespread development of heat pump technology is more recent, pursuing the use of thermal properties of subsoil and groundwater flow for building heating and/or cooling, based on climate and subsoil characteristics. In these cases, the low-enthalpy geothermal energy is valorised, significantly contributing to the reduction of environmental impacts [12,13].

In Italy, the lack of economic incentives for direct use of geothermal heat has delayed the widespread development of low enthalpy solutions [14]. Heat pumps technology started to grow after 2000 [11].

The Apulia Region (Figure 1a) shows a high potential for the application of heat pump technology exploiting groundwater or ground with coupled temperature solutions at shallow depth (10–200 m depth). The Cesine Wetland was selected for an experimental experience, supported by EU, with the purpose to realize a low-enthalpy geothermal power plant (closed loop) in the complex hydrogeological context of a GDE. The scope was to verify the environmental compatibility, the validity of the study approach, assessing the effects on the physical ecological system.

The preliminary study and modelling of the plant and the validation after two working years of the plant were realized in a selected wetland system, a groundwater dependent ecosystem of high value, where the delicate hydrogeological and eco-systemic equilibria, characterized by abiotic and biotic factors, allow the presence of animal and plant species of particular Community importance.

This low-enthalpy geothermal plant is the first one realized in a coastal wetland, which includes a monitoring system for the assessment of the effects of the geothermal plant on the wetland, mainly focusing groundwater temperature variations.

2. Study Area

The Cesine Wetland, located in the Apulia Region along the Adriatic coast (Figure 1a), represents one of the most famous Italian Groundwater Depending Ecosystem (GDE) (Figure 1b,d).

Reclamation activities started at the end of 19th century and finished in 1950. At the beginning of the Fifties, a national law imposed the “Land Reform” establishing relevant modifications in the ownership of the agricultural land. Under this law, roughly 350 hectares of land were entrusted to the Apulia Region whilst the rest, of negligible environmental value, roughly 300 hectares, was given to the reclamation workers to promote their agricultural incomes.

It was recognized Wetland of International Interest in 1971 and became National Natural Park in 1980. The park is managed by WWF.

The protected area extends nowadays over 620 hectares (Figure 1b). It shows a narrow and elongated shape that follows the coastline and is crossed by numerous artificial channels, some of which represent, to the west and southwest, the inland boundary (Figure 1b).

The eastern sector, extending roughly about 350 hectares, is the most environmentally valuable.

It includes the brackish water marshes, wooded areas and the Mediterranean maquis that host the habitats of priority interest from an ecological and environmental point of view.

In this area, there are two perennial coastal lagoons, “Li Salapi” (14 hectares) and “Pantano Grande” (68 hectares), connected by a wide channel (roughly 8 m wide) (Figure 2a). The water lagoon depth is 0.8 m on average, widespread reduced by silting processes up to 0.3 m. The lagoons are connected to marshy areas and to other small retrodunal basins by short dug channels.

In this zone, there is the building “Masseria Cesine” (Figure 1c), an ancient farm used as a visitor centre where the closed loop geothermal plant was designed, preliminary verified for environmental sustainability and realized in 2014. The plant works from September 2014, exchanging heat in the whole depth, down to 200 m below ground surface.

3. Materials and Methods

The detailed geological and hydrogeological study included deep boring and the hydrogeochemical characterization of sea water, surface water, including channels and lagoons, and groundwater. All data concerning existing wells of the study area were registered, improving the stratigraphic and hydrogeological geodatabase together with the realization of piezometric periodic surveys. The topographical survey of wells was performed to improve the accuracy of piezometric head measurements. Boreholes were used for aquifer testing, for temperature logs and for water sampling.

Groundwater temperature logs were carried out before and after the plant construction and activation. The pre-activation thermal logs were carried out both in winter and in summer, without observing relevant thermometric variations, apart from a very narrow depth, for a few meters below the piezometric head of the shallow aquifer.

The aquifer rock thermal properties were determined sampling rock during boring: six laboratory tests were performed. Two in-situ Ground Response Test (GRTs) were also performed using each thermal probe. It consists of a thermal stimulation of geothermal heat exchanger able to measure indirectly the thermal properties of the exchanger and of adjacent soils and rocks, as thermal conductivity, thermal diffusivity and heat capacity [15]. The heat injection and extraction method was used [16].Surface water, groundwater and seawater samples were collected in February 2017 and July 2017 involving six sampling points. The location of water sampling points is shown in Figure 2a.

pH, electrical conductivity (EC), temperature (T), Dissolved oxygen (DO) and redox potential (Eh) of water samples were measured on site with a multiparametric probe together with HCO₃⁻ concentration, determined using volumetric titration with HCl.

The water samples were filtered using a 0.45 μm membrane filter and then collected in double capped 500 mL polyethylene bottles. The samples for cation analysis were acidified by the addition of HNO₃ to pH < 2, while the samples for anion analysis were unacidified. Cations (Li⁺, Na⁺, K⁺, Ca²⁺, Mg²⁺, Sr²⁺) and anions (F⁻, Cl⁻, Br⁻, SO₄²⁻, NO₃⁻) were analysed by means of ion chromatography (IC) methods at the CNR-IRPI Laboratory. The electrical balance was <5% for surface and ground waters. Results on sampling and analyses of the water classification are also shown in Tables 4 and 5 and Figure 4a respectively.

Isotopic determinations concern ³H, δD, δ¹⁸O and δ¹¹B.

³H concentration is the radioactive isotope of hydrogen that decays with a half-life of 12.43 years; it could be useful to assess groundwater age or to assess differences between groundwater paths. It was determined at the CNR-IGG Laboratory (Italy) by means of liquid scintillation spectrometer Quantulus, for counting of samples with a low level of radiation. Analytical error is ±0.3 TU.

Water stable isotopes (δD and δ¹⁸O, reported as permil relative to Vienna Standard Mean Ocean Water) and Boron isotopic composition (δ¹¹B, reported as permil relative NIST SRM951 based on the 11B/10B ratio determination) measurements were carried out at the CIRCE (Centre for Isotopic Research on Cultural and Environmental heritage) laboratory of the Department of Mathematics and Physics of the University of Campania using a TC/EA-CF-IRMS system (Delta V Thermo Fisher; precision was 0.2‰ and 1‰, for δ¹⁸O and δD respectively) and a high resolution multicollector ICP-MS (Thermo Fisher Scientific Neptune Plus; instrumental precision is roughly 0.15‰) respectively.

Samples for δD and δ¹⁸O determinations were collected filling 50 mL high-density polyethylene bottles leaving no headspace to avoid contact with air and horizontally stored at 4 °C [17]. Boron has two stable isotopes, ¹⁰B and ¹¹B with natural abundances of 19.9 and 80.1 atom percent, respectively. Boron is very soluble in natural waters (seawater, groundwater and surface water) and it is a natural constituent of the continental crust. The isotopic composition of dissolved boron in water can be used to determine the natural (e.g., seawater, sea spray and geogenic origin released from rocks and soils) and anthropogenic sources of boron (e.g., fertilizers, detergents) and to investigate the fate of boron co-migrating contaminants (e.g., nitrates, pharmaceuticals). Boron isotopic values are expressed as δ¹¹B permil relative to NIST SRM951 (based on the ¹¹B/¹⁰B ratio determination).

Water for boron isotopic measurements was collected in 2 L polyethylene bottles leaving no headspace. Boron concentration in water was determined in the laboratory by the Azomethine-H spectrophotometric method [18] in order to quantify the volume of sample containing about 50–100 µg of dissolved boron, which is the amount necessary for the isotopic measurement.

The chemical procedure for boron extraction, using only plastic labware, is based on the GAMA method [19]. In particular, all samples were filtered on nylon membrane filters (0.45 µm) to remove particles, adjusted to pH ~ 9 and loaded onto a boron-specific ion exchange resin, Amberlite IRA743 [20]. Before loading the volume of sample containing 50–100 µg of dissolved boron, the resin was crushed (100–200 mesh size), washed and conditioned (Milli-Q water, HNO₃ 3% and NH₄OH 2%) to have a basic pH. Finally, the boron collected on the resin was eluted by nitric acid (3%) and analysed. In each procedural batch, about 20 boron samples, 2 procedural blanks, 2 boric acid standards (NIST SRM951), 2 Sigma Aldrich boron as a Quality Check (QC) sample and 1 unknown sample were prepared.

This procedure involves data calibration through procedural Reference Materials (RMs) and internal standards measurements. Both processed and unprocessed Quality Check (QC) samples were measured in order to evaluate quality of measurement batches by applying Quality Assurance (QA) rules.

Operating with these methods, the hydrogeological conceptualization was completed to support numerical modelling of the thermal plant working. FEFLOW [21] (Finite Element subsurface FLOW system) is a software for simulating groundwater flow with mass and heat transfer in different type of aquifers. It solves the groundwater flow equation in saturated and unsaturated conditions as well as mass and heat transport, including fluid density effects, as particularly useful in case of salt water intrusion. FEFLOW was used representing all the different stratigraphic layers with relative material properties and distributing values of hydraulic, thermal and chemical parameters (water density mainly) in the 3D domain.

A new low enthalpy closed-loop geothermal plant was designed and realized, after the numerical preliminary check of the environmental compatibly.

It consists of a pair of double-U pipe heat exchangers (S1 and S2), bored with a diameter of 0.180 m and deepened until 200 m from the ground surface that intercept three aquifers separated by aquiclude levels with different hydrogeological properties. In each borehole, there is one piezometer tube with a diameter of 0.051 m, screened from 35 to 110 m below ground surface.

Two monitoring wells (W1 and W2) were bored with a diameter of 0.180 m and they reach a depth of 200 m and 25 m respectively. Both the boreholes host two piezometers with a diameter of 0.051 m, screened from 155 to 197 m and from 0 to 25 m below ground surface respectively. Boreholes sections were cemented where unscreened piezometric tubes were installed.

The distance between the geothermal probes is 8.5 m; the distance between the S1 probe and W1 is 6.4 m, while the distance between the S2 probe and W2 is roughly 2.2 m (Figure 1c).

4. Conceptualisation and Discussion

4.1. Geological and Lithological Features

The Cesine Wetland is situated on the eastern part of the Salento Peninsula, the southernmost emerged part of the Apulian carbonate platform, consisting of a thick sequence of Mesozoic carbonate, formed by the Bari and Altamura limestone Formations, covered by deposits from Tertiary to Pleistocene [22].

The geological framework could be conceptualized as a multilayer system of calcarenite and fine-grained limestone, partially affected by significant karst and deep fracturing (Figure 2). Some thin intervals of fine grained soils or massive rocks can be distinguished, with relevant effects on the hydrogeological conceptualisation due to their prevailing impervious nature.

The Mesozoic carbonate basement, which extensively outcrops in the inland area (Altamura limestone Formation), is lowered in the coastal sector due to faults; the Mesozoic carbonate top is roughly 162 m below the ground surface in the study area. Four sedimentary cycles can be distinguished from Tertiary to Quaternary. From the oldest, each cycle corresponds to: the “Pietra Leccese” Formation and the overlying Andrano calcarenite Formation (Miocene); the “Leuca” Formation (Lower Pliocene); the “Uggiano La Chiesa” Formation and the “Trubi” Formation (Pliocene); the Salento calcarenite Formation [23,24,25]. The Pleistocene units are Marine terraced Deposits [26]; they include, somewhere but close to the coast, Aeolian deposits, arranged in continuous dune belts [27]. Recent thin marshy deposits can be observed in the coastal back side of dunes.

The Mesozoic limestone (Altamura limestone Formation) is thick thousands of meters and formed by an alternation of white limestone and dolomitic limestone. This limestone is widespread affected by karstic and fracturing phenomena.

The “Pietra Leccese” Formation is composed by very fine, silty and dark calcarenite which is tens of meters thick. At the bottom, a layer of brown argillite is observed which marks the passage to the Mesozoic limestone.

The Andrano calcarenite Formation is constituted by laminated limestone, marl limestone, grey marl, in which are intercalated shelly levels of gastropods (the passage at the bottom to the “Pietra Leccese” Formation is highlighted by a characteristic level, roughly 4 m thick, rich in fossils, which was observed in the W1 coring).

The Mesozoic and Miocene formations outcrop inland, far from the study area. Above these, an irregular alternation of violaceous limestone, grey silty sand and calcarenite, thick up to 25 m, can be distinguished, corresponding to the “Leuca” Formation.

Moving up to the Pliocene, the deposits of the “Uggiano La Chiesa” Formation are represented by yellowish coarse calcarenite, rich in fossils, in which very fine sand and clay are intercalated. This rock shows generally homogeneous and massive texture and low cohesive clay. The thickness should be less than hundred meters (84 m was the peak value observed in the study area). The “Trubi” Formation is characterized by alternation of violaceous limestone, clayey sand, and calcarenite.

Moving up to the Pleistocene, an alternation of calcarenite sandstone, sand and gravel with clayey matrix, rich in fossils, corresponds to the Marine Terraced Deposits [28]. Outcrops of Aeolian and marshy deposits can be observed close to the coast, over the Marine Terraced Deposits.

The detailed stratigraphy was verified with the continuous coring of the borehole W1, roughly 200 m deep (Figure 2), located very close to the plant location.

Moving from geological features to the lithological characteristics, defined due to boring and in situ observation of outcrops, five lithological unit where distinguished: the fine sand and silty clay unit, including Marshy deposits, part of the “Uggiano La Chiesa” and “Pietra Leccese” Formations; the calcarenite and sand unit, including Aeolian deposits and the rest of the “Uggiano La Chiesa” Formation; the whitish limestone and calcarenite unit, including the Trubi, Leuca and Andrano calcarenite Formations; the fine calcarenite unit, including the main part of the “Pietra Leccese” Formation; and the Mesozoic limestone unit, corresponding to the Altamura Limestone Formation.

4.2. Hydrogeological Conceptualisation

The groundwater flow in the Salento peninsula can be schematized considering two different types of aquifers [29]. The former is represented by a unique aquifer, hundreds of meters thick, which correspond to the intensely fissured and karstic carbonate Mesozoic basement; it extends over the whole Salento Peninsula but does not outcrop in the study area; it is hereinafter called deep aquifer [30,31]. This aquifer is affected by seawater intrusion by the Adriatic Sea to the Ionian Sea [32]. The latter aquifer type is constituted by shallow aquifers: they are spatially discontinuous if the whole peninsula is considered; they are generally outcropping but can be locally overlapped and vertically bounded by impervious strata [30,33].

Five hydrogeological complexes were distinguished on the basis of the geological, lithological and geometrical information together with the observed hydrogeological characteristics and the tests realised during the surveys. They are, from the bottom: limestone, fine calcarenite and clay, calcarenite and limestone, fine sand and clay, calcarenite and sand and silt and clay. As an effect of this conceptualisation, three overlapped aquifers were distinguished in the study area. They are, from the bottom (Figure 2b,c): the confined deep aquifer (DA), the semiconfined intermediate (IA) aquifer and the shallow aquifer (SA). DA corresponds to the aquifer of the carbonate Mesozoic basement, which extends in the whole Salento peninsula.

The limestone hydrogeological complex corresponds to the deep aquifer. It outcrops inland; this portion of the Salento Peninsula is recharged only by direct rainfall infiltration [29]. This aquifer is confined in the study area and show hydraulic conductivity from medium to high fracturing and karst. The DA top, located at 162.5 m from surface in W1, corresponds to the bottom of the fine calcarenite and clay hydrogeological complex. This complex includes fine calcarenite and dark clays which show globally low hydraulic conductivity.

The semiconfined intermediate aquifer corresponds to the calcarenite and limestone hydrogeological complex, due to a lithological succession constituted by calcarenite and whitish limestone, in which secondary sandy and silty levels, able to divide locally the groundwater flow, can be distinguished. The hydraulic conductivity is from low to medium. The IA top corresponds to the bottom of the fine sand and clay hydrogeological complex, roughly 5 m thick, showing very low hydraulic conductivity. IA is generally slightly confined and some tens of meters thick. This aquifer is of negligible relevance as source of usable water. It seems realistic this aquifer lacks direct recharge and receives inflow from the deep aquifer along the upward tectonical boundary, by deep aquifer top and shallow aquifer bottom in decreasing order of relevance.

The groundwater outflow of both the deep and intermediate aquifers happens underwater, and, probably, close to the coastal area.

The shallow aquifer is housed in the calcarenite and sand hydrogeological complex and is phreatic. The hydraulic conductivity, due mainly to primary porosity, is high. The thickness is of some tens of meters. It is recharged by direct rainfall infiltration but exchanges water with the drainage network and steadily feeds the back-dune basins and the wetland system with fresh water.

As all groundwater outflows happen in the wetland or in the very close marine coastal water, the global groundwater quality and quantity affects the ecological equilibria of the wetland ecosystem [34].

Table 1. Mean or statistical values of data measured during surveys from November 2014 to July 2017. K = hydraulic conductivity, TDS = total dissolved solids or salinity.

Aquifer	K [m/s]	TDS [g/L]		Piezometric Head [m asl]		Piezometric Gradient [%]
		Min	Max	Min	Max	Min	Max	Mean
Shallow (SA)	1 × 10−3	0.50	5.00	0.40	1.20	0.164	0.060	0.112
Intermediate (IA)	3 × 10−5	34.10	34.10	0.40	0.96	0.040	0.097	0.060
Deep (DA)	4 × 10−4	35.00	38.00	0.39	1.16	0.039	0.118	0.079

summarizes statistical characteristics of the three aquifers as observed during the surveys conducted from 2014 to 2017.

The preliminary hydrogeological conceptualisation was refined and discussed considering the entire approach, focusing on the hydrogeochemical characteristics of each aquifer/groundwater and superficial water body, including the effect of the seawater intrusion.

4.3. Thermal Rock and Groundwater Features

Thermal conductivity, thermal diffusion and volumetric heat capacity of rocks were determined by laboratory tests using rock samples of W1 (

Table 2. Rock thermal results of laboratory tests.

Sample Depth (m)	Lithological Unit-Sample Description	Thermal Conductivity (W/mK)	Thermal Diffusion (m2/s)	Heat Capacity (MJ/m3K)
14	Calcarenite and sand-yellow tufaceous calcarenite	0.7344	0.4991 × 10−6	1.4715
33	Calcarenite and sand-white tufaceous calcarenite	0.7451	0.4364 × 10−6	1.7075
63	Calcarenite and sand-organic limestone	1.8898	1.0915 × 10−6	1.7314
136	Fine calcarenite-organic calcarenite	1.0703	0.6419 × 10−6	1.6674
157	Fine calcarenite-silty calcarenite with organic matter and fossils	1.1034	0.6527 × 10−6	1.6905
163	Mesozoic limestone-micritic limestone with rudist	1.9771	1.1787 × 10−6	1.6774

). The results of thermal conductivity and thermal diffusion measurements were determined by performing a statistical analysis of 150 measurements for each parameter and sample. The laboratory results are coherent with literature data [16,35,36]. The best thermal properties for heat exchange, thermal conductivity and diffusion, are shown by limestone samples; the heat capacity is little variable between samples.

Moving to the in-situ tests, two GRT tests were realised (

Table 3. GRT tests results. Qa = average discharge rate used during the test, ΔT = thermal variation between probe inflow and outflow (absolute value), λeff = effective thermal conductivity, Rb = hole thermal resistance.

N.	Season	Qa (m3/h)	Mean ΔT (°C)	λeff (W/m/K)	Rb (k/w/m)
1	summer	1.54	6.3	5.4	0.108
2	winter	1.55	4.5	8.4	0.126

). The first GRT was carried out with heat input, simulating the summer plant working. The duration was 72 h; the specific thermal yield was 37 W/m in average, using a total probe power of 7.4 kW. As two identical probes were concurrently used, the output power of probes for the summer test was approximately 14.8 kW (Figure 3).

The second GRT was carried out with heat extraction simulating the winter plant working. The duration was 74 h; the specific thermal yield was −27 W/m, with a total probe power of −5.4 kW. The total output power of probes for winter test was roughly 10.8 kW.

In addition to the results cited in the text, in Table 3 are shown the values of the other parameters derived from the GRT test:

The thermal logs performed before the plant activations were used to define the initial or natural groundwater temperature used as initial conditions for the numerical simulation of heat transport. The schematised initial thermal log corresponds to a homogenous temperature roughly equal to 19.5 °C in the whole surveyed depth. This assumption is coherent with the survey results and the literature experiences. The geothermal gradient is very low in the Apulian region and widely negligible at the considered depths in the Salento Peninsula [37,38], showing the groundwater flow almost cancels the effects of the geothermal gradient [30,39].

Finally, to monitor the thermal perturbation generated by the geothermal probes and to compare the modelling results, thermal logs were realised after the plant activations.

4.4. Geochemical Features

Survey and sampling were performed during wet (February 2016, distinguished with * in

Table 4. Water chemical and isotopic composition. Electrical conductivity (EC) at 25 °C, temperature (T), dissolved oxygen (DO) and redox potential (Eh). (n.d.) not determined; (n.p.) not present; (*) First survey, winter period (February 2016); (**) Second survey, summer period (July 2016); SA, IA and DA for shallow, intermediate and deep aquifers respectively.

Sample Identifier		Water Body	EC	T	pH	DO	Eh	δD	δ18O	3H	B	δ11B
Location	Survey		mS/cm	°C		mg/L	mV	‰VSMOW	‰VSMOW	TU	µg/L	‰NBS951
1	*	SA	1.12	18.13	7.43	0.78	171	−31.88	−6.22	1.2	<50	32.28
2	*	IA	42.4	16.66	6.73	1.63	124	−4.24	−2.45	0.0	2528	39.82
3	*	DA	48.6	18.79	6.77	0.67	−23	1.21	−2.08	0.4	2575	40.99
4	*	channel	1.03	9.00	8.2	6.85	139	−48.25	−8.56	2.4	<50	32.55
5	*	lagoon	8.1	9.85	8.75	3.66	228	−13.99	−3.92	3.3	1785	38.00
6	*	sea	59.1	10.52	7.88	4.23	215	6.22	0.29	1.2	2839	42.34
2	**	SA	0.72	19.23	7.66	1.91	−64	n.d.	n.d.	n.d.	59	34.6
3	**	IA	56.6	21.53	6.79	1.77	12	n.d.	n.d.	n.d.	2527	43.4
3	**	DA	53.4	20.01	6.87	0.77	−30	n.d.	n.d.	n.d.	2521	41.3
4	**	channel	16.3	22.36	7.66	1.05	−96	n.d.	n.d.	n.d.	1128	36.5
5	**	lagoon	15.3	23.13	8.70	6.33	101	n.d.	n.d.	n.d.	1114	37.5
6	**	Sea	60.9	23.68	7.60	6.35	216	n.d.	n.d.	n.d.	2712	42.1

and

Table 5. Chemical composition of water samples. (n.p.) not present; (*) First survey, winter period (February 2016); (**) Second survey, summer period (July 2016); SA, IA and DA for shallow, intermediate and deep aquifers respectively.

Sample Identifier		Water Body	HCO3	F	Cl	Br	NO3	SO4	Li	Na	K	Mg	Ca
Location	Survey		mg/L	mg/L	mg/L	mg/L	mg/L	mg/L	mg/L	mg/L	mg/L	mg/L	mg/L
1	*	SA	475.0	0.21	85.5	0.2	5.2	32.4	0.01	58.8	4.0	13.04	165.6
2	*	IA	250.1	0.00	13980.9	41.8	n.p.	1493.6	0.00	7713.9	217.3	678.3	590.9
3	*	DA	231.8	0.00	17417.5	46.9	n.p.	1892.8	0.32	9301.2	245.2	820.7	466.4
4	*	channel	250.0	0.21	186.5	0.6	0.98	24.8	0.01	125.0	2.9	6.0	80.9
5	*	lagoon	280.6	13.99	2424.8	25.5	65.4	48.4	0.04	1296.6	46.5	132.0	99.4
6	*	sea	146.4	2.38	22378.2	59.0	n.p.	2503.6	0.26	12805.7	340.1	1602.5	231.1
1	**	SA	265.0	0.2	64.0	n.p.	n.p.	13.5	n.p.	48.0	2.1	10.9	58.3
2	**	IA	195.0	1.62	18742.8	n.p.	n.p.	1780.6	n.p.	10135.2	336.5	762.5	708.5
3	**	DA	183.1	0.62	19100.0	62.80	5.45	2131.5	0.12	9766.2	284.9	851.0	694.9
4	**	channel	341.7	n.p.	4817.1	20.8	n.p.	460.9	n.p.	2672.7	100.7	134.7	198.0
5	**	lagoon	281.0	0.7	4470.0	17.7	n.p.	440.9	n.p.	2527.4	98.2	167.7	150.6
6	**	sea	183.2	1.6	23520.3	60.2	n.p.	2405.5	0.06	13568.0	403.03	904.0	120.3

) and dry season (July 2016, distinguished with **).

EC of groundwater ranges between 1.12 and 48.6 mS/cm (at 25 °C) during the wet season survey and between 0.72 and 56.6 mS/cm during the dry season survey (Table 4). Groundwater of intermediate and deep aquifers show the highest EC values, almost close to EC values in the sea and higher values in the dry period, as effect of enhanced seawater intrusion due to the summer well overexploitation. The shallow groundwater is fresh water and generally fresher than surface waters, intermediate and deep groundwater and sea water.

The pH groundwater variability is low between aquifers and seasons; the range is between 6.73 (wet minimum) and 7.66 (dry peak).

The concentration of major ions shows relevant differences between shallow groundwater (samples 1) and the remaining waters (Table 4 and Table 5). Shallow groundwater shows the chemical characteristics of fresh water with highest concentrations of Ca²⁺ and HCO₃^- and can be classified as Ca-HCO3 type. Deep groundwater (samples 3) and intermediate groundwater (samples 2), are characterized by a predominance of alkaline ions (Na⁺-K⁺) and also high content of chloride ion, a typical chemical characteristic of fresh groundwater mixed with by seawater intrusion (Figure 4a,b). Deep groundwater and intermediate groundwater can be classified as Na-Cl rich.

The lagoon samples (location 5) are characterized by high values of EC and high concentrations of alkaline ions (Na⁺-K⁺) and chloride ions and show typical chemical characteristics of fresh groundwater mixed with by seawater intrusion. The effect of seawater intrusion is particularly evident during the dry season (second survey) together with evaporation. The channel samples (4) show the highest variability in terms of EC, from low winter values, which are almost similar to shallow groundwater, to almost high values, which almost similar to lagoon sample. The sea samples (6) were collected tens of meters from the shoreline; the sample 6** shows the highest chloride content, so high to be considered completely free from mixing with continental fresh water.

The geochemical composition of groundwater can be considered to be a mixture of seawater and pure fresh [40] groundwater, as shown by Figure 4 and confirmed by literature. To calculate the mixing ratio with saline water, an inland well of the deep aquifer, was used to sample a representative fresh groundwater (PFG), that was selected as a reference (Well 5 in [41]), to be compared with seawater.

The fraction of seawater (f_sea) in each sample was calculated from the concentration of chloride ions (mmol/L), which is considered a conservative ion in the mixing process [42], using the sample 6** as sea reference:

$$f_{sea}=\frac{m_{C{l}^{-},sample}-m_{C{l}^{-},fresh}}{m_{C{l}^{-},sea}-m_{C{l}^{-},fresh}}$$

(1)

The samples of deep (3) and intermediate (2) aquifers show high rate of mixing, increasing from the wet to the dry season: 74.0% and 81.2% for the former aquifer and 59.4% and 79.7% for the latter, variations due to the seasonal effects of decreasing recharge and increasing well tapping. The shallow groundwater samples (1) show a very low and negligible mixing ratio (not greater than 0.3%).

The samples of channel (4) and lagoon (5) show a mixing ratio with seawater between 0.7% and 20.4% and 10.3% and 19.0% during the rainy and dry season respectively. These results confirm the remarks defined on the basis of EC values from a side and highlight the chemical peculiarities of wetland waters, due to the almost steady outflow of fresh groundwater in the lagoon.

The δ¹⁸O and δD values of groundwater range between −2.08‰ to−6.22‰ and from 1.21‰ to −31.88‰, respectively (Table 4).

Figure 5 shows the δ¹⁸O/δD values compared with the Global Meteoric Water Line (GMWL; [43]) and the Mediterranean Meteoric Water Line (MMWL; [44]). The isotopic composition of water samples was confronted with the value of δ¹⁸O/δD of a rainwater sample (δ¹⁸O = −6.64 and δD = −31.8; [45]) and PFG representative of the deep aquifer, selected as reference (sample 5 in [46]).

All water samples fall between the MMWL and GMWL lines, thus suggesting they are mainly due to meteoric water in a temperate climate, causing infiltration and runoff.

Tritium can be used to determine whether groundwater is modern (roughly less than about 50 years in age for values greater than 1 TU) or pre-modern (older than 50 years in age for values less than 1 TU) [17].

Groundwater tritium ranged from 0.4 to 1.2 TU. The deep groundwater shows very low tritium content (0.4 ± 0.4 TU), which is less than the present seawater tritium (1.1 ± 0.5 TU). The intermediate aquifer shows tritium equal to 0 ± 0.4 TU, a bit lower than deep groundwater value, confirming the slowness of flowing and low recharging path of this aquifer. The highest tritium values (1.2 ± 0.5 TU) were shown by shallow groundwater, confirming the quick recharging process of this aquifer.

These results suggest a very long residence time (>50 years) of the water of the deep and intermediate aquifers (pre-modern groundwater according to [17]) while the water of the shallow aquifer is modern groundwater.

The lagoon and channel waters are characterized by the highest values of tritium (2.4 ± 0.5 TU and 3.3 ± 0.6 TU respectively) suggesting that the contribution of direct rainfall and runoff is a not negligible contribution, especially during the wet season.

According to [47], there are several potential sources of boron in an aquifer and each presents a range of isotopic signatures. Those possible in this study case were boron supplied by seawater intrusion with stable δ¹¹B roughly equal to 39.5‰; boron provided by the modern marine carbonate sediments (δ¹¹B between 14.2‰ and 32.2‰, [48]) and boron supplied by marine aerosols in the recharge water (up to 30‰). In addition, sorption processes can significantly modulate both concentrations and isotopic signature [49].

Shallow groundwater shows low B concentration (location or samples 1); marine aerosols contributing to SA recharge can explain this isotopic signature [50]. In contrast, in high-concentration samples of IA and DA (location 2 and 3), an influence of seawater can be confirmed. The sample of the intermediate aquifer is most influenced of seawater intrusion during the dry season.

The boron concentration values are compatible with the water in the coastal aquifer (Figure 6a). In addition, it is possible to assume that no Sorption-Desorption occurs (Figure 6b).

5. Numerical Model

The main purpose of the model was to define the maximum impact of geothermal plant inside the particular context of the Cesine Wetland protected territory. A 10-year scenario of plant working was defined as a tool, to be numerically simulated, to assess the worst possible effects on the environment from spatial and temporal evolution of the temperature plume. This tool was applied before the plant design to obtain results able to justify the plant operation and after two years of plant working. The results of the second application are here described.

FEFLOW [21] model domain was extended on a 2.4 km² area and a total depth of 300 m. The overall volume was divided in 29 layers, each consisting of 16,700 elements, refining the mesh near the heat exchangers S1 and S2 and the wetland boundaries (Figure 7).

First type hydraulic boundary conditions (Dirichlet BC) were imposed to the southwest or upward limit and to the northeast or downward limit, in order to reproduce a steady flow field based on piezometer gradient and hydrogeological parameters measured during surveys (Table 1 and

Table 6. Numerical model parameters values.

Layer	Hydraulic Conductivity (m/s)	Effective Porosity (-)	Thermal Conductivity of Solid Phase (W/mK)	Volumetric Heat Capacity of Solid Phase (MJ/m3K)	Longitudinal Heat Dispersivity (m)	Transverse Heat Dispersivity (m)
Shallow aquifer	1 × 10−3	0.2	0.7344	1.4715	5	0.5
First separation	1 × 10−8	0.2	0.7344	1.4715	5	0.5
Intermediate aquifer	3 × 10−5	0.2	0.7451–1.8898	1.7075–1.7314	5	0.5
Deep separation	1 × 10−5	0.2	1.0703–1.1034	1.6674–1.6905	5	0.5
Deep aquifer	4 × 10−4	0.2	1.9771	1.6774	5	0.5

). The former and the latter limits roughly correspond to a potentiometric contour line (from 1.12 to 1.56 m a.s.l. as piezometric head of the three aquifers) and to the coastline respectively (0 m a.s.l.). Results of the steady simulation obtained after the calibration were assumed as initial conditions during the transient simulation of the 10-year scenario. The natural daily/seasonal temperature variations observed a few meters below the ground surface were neglected. The effect of submerged areas of the Cesine Wetland were simulated with a simplification that permits to use active cells in these spatial domains: hydraulic conductivity was assigned in these mesh elements a value much higher than the rest (K = 100 m/s), as suggested by USGS [51] and successfully applied in more published study cases (more details are available in [51]). This choice is suggested in the case of low data availability of lake/lagoon bottom, as in our case, and/or to lighten highly the transient calculation simplifying the boundary conditions. This choice seemed particularly advantageous in a portion of the model domain, the outflow area, where this simplification was not able to modify the spatial variability of thermal effects of the plant, without any practical effect upward lagoon, which is at the core of numerical simulation purposes. Due to this simplification, the numerical groundwater outflow happens totally along the coastline instead to happen partially in the lagoon. The reliability of this simplification was verified with simplified numerical tests in steady state conditions.

Boundary conditions for groundwater temperature and salinity were defined, as for piezometric conditions, on the basis of survey results.

The preliminary steady state simulation generated piezometric head, temperature and salinity values in each cell. The steady state results were used for the next transient simulations.

S1 and S2 heat exchangers were simulated, as realised, as closed loop. The specific BHE (Borehole Heat Exchanger) double-U shape module was used, assigning tube diameter, wall thickness and other constructive properties to the 200 m depth boreholes.

In order to keep the most conservative results, the extreme values of operating temperature (

Table 7. Simulated geothermal plant annual activity.

Period	Duration (days)	Daily Activity (h)	Inlet Temperature (°C)	Delta T (°C)
Summer	138	8	50	+30.5
Autumn	58	-	-	-
Winter	110	10	4	−15.5
Spring	59	-	-	-

) were considered for both heat exchangers S1 and S2, ensuring a relevant overestimation of heat exchange effects on groundwater. The working schematisation of Table 7 was applied for each year of the simulation.

In each time step, the numerical results show the maximum temperature variations occur in the shallow aquifer, as confirmed by monitoring activities. It seems the reasonable effect of the highest advection contribution due to the highest hydraulic conductivity, as already observed by previous experiences [52,53].

The numerical results are hereinafter shown and discussed as thermal variations respect to the initial undisturbed groundwater temperature, calculated in steady state conditions (roughly 19.5 °C), mainly focusing on SA. The absolute value of 0.5 °C of thermal variation was considered as benchmark or threshold of a significant thermal variation (a global or generalised definition of threshold or limit of thermal impact does not exist but some Italian regional regulations define 0.5 °C as limit in the case of heat pump systems).

The extreme temperature variations were calculated very close to S1 and S2: the range is from −1.3 to +2.3 °C. The plume of maximum thermal variation was plotted with the latest period of summer activity; it is located to a maximum distance of about 60 m from the heat exchangers, with a lateral width of about 30 m, with a thermal variation range between 0.5 and 2.3 °C (Figure 8).

The geometrical maximum extension (length and width) and distance from the plant of a thermal variation plume and the extreme values of variations were distinguished at the end of each winter season, without any relevant change in the latest six years of the scenario (Figure 9).

The first summer hot plume is still visible at the end of the first winter period, which corresponds to the end of the annual cycle, downgradient of the plant and of the cold or winter plume (Figure 9, top left image).

The comparison of Figure 8 and Figure 9 permits interesting observations. The 10th cold plume is significantly less extended than the 10th warm one: the maximum absolute value of the thermal variation is less of 1.0 °C (the winter minimum is −1.3 °C) and the anomaly is negligible at a distance of 40 m. High is the attenuation of summer thermal variations, due to advection and dispersion in the aquifers. These characteristics are able to reduce the maximum summer hot plumes of previous years (Figure 9, starting from year 2) and cause the disappearance of each annual winter plume. In the former case, the peak variation is reduced from 2.3 to 1.1 °C at roughly 80 m from the site (Figure 9 and Figure 10).

The 10-year simulation shows that temperature variation produced by the geothermal plant, based on the 0.5 °C contour line, is completely attenuated within the maximum distance of 100 m from S1 and S2 borehole heat exchangers (Figure 8, Figure 9 and Figure 10). If the thermal variation of 0.2 °C is considered, the maximum distance by the plant of a thermal plume is roughly 320 m, at a minimum distance from the lagoon of 150 m (Figure 10, bottom left or full view image).

Relative dimensions and magnitudes of generated plumes strongly depend on simulated plant activity cycles, in terms of heat exchanger temperatures and operation time. During the two years of plant working the true plant regime, automatically monitored, was generally very lower than the plant potential maximum, which was applied as steady condition in each season of the modelled 10-year scenario (Table 7). The thermal effects of the working plant shown by thermal logs were very significantly than assessed by modelling, after two years of working and practically negligible.

Vertical cross section views show results on the entire hydrogeological system along the main groundwater flow path (Figure 10). Significant thermal variations at distance from plant were observed in Figure 10b: the temperature increase of the 10th summer is still visible at a maximum distance of 100 m at the end of 10th winter; in this step, the maximum variation of the 9th summer was lower than 0.5 °C (Figure 9, middle right image); the whole effect of previous summers is not significant at a distance greater than 320 m, in any case far from groundwater outflow area.

Heat transfer processes mainly affect the most conductive shallow aquifer that receives the largest amount of the overall thermal impact coming from borehole heat exchangers. Deeper aquifers show a remarkably reduced effect, because S1 and S2 inlet temperature decreases with depth after superficial layer interaction and both the intermediate and the deep aquifers have lower values of hydraulic conductivity. These results show that the thermal plumes extension is highly dependent on aquifer hydraulic conductivity, directly influencing both heat transport velocity and temperature attenuation rate.

The long-term simulation revealed a low impact of the geothermal plant on groundwater system and negligible negative effects on the natural system, as the thermal variation are mainly confined at the shallow aquifer, where thermal effects were calculated at the maximum distance of 320 m from the borehole heat exchangers; if the threshold or limit of thermal impact of 0.5 °C is considered to bound negative geothermal plant effects, this distance decreases at 100 m from the plant.

The resulting thermal affected zone is restricted to a very small portion of the overall domain (Figure 9, bottom left image), notwithstanding the modelled 10-year scenario was defined applying the worst hypothesis of plant regime: in each season, the absolute maximum of heat exchange. This modelling choice guarantees a thermal impact assessment certainly approximated to advantage of caution and justifies the observed thermal variations, very lower than simulated values in each aquifer.

6. Conclusions

The research experience was focused on the definition of a multi-methodological approach to verify the environmental sustainability of a low enthalpy geothermal plant in a GDE of high ecologic relevance, recognized Wetland of International Interest and National Natural Park since 1980.

Due to the intrinsic complexity of the local water system, in which the coastal wetland is due to the overlapped effects of a drainage network and the coastal groundwater outflow of more aquifers, discussion of geological, stratigraphical, hydrogeological and geochemical data was necessary and exhaustive to define the conceptualization of the global water system.

Groundwater salinization due to seawater intrusion is the major process controlling water chemistry of both deep and intermediate aquifer that can be classified as Na-Cl type. The shallow groundwater, classified as Ca-HCO₃ type, is fresh, generally fresher than the other types of water (surface water, other groundwater and seawater).

The isotopic investigation, used to define the hydrological relationships between water bodies, confirmed the existence of separated groundwater paths between the three aquifers, the existence of exchange between aquifers, the groundwater outflow to the wetland system and different sources of water salinization in the study area of the Cesine Wetland. The tritium results shown the existence of pre-modern groundwater with very long residence time (deep and intermediate aquifers) and modern groundwater (shallow aquifer). The tritium results confirmed the hypothesis of negligible direct recharge and a low flow rate for the intermediate aquifer. The δ¹¹B values exhibit two trends, reflecting different B sources. One is clearly related to seawater intrusion, for the deep aquifer and intermediate aquifer, while the second one is related to marine aerosols or to geogenic contribution, for the shallow aquifer. It seems reasonable to consider negligible the role of anthropogenic contributions to the boron content, on the basis of available knowledge.

The obtained results show the importance of the modelling tool for the valuation of the effects induced in the subsoil from the low enthalpy geothermal plant in a complex geological and hydrogeological environment. The results demonstrate also the utility of numerical simulation both in the plant design and in subsequent monitoring phases of its impacts on groundwater, which is at the core of the wetland ecosystem equilibria. The simulation shows, in a scenario after 10 years of operation, a limited thermal perturbation within the natural system, concentrated almost exclusively in the shallow aquifer. In the summer period, the thermal perturbation extends up to a maximum of 60 m from the geothermal probes field, following the main flow direction of shallow aquifer, towards NE, with a lateral expansion of about 30 m. In the winter period, the thermal perturbation completely disappears within 120 m from the source. The thermal effects generated by the geothermal probes are significantly lower in the remaining two aquifers: if the maximum effects are observed in the shallow aquifer, the minimum are shown by the intermediate aquifer.

The extension and maximum thermal variation of the thermal plume, is most affected, in terms of heat transport, by the hydraulic conductivity. It affects both the extension of the thermal perturbation and the attenuation speed of the thermal anomalies generated by the geothermal probes.

In conclusion, the methodological approach used in this study, provides the necessary elements to estimate the thermal impact of a closed-loop geothermal system in a GDE protected area. It can be also used in other protected areas characterized by a high vulnerable environmental equilibrium where it is important to predict even the very low variations of environmental parameters, not only temperature, for their preservation.

Further research efforts will be conducted in test site of the Cesine Wetland to increase the spatial accuracy of the hydrogeological conductivity with new hydrogeological surveys, from a side and to clarify the existence of anthropogenic sources of pollution.