Simulation of coastal groundwater remediation: the case of Nardò fractured aquifer in Southern Italy

Environmental Modelling & Software 21 (2006) 85–97 www.elsevier.com/locate/envsoft

Simulation of coastal groundwater remediation: the case of Nardo` fractured aquifer in Southern Italy Costantino Masciopinto) Consiglio Nazionale delle Ricerche – Istituto di Ricerca sulle Acque, Reparto di Chimica e Tecnologia delle Acque1, Via Francesco De Blasio, 5, 70123 Bari, Italy Received 11 July 2003; received in revised form 21 July 2004; accepted 29 September 2004 Available online 8 January 2005

Abstract A new theoretical approach for evaluating the sharp interface position in a fractured aquifer was applied to the Nardo` aquifer (Southern Italy). The results, based on Dupuit and Ghyben–Herzberg approximations, clearly show both the extent of seawater intrusion and how the latter can be reduced by means of artificial recharge. The new methodology describes a simplified mathematical formulation based on the sharp freshwater–saline interface in a fractured aquifer, rather than solving complicated density dependent partial differential equations of flow and solute transport in porous media. Application of the method requires, in a preliminary stage a knowledge of the freshwater flow rates in the fractured media. This can be derived by means of specific computational codes, under steady conditions. In the second stage, the extent of seawater intrusion can be evaluated in vertical planes of a three-dimensional domain. The methodology was applied to the Nardo` fractured rock aquifer in order to estimate the extent to which seawater had intruded into it. The procedure was easily implemented in a job file of EXCEL Microsoft Office and performed well. Results indicated a good correlation between the coastal spring zones of the model with Landsat Thematic Mapper sensor photos and field investigation. Ó 2004 Elsevier Ltd. All rights reserved. Keywords: Fractured aquifers; Mathematical models; Seawater intrusion; Coastal springs

1. Introduction The intrusion of seawater into coastal aquifers is a common problem in coastal zones of the world where increasing water requirements (Masciopinto et al., 1999) and arid climate have induced over exploitation of groundwater (Park and Aral, 2004). Although many seawater intrusion studies have been carried out both around the world (Oude Essink and Boekelman, 2000; Oude Essink, 2001) and in the Salento Peninsula (Southern Italy) (Cotecchia et al., 1974), very few of ) Tel.: C39 80 5829537; fax: C39 80 5313365. E-mail address: [email protected] 1 National Research Council – Water Research Institute, Chemistry and Technology Section. 1364-8152/$ - see front matter Ó 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.envsoft.2004.09.028

these studies have dealt with mathematical developments of seawater intrusion in fractured aquifers. On the Ionian coast of Southern Italy, a European experimental project (Polemio and Gallicchio, 2003) tested the efficacy of experimental grouting barriers, based on controlled sulphate crystallisation (Gomis-Yagu¨es et al., 1999), in order to stop seawater intrusion into the coastal aquifer. One well-documented system aimed at minimizing seawater intrusion consists in the injection of reclaimed municipal wastewater (Asano, 2002; Masciopinto et al., 1991) into the coastal Salento aquifer (Masciopinto and Carrieri, 2002). However, it should always be remembered that a wide spectrum of technical and health challenges must be carefully evaluated prior to undertaking such a project (Asano and Cotruvo, 2004). For instance, Li et al. (2002)

86

C. Masciopinto / Environmental Modelling & Software 21 (2006) 85–97

Nomenclature bi bm B P

bi

g H(x) Hs i,j K L Lid

Li M Nf n Qi Q0 Qij Reij T Xi

spatially-dependent local fracture aperture [L]; mean aperture of the parallel set of horizontal fractures [L]; maximum saturated aquifer thickness [L]; sum of the all horizontal apertures in the vertical aquifer column with unitary horizontal area (1 ! 1 m2) and thickness B [L]; gravity acceleration (Lt
2); depth of the sharp interface below the sea (i.e., freshwater thickness) [L]; depth (below the sea surface) of the sharp interface at the outflow section [L]; extremity nodes of the generic channel of computational grid (
); hydraulic conductivity of the set of parallel fractures [Lt
1]; distance of coastline from the Ghyben– Herzberg interface toe position [L]; estimated position where the freshwater outflow takes place (i.e., Lid Z L if the intrusion extent is equal to 0) [L]; extent of the seawater intrusion ( Z L
Lid) with respect to the coastline [L]; number of the converging flow rates Qij in every node of the computational grid (
); total number of the fractures of the parallel set (
). aquifer effective porosity (
); flow rate of single fracture of the parallel set per unit of coast length [L2t
1]; freshwater discharged into the sea per unit of coast length [L2t
1]; grid channel flow rate [L2t
1]; hydraulic conductance of grid channel [L2t
1]; hydraulic transmissivity of the parallel set of fractures [L2t
1]; position of the seawater intrusion with respect to y axis [L].

Greek letters d freshwater/(seawater–freshwater) specific weights ratio (
); D4 piezometric head difference [L]; Dx, Dy step size of computational grid [L]; 4, 40, groundwater piezometric head [L]; 4i , 4j gf specific weight of freshwater [ML
2t
2]; gs specific weight of seawater [ML
2t
2]; mf dynamic viscosity of freshwater [ML
1t
1].

observed that watertable accretion enhances beach erosion. However, artificial recharge, as recommended by the Water Framework Directive (ECD, 2000), is a useful supplementary measure in preventing and controlling groundwater pollution. In order to evaluate the extent of seawater intrusion, the theory solves the three-dimensional (3D) flow and transport equations of two density dependent miscible fluids, in saturated porous media (Iribar et al., 1997; Paniconi et al., 2001; Gomis-Yagu¨es et al., 1997). Several computational codes such as SUTRA (Voss et al., 1997), HST3D (Kipp, 1986), SWICHA (Huyakorn et al., 1987), METROPOL (Sauter et al., 1993), FEFLOW (Diersch, 1996) and MOCDENS3D (Konikow et al., 1996) have implemented the above mentioned equation solutions (Putti and Paniconi, 1995) including a code for fractured aquifers (SWIFT (Ward, 1991)). Nevertheless, according to Naji et al. (1999), it is usually very difficult to implement this theoretical approach for practical applications, as there is generally insufficient data available regarding 3D hydrodynamic aquifer parameters (i.e., hydraulic conductivity, boundary and initial conditions, to name but three) and salinity. In addition a 3D representation usually requires very complicated media idealization, with a high number of grid nodes, which involves time-consuming computational procedures. Therefore, the solution of these equations, known in the literature as the Henry solution (Henry, 1964), is usually simplified into a vertical (i.e., two-dimensional) cross-section of the aquifer (Voss et al., 1997; Bear and Bensabat, 2003) with a thickness ranging from 50 to 100 m, which does not give an immediate visualisation of the extent of seawater intrusion. Taking into account these problems, in a recent intrusion study carried out in the Puglia Region (Southern Italy) to evaluate the landward extent of seawater intrusion, scientists (Fidelibus and Tulipano, 1996) preferred experimental monitoring of coastal water quality to the application of mathematical models. The sharp interface approach, on the contrary, facilitates regional scale studies of coastal areas, as suggested by Essaid (1990). Moreover, on a regional scale, the aquifer thickness is minimal compared with the horizontal domain length (usually over kilometres). Naji et al. (1999) suggested a theoretical simplified approach to solve seawater intrusion equations, based on the theory of the stationary position of the sharp seawater–freshwater interface and Boundary Element Methods (Masciopinto et al., 1994). The sharp interface position can easily enhance critical conditions for groundwater well positions (Bower et al., 1999) due to inland seawater movement. A similar theory was addressed here to locate the stationary position of the sharp interface in a fractured aquifer rapidly, where the freshwater flows over the seawater which is considered stagnant. Seawater intrusion in the coastal Nardo`

C. Masciopinto / Environmental Modelling & Software 21 (2006) 85–97

87

(a)

wells

sinkhole

(b)

Sea

(c) Limestone (Cretaceous)

Sandstone (Pleistocene)

Sand (Holocene)

Fig. 1. Study area: Nardo` aquifer (a) in the Salento Peninsula (Southern Italy); (b) geological map at Torre S. Isidoro springs (Magri and Troisi, 1969; Carlin et al., 1973); (c) vertical cross-section.

aquifer of Southern Italy (Fig. 1) was evaluated using a number of vertical planes of the three-dimensional domain. The results of this procedure allowed the extent of the intrusion to be assessed.

2. Theory The fractured aquifer was idealized in a layered model (Essaid, 1990; Kazemi and Gilman, 1993: 272) made by several horizontal fractures bounded by impermeable sediments (Fig. 2). The Nardo` idealization was confirmed by tracer (Rhodamine-Wt and Iodine) tests carried out in the same aquifer, under pumping (Troisi et al., 1985) and natural pressure gradient (Masciopinto et al., 2001). Each fracture is characterized by variable apertures which can be derived from the aquifer transmissivity values estimated in 49 wells of the study area, by mean stochastic methods (Moreno et al., 1988; Dagan, 1989; Brown, 1987). In each fracture, we can define the interface position by considering the hydrostatic equilibrium of freshwater/saltwater pressures, based on the Ghyben–Herzberg equation (see Fig. 2). In addition, using the Dupuit assumption, we assume that inside the fractures freshwater flows in

a horizontal direction (Bear and Verruijt, 1987: 196– 206). In the proposed conceptual scheme all the fractures were assumed to have hydraulic connections between them and to have the same mean aperture bm (mm), whereas the sharp interface approximates the 50% salt concentration contour line in the water. The analytical solution of the interface position can be derived from integration of the Laplace (i.e., the continuity) equation in the vertical plane, where the total freshwater outflow per unit of seacoast length Q0 (m3/s/m) (Bear and Verruijt, 1987: 205), can be derived from the solution of Navier–Stokes’ equation for flow in a single fracture bounded by two parallel plates (Bear, 1972: 164): QðxÞZ


b2m gf vfðxÞ ZconstZQ0 nHðxÞ vx 12 mf

ð1Þ

where bm Z mean aperture of the horizontal fracture (mm); x Z coordinate along the fracture length into the sea direction (m); gf/mf Z density/viscosity ratio of freshwater (Z 107 m
1 s
1 at 20  C);

88

C. Masciopinto / Environmental Modelling & Software 21 (2006) 85–97

z x

Spring at coastline

Lid = available length estimated by computational code (step 1)

Lid Inland spring

Offshore spring

Hs

Set of parallel fractures with the same mean aperture bm

Li =L-Lid

Freshwater

= seawater intrusion extent

Ghyben-Herzberg sharp interface model X= L

L= minimum extent required to avoid seawater intrusion given by Eq. (8) (step 2)

Fig. 2. Idealization of the parallel set of fractures and Ghyben–Herzberg sharp stationary interface under Dupuit assumption (horizontal motion).

f(x) Z piezometric head of freshwater in x direction (m); H(x) Z depth of the sharp interface below the sea (i.e., freshwater thickness) (m); n Z effective aquifer porosity (dimensionless). Based on the analogy of Eq. (1) with Darcy’s formula we can define the hydraulic conductivity of the aquifer as: KZ


b2m gf n 12 mf

ð2Þ

As known, n defines the uniform ratio of the void–space per unit volume of aquifer (Bear and Verruijt, 1987: 19– 30; Borgia et al., 2002) and at section x Z 0 it follows that: Nf P

bi nZiZ1 ð3Þ B P where bi (m) is the sum of all the horizontal apertures in the vertical aquifer column with unitary horizontal area (1 ! 1 m2) and thickness B (m), while Nf is the total number of the fractures of the parallel set. In addition, as all fractures have been assumed to have the same mean aperture bm (mm), it follows that, on average: Nf P

bi Q0 ZNf Qi ZiZ1 Qi bm

ð4Þ

where Qi is the flowrate of the single fracture of the parallel set per unit of coast length. Following the Ghyben–Herzberg formula, in every cross-section at distance x, for 0 % x % L, where L is the distance of coastline from the Ghyben–Herzberg interface toe position (see Fig. 2), we can write the freshwater piezometric head H as: HðxÞZf

gf Zdf gs
gf

ð5Þ

and Eq. (1) can be written as: Q0 vxZ
K

HðxÞ vHðxÞ d

ð6Þ

Eq. (6) is a first order differential equation in x and H, which can easily be integrated between the section x Z 0, at H Z B, and the generic vertical cross-section x at H Z H(x): 2

Q0 xZ
K

ðB2
HðxÞ Þ 2d

ð7Þ

At the outflow section, i.e., at x Z L, we have H(L) Z Hs and: 2

Q0 LZ
K

B2
H2s b2 g ðdf0 Þ
H2s Z
n m f 2d 12 mf 2d

ð8Þ

where Hs (m) is the depth (below the sea surface) of the sharp interface at the outflow section. Obviously, Hs Z 0 can also be set for zero discharge. The Ghyben– Herzberg interface is an approximate theory which

89

C. Masciopinto / Environmental Modelling & Software 21 (2006) 85–97

discretized in a two-dimensional grid. The water flow in each grid element (or channel) of the fracture, with constant aperture bi and bounded by two parallel plates, is governed by the solution of Navier–Stokes’ equation (Brown, 1987):

considers that at outflow section is H(L) Z 0. This is of course not possible, and in the real situation we have an outflow face (see Bear and Verruijt, 1987: 203–205), where the water pressure almost abruptly decreases to zero, at the outflow section. Then, H(L) Z 0, does not mean that we do not have spring outflow but defines where the outflow takes place: i.e., (1) inland; (2) along the coastline or (3) offshore, below the sea surface. Eq. (8) allows the sharp seawater/freshwater interface in the three-dimensional domain of the study area to be defined when the freshwater discharge at the coastline Q0 is known. The reader should note that the origin (x Z 0) of L measurement is defined at the position where f Z f0 (or H Z B).

QZ


1 3 gf b Vf Dy 12 i mf

ð9Þ

where Vf (m) is a function of the difference of piezometric heads (unknown) (fi
fj) at the extremities of the channel with constant length Dx and constant cross-section Dy ! bi in x direction. In order to guarantee the mass conservation at each grid node (without sinks and sources), the following equation can be used: M X

2.1. Freshwater flow equations and numerical formulation

Qij Z

M X fi
fj

1

At the first step of the computational procedure a FORTRAN code is required to estimate the freshwater flow rate in the fractured aquifer under steady conditions (Masciopinto, 1999). The groundwater flow is assumed to take place in a set of parallel and horizontal 2D-fractures, with variable apertures which have the same mean value. The aperture values derive from field transmissivity estimations carried out in the study area (Masciopinto et al., 2001). For each realization of the aperture field, every fracture was

1

Reij

where the converging grid channels, M, are three, instead of four, along the borders of the studied area (Fig. 3) and 1/Reij (m2 s
1) is the conductance of the grid channel i / j (see Fig. 3). This means that no flux boundary conditions are fixed in the borders where piezometric heads not are imposed and that water can enter or leave the aquifer system from the nodes at prescribed heads. Moreover, as the fracture thickness varies abruptly between adjacent nodes, the linear interpolation of the gradient gives a channel conductance,

Qij

Generic border node where piezometric head is not imposed j

Qij

ð10Þ

Z0

i

7600 m

8400 m

Border of domain with fixed head

Generic grid node where piezometric head is not imposed Border of domain with fixed head

Sea

Border of domain with no flux Fig. 3. Diagram of computational scheme used in the FORTRAN code showing the grid nodes of domain subject to mass balance Eq. (10) and the type of boundary condition imposed during simulations.

90

C. Masciopinto / Environmental Modelling & Software 21 (2006) 85–97

1/Reij, that is proportional to the harmonic mean of the nodal values of apertures. By applying Eq. (10) to all grid nodes, a system of equations is obtained; its solution can be found after imposing suitable boundary conditions (Masciopinto, 1999). The piezometric heads in wells and border nodes were used here as Dirichlet boundary conditions with fracture apertures (i.e. flow resistances) and a geometrical description of domain (number of grid lines, borders length and grid size). Other parameters, such as vertical leakages, should be evaluated in advance (Masciopinto et al., 1994) and the computational procedures can only be applied for steady-state conditions. At the coastline the effect of tidal waves on the groundwater flow (Su et al., 2003; Li et al., 2002) was neglected and the periodic seawater level was taken coincident with mean sea level. The Successive Over Relaxation (SOR) method was easily implemented in the code to solve the equation system. System solutions (i.e., the piezometric heads in each grid node) were obtained using a tolerance equal to 10
6 m. The code was successfully tested at the Bari-CNR (Southern Italy) fractured aquifer, under a groundwater radial flow forced by a pumping well, from 4 border wells to the central one (Masciopinto et al., 1997, 2000). The code also showed good performances when compared with MODFLOW (Ruskauff, 1998) results at the experimental site of Bari-Gasometro (Cherubini et al., 2004). The output of the computational code provided flow velocities and flow rates in each fracture domain of the study area. At the initial step of this computation the freshwater model simulation was referred to a zero extent of the seawater intrusion Li (m) (i.e., Lid Z L; see Fig. 2). The model results also produced the groundwater discharges along the coast (Q0) and the contour lines of the piezometric heads. 2.2. Extent of seawater intrusion On the basis of the above calculations, the extent of Lid (see Fig. 2) in x direction can be estimated. Though in the case study Lid had a value close to 3000 m for f0 Z 1 m, on average, the variability of the contour head f Z f0 in the horizontal domain XY usually leads to different Lid values, i.e., one value for each vertical plane ZX. In the same way, Eq. (8) gives a L value for each vertical plane ZX. When the distance Lid is less than that required by Eq. (8), freshwater is discharged into the cross-section located at a distance x Z Lid (!L) (inland spring) (see Fig. 2). Then the extent of seawater intrusion Li can be calculated from the difference L
Lid in each plane ZX. In other words, Li is defined by the difference between the minimum extent required to avoid seawater intrusion according to the Ghyben– Herzberg model and the available extent due to freshwater flow estimated by code. The extent of Li (or L) can be reduced by increasing freshwater discharge Q0

(see Eq. (8)) with artificial groundwater recharges or, by reducing groundwater extraction. On the other hand, negative values of L
Lid, prove the presence of spring outflows, below the sea surface (see Fig. 2). The computational procedure evaluates the length of L
Lid in all planes ZX, at a distance Dy to each other, in order to obtain a three-dimensional representation of the seawater intrusion into the freshwater and the zones with coastal springs (see Fig. 2). Finally, once all locations of the outflow cross-sections along the coast have been evaluated by means of Eq. (8), one last computation is required in order to evaluate the effective freshwater discharged into the sea.

3. The study area The Salento Peninsula consists mainly of sandstone, limestone and dolomite deposits. The Nardo` aquifer is located approximately at 8 km from the Ionian Sea coastline (see Fig. 1). Geological studies (Grassi et al., 1975) of the Nardo` aquifer show that the limestone rock formations are significantly fractured and very permeable. As illustrated in Fig. 1, Pleistocene deposits known as Calcareniti di Gallipoli (sandstone), with spatially variable thickness ranging from 5 to 7 m, are nearest to the ground surface. Below there is a formation with average thickness 30 m, consisting of Calcare di Altamura (limestone) (Upper Cretaceous) intercalated by lenses of terra rossa (calcspar) and loamy sand. The underlying deposits are mainly limestone and dolomite. The watertable is approximately 32 m below the ground surface, close to the sinkhole (see Fig. 1). Freshwater floats over saline water originating from the Ionian Sea due to intrusion and it is confined within the fractured limestone. Fractures are interconnected and partly filled with calcspar or terra rossa. Groundwater flows along preferential pathways within the aquifer under a piezometric head of approximately 3 m above the mean sea level. The preferential pathways are mainly horizontal conduits. The Pleistocene sediments are also fractured but the interconnected fractures are orientated along the local bedding planes. The Nardo` aquifer supplies freshwater to numerous households along the Ionian Sea coastline. Consequently, the water quality of the Nardo` aquifer should be maintained at acceptable levels. The presence of organic compounds (revealed by an oxidation level over 0.5 mg/L as O2) together with an elevated concentration of chloride due to marine intrusion, makes the groundwater unfit for drinking purposes. Maps of the salinity contour line show the improvement of groundwater quality subsequent to injection into the sinkhole (Fig. 4), which was started in 1991. In 1999 the ground water salinity increased slightly with depth below the watertable with values in the range

C. Masciopinto / Environmental Modelling & Software 21 (2006) 85–97

91

1 g/L (1999)

0.5 g/L (1999)

0.5 g/L (1999)

Pendinello Torre S. Isidoro

Oliveto

Torre Termide 12 RF

1 g/L (1969)

Nardò 2 g/L (1969) wells 0

1

2 km

Fig. 4. Salinity contour map (solid lines) in 1969 (i.e., before injection into the sinkhole) reported by Cotecchia (1977) and contour salinities (broken lines) of water samples in wells in 1999 (after 8 years of the injection in the sinkhole).

from 0.2 to 1.5 g/L. This information, together with water temperature, pH and dissolved oxygen, was obtained in a well by means of the Ocean Service (Idronaut S.r.l.) probe, which was able to draw these water constituent profiles up to 50 m below the watertable (Fig. 5). The sinkhole has been loaded in the past 13 years, by the Asso Channel which outflows municipal effluent from treatment plants and local surface drainages (Balice et al., 1989) at a constant rate of about 140 L/s. In the coastal aquifer, the groundwater was sampled in monitoring wells located at different distances from the sinkhole, before and after 1991 when injection started (Table 1). The water samples were collected after ten minutes of pumping to ensure the removal of stagnant water in the pipelines. In wells without pumps, the water was collected by means of samplers located at 0.5–1 m below the watertable. The analytical methods utilized are reported in the Standard Methods (1995) and CNR-IRSA methods (1994). Ten wells (3 C 7) were monitored from November 1953 to February 1954, 8 wells (6 C 2) during winter 1979/1980 and 14 wells (7 C 7) from November 1998 to March 1999. Water sampled from these wells confirmed a reduction of groundwater salinity owing to the injection of treated wastewater (see Table 1). However, a general improvement in groundwater quality due to the reduction of

ammonia and nitrites, was also observed (Fig. 6) coincident with increased nitrates, chemical oxygen demand (COD) and Total dissolved Organic Carbon (TOC) in the wells near the sinkhole (see Table 1). Due to the presence of dissolved organic compounds and pathogens (Masciopinto et al., 2003), the groundwater requires treatment to make it safe for drinking. 3.1. Simulation results The computational procedure was applied in a threedimensional XYZ domain 8400 m ! 7600 m ! 30 m of Nardo` (Lecce) fractured aquifer. In this area the groundwater flow is affected by withdrawals for drinking and irrigation purposes. The study of groundwater flow refers to the years 1980 and 1999. Due to the arid climate, precipitation usually occurs consecutively for a few days (5–6) and 3–4 times a year during winter/spring seasons. In these periods the runoff abruptly increases the Asso Channel flowrate to 450 L/s. The model, however, simulated the steady conditions which occurred during a year in more frequent drought periods. Moreover a constant recharge flowrate of 140 L/s (or 4.4 million of m3/y) into the sinkhole was considered in the year 1999. The piezometric heads measured in the Nardo` wells in 1999 were imposed as boundary conditions. As the groundwater

92

O2 (mg/L)

EC (µS/cm)

15.7

15.6

15.4

15.3

15.2

15.1

15

0

14.9

7

6

5

4

3

2

0

1

1200

1000

800

600

400

200

0

C. Masciopinto / Environmental Modelling & Software 21 (2006) 85–97

T (°C)

10

20

30

(m)

O2 (mg/L)

EC (µS/cm)

17.725

17.72

17.7105

17.71

17.705

17.7

17.695

12

10

8

6

4

2

0

2500

2000

1500

1000

(a)

500

0

40

T (°C)

1

2

3

4

5

6

(m)

(b)

Fig. 5. Vertical profiles (Ocean Service probe) into Oliveto well, 600 m from the sinkhole (a) (and 7600 from the coast) and Torre Termide 12 RF well (b) at 5800 m from the sinkhole (2600 m from the coast), during winter 1999.

Table 1 Ground water constituents (mean of the average values in every well) during winter in wells near the coast (3–4 km) and the injection sinkhole (7–8 km from the coastline) before and after implementation of the artificial recharge in 1991 N-NH3 (mg/L)

N-NO
2 (mg/L)

Area 1 – distance from injection: 3–4 km (7–8 km from the coastline) 1953/1954 3 2093 7.5 1979/1980 6 2127 7.1 12.3 1998/1999 7 1418 7.0 17.4

0.1 0

0.079 0.002

Area 2 – distance from the injection: 7–8 km (3–4 km from the coastline) 1953/1954 7 2352 7.4 1979/1980 2 2820 6.8 23.1 1998/1999 7 1633 7.5 20.4

0.1 0

0.14 0

Sampling period

No. of wells

Electrical conductivity (mS/cm)

pH

COD (mg/L as O2)

N-NO
3 (mg/L)

NaCl (mg/L)

TOC (mg/L)

2.2 5.7

513 506 396

4.8 6.7

1.9 10.2

584 711 461

4.5 5.9

C. Masciopinto / Environmental Modelling & Software 21 (2006) 85–97

3500

µS/cm

93

Torre Termide 12 RF Pendinello

3000 2500

mg/L

2000

100

1980 1999

1500

10

1000 1

500 0

Torre Termide 12 RF

0.1 1954

1980

Pendinello

1999

Fig. 6. Water Electrical Conductivity (EC) (on the left) and N-NO3 concentration (on the right) vs. time in two wells (Torre Termide 12 RF and Pendinello) at respectively, 2600 m and 3200 m from the coast (i.e., 5400 and 5800 m from the sinkhole).

salinity in the wells (see Figs. 4 and 5) ranged form 0.2 to 1.5 g/L, the water specific weight variations (!2&) did not affect piezometric measurements. Prescribed piezometric heads were also imposed along the border nodes of domain in both the 1999 and 1980 simulations (Fig. 7). The piezometric head in the grid node which simulated the injection sink was derived from the analytical solution of the equation for a steady radial flow to a well in a confined aquifer, under a constant injection rate of 140 L/s. Pumping tests carried out in the study area showed that the aquifer transmissivity ranged from 0.09 to 0.005 m2/s. The average aperture of the y

parallel set bm (1.05 mm) was estimated by using an average value of 0.0475 m2/s (Masciopinto et al., 1997). In addition, tracer tests allowed the total fracture apertures in the verticalPcross-section of the saturated freshwater aquifer, i.e., bi Z 0.04 m, to be evaluated. The latter value suggests a total number of vertical fractures of the coastal aquifer affected by freshwater flow equal to 40, on average, if all the fractures are considered as having the same aperture of one mm. The computational code results showed a watertable mound due to artificial recharge, which has an extension of 3 km around the sinkhole (see Fig. 7). The average rise in the

x

0

3.5

Sinkhole

3.0 0 140 L/s

2.30

Asso Channel

Ionian Sea

1.50

0 0.1

00

2.

1.00

1.50

0.50

1.00

0.60 0.40 0.30

0.20

N

NARDO’

Scaled vectors proportional to water velocities (1999) Piezometric contour heads: 1980 1999

0 Monitoring wells

1

2 km

Fig. 7. Piezometric heads (in meters above the sea) before (1980) and after (1999) the artificial recharge of 140 L/s started in 1991 and the 1999 groundwater velocities.

94

C. Masciopinto / Environmental Modelling & Software 21 (2006) 85–97

piezometric head was almost 1.5 m, in a region of about 20 km2, around the sinkhole, with respect to the piezometric heads before the injection reported by PRA (1983).

freshwater velocity was 20–50 m/d (with higher values close to the sinkhole) and the reduction in the extent of seawater intrusion was 2 km, on average. This reduction increased the groundwater availability by 0.8 ! 106 m3/yr due to the volume of groundwater (Z 2000 ! 8000 ! 10 ! 0.005, where 0.005 is the aquifer storativity of the average 10 m saturated thickness of the aquifer near the coast). Low salinity was restored in 20% of the total groundwater stored in the subsoil. Moreover, the recharge directly increased the storage by 1.2 ! 106 m3/yr, which corresponds to the portion of the freshwater injected which was not discharged into the sea. These effects were confirmed by a reduction in the salinity of sampled water in wells monitored in comparison with values monitored before recharge (see Fig. 6).

3.2. Reduction in the extent of seawater intrusion Artificial recharge has increased the piezometric heads and reduced seawater intrusion. By using the above mentioned computational procedure (see Section 2), the extent of seawater intrusion was evaluated both with and without artificial recharge. Eq. (8) was applied to each of the 39 vertical planes, at distances of 200 m. The methodology was implemented in an EXCEL Microsoft Office Premium (2000) file (Fig. 8) in order to draw the intrusion line development for the entire region. Due to variability of the estimated freshwater outflow along the coast, the extent of seawater intrusion is not uniform and shows a smoothed fingered shape, along the coast line (see Figs. 8 and 9). The calculations were carried out by using a fixed ratio of B/f0 equal to 30 (Cotecchia et al., 1974). With injection, the total freshwater outflow in the Ionian Sea was increased to 208 L/s from an initial value of 100 L/s. The average

….. ……

……

..

4. Validation of the theory for the saltwater/ freshwater interface placement The shortage of salinity data in the study area (see Fig. 4) did not allow representative contour lines of groundwater salinity to be interpolated, in order to

…

……. ……

(a) y

8000

Li (1999) m Li (1980) m

7000 6000 5000

(b)

4000 3000 2000

5000

(c)

4000 3000 2000

Xcoast

1000

0 -2000

Xi (1999) m Xi (1980) m

6000

1000 -5000

y

7000

Distance parallel y (m)

Distance parallel with the coastline (m)

8000

1000

x 4000

Distance from coastline (m)

-5000

-2000

0

1000

x 4000

Distance from y (m)

Fig. 8. Microsoft Excel Worksheet showing (a) the computation of the intrusion extent Li, (b) its plot with respect to the coastline by using Eq. (8) (see Sections 3.1 and 3.2 for values) and (c) its transformation to the position Xi by means of the coordinates of the coastline (Xcoast).

95

C. Masciopinto / Environmental Modelling & Software 21 (2006) 85–97

Torre S. Isidoro 1 1

1

2 2

3

3

4

4 4 5 6

5

5

(a)

6 6

1980

1999

(b)

(c) y

0

3 km

o

x

Fig. 9. Comparison of the freshwater outflow zones at Torre S. Isidoro: (a) springs monitored by Cotecchia (1977: 145–158); (b) springs visualized by Landsat Thematic Mapper photo (1997); (c) spring zones defined by computational procedure (1980 and 1999).

define an accurate 50% saltwater contour line along the coastline and to validate the theory here developed. Indeed as mentioned before (see section 1), the availability of such data is not very frequent due to the complexity of measuring and collecting information (Melloul and Goldenberg, 1997). Despite data shortage, the validation of the simulated placement of the seawater/freshwater interface at Nardo` aquifer was carried out by comparing the coastline images and studies that monitored the outflow of springs present on the Torre S. Isidoro coast (Fig. 9). According to some authors (Magri and Troisi, 1969; Cotecchia, 1977) the Torre S. Isidoro springs have flowrates of up to 20 L/s (see zone 1, in Fig. 9). These coastal spring outflows can also be pinpointed examining a photo taken of the area by a Landsat Thematic Mapper (TM) sensor (Uricchio and Masciopinto, 1999) thanks to the different reflectance of freshwater with respect to saltwater (Leone and Menenti, 2003) (Fig. 9b). The satellite image taken in June 1997 shows a possible correlation between spring zones 1, 2, 3, 4 and 6, and the results of the computational procedure used, referred to 1999. On the contrary, zone 5 does not show a correlation, even though Carlin et al. (1973) monitored springs in this zone (Fig. 9a). This is probably due to the depth of the freshwater discharges below the sea surface, which prevented the landsat image from detecting any changes in color due to the freshwater discharges. In this zone, following the proposed procedure, the outflow crosssection is 3 km offshore. On the contrary the absence of the coastal springs in zone 2 (Fig. 9a) agrees with the

inland placement of the interface given by model simulation referred to 1980 (see Fig. 9c).

5. Conclusions A new theoretical approach has been applied in the Nardo` coastal aquifer in order to define the stationary seawater–freshwater interface in fractured media and estimate by how much the artificial injection of 140 L/s has reduced the extent of seawater intrusion since 1991. The computational procedure allowed a quasi threedimensional visualization of seawater intrusion into freshwater in the whole region to be obtained. The latter was validated by means of satellite images and the results of spring outflow monitoring in a previous study. Almost all coastal zones (zones 1, 2, 3, 4 and 6) showed possible correlations, enhanced by means of computational procedure, with the spring outflows on the coast. In zone 5, the depth of the spring below the sea prevented it from being visible in the Landsat TM image. The methodology applied allows a rapid placement of the seawater/freshwater interface, following Ghyben–Herzerg and Dupuit approximations, over a large coastal area, when the coastal groundwater discharge and piezometric heads are known. The methodology, applied to the Nardo` aquifer, estimated that the sinkhole allowed groundwater remediation, on average, of 20% of the groundwater volume. This was due to a reduction in the extent of seawater intrusion and yielded an increase of 1.2 ! 106 m3/yr in the total

96

C. Masciopinto / Environmental Modelling & Software 21 (2006) 85–97

water stored. These volumes could be very useful for irrigation purposes, although treatment of pumped water from wells is required, before drinking. The improvement achieved in the Nardo` area should encourage the local Authority to implement new recharge projects into these coastal and arid regions. In which case, precautionary disinfection of pumped water must be carried out in order to avoid outbreaks of waterborne diseases due to incomplete pathogen inactivation in the subsoil.

Acknowledgments Thanks are due to anonymous reviewers for constructive comments which were very useful in improving the descriptions of theory and results of the computational procedure reported in the manuscript.

References Asano, T., 2002. Multiple uses of water: reclamation and reuse. GAIA 11 (4), http://www.gaia-online.net/. Asano, T., Cotruvo, J., 2004. Groundwater recharge with reclaimed municipal wastewater: health and regulatory considerations. Water Research 38, 1941–1951. Balice, V., Vurro, M., Di Fazio, A., Soria, B., 1989. Qualita` dell’acqua sotterranea in una zona della penisola salentina. Ingegneria Ambientale 18 (2), 93–98. Bear, J., 1972. Dynamics of Fluids in Porous Media, vol. 126. American Elsevier Publishing Company, Inc., New York, pp. 164–165. Bear, J., Bensabat, J., 2003. Analysis of grouting performance. Proceedings of the IV Meeting on Crystallization Technology for Prevention of Salt Water Intrusion, September 26–29, 2002 – Scanzano Joinico (Milano, Italy). GEO-GRAPH S.n.c. Publishing Company, Segrate, Matera, Italy, pp. 32–38. Bear, J., Verruijt, A., 1987. Modeling Groundwater Flow and Pollution, D. Reidel Publishing, Dordrecht, Holland, pp. 196–206. Borgia, G.C., Bortolotti, V., Masciopinto, C., 2002. Valutazione del contributo della porosita` effettiva alla trasmissivita` di acquiferi fratturati con tecniche di laboratorio e di campo. IGEA-Groundwater Geoengineering 17, 31–44. Bower, J.W., Motz, L.H., Duden, D.W., 1999. Analytical solution for determining the critical conditions of saltwater upconing in a leaky artesian aquifer. Journal of Hydrology 221, 43–54. Brown, S.R., 1987. Fluid flow through rock joints: the effect of surface roughness. Journal of Geophysical Research 92 (B2), 1337–1347. Carlin, F., Magri, G., Mongelli, F., 1973. Temperatura delle acque sotterranee della penisola salentina. Geologia Applicata e Idrogeologia III (part II), 155–165. CNR-IRSA, 1994. Metodi Analitici per le Acque, Quad. IRSA, 100, Roma (I). Cherubini, C., Giasi, C.I., Masciopinto, C. Stochastic fracture flow modeling and influence of the computational procedures on the pollutants migration in a case study. Proceedings of the International Conference on Finite Element Model, MODFLOW and More 2003 Solving Groundwater Problems. Karlovy Vary (Czech Republic), 13–16 September 2004, 189–192. Cotecchia, V., Tadolini, T., Tulipano, L., 1974. The results of research carried out on diffusion zone between fresh water and seawater

intruding the land mass of Salentine Peninsula (Southern Italy). Proceedings of the International Symposium on Hydrogeology of Volcanic Rocks. Isole Canarie, Lanzarote, Spain, pp. 1–15. Cotecchia, V., 1977. Studies and investigations on Apulian groundwater and intruding seawaters (Salento Peninsula). Quaderni IRSA 20, 145–158, Rome (I). Dagan, G., 1989. Flow and Transport in Porous Formations. Springer-Verlag, Berlin, Germany, pp. 465. Diersch, H.J.G., 1996. Interactive, graphics-based finite-element simulation system FEFLOW for modeling groundwater flow, contaminant mass and heat transport processes. User’s Manual v. 4.5, Institute for Water Resources Planning and System Research, Ltd., The Netherlands. Essaid, H.I., 1990. A multilayered sharp interface model of coupled freshwater and saltwater flow in coastal systems: model development and application. Water Resources Research 26 (7), 1431–1454. European Community Directive (ECD), 2000. Water Framework Directive 2000/60. European Parliament and of the Council of 23 October 2000. Establishing a frameworks for Community action in the field of water policy. Official Journal of the European Communities, (22-12-2000). EXCEL Microsoft Office Premium, 2000. Microsoft Corporation. http://support.microsoft.com/. Fidelibus, M.D., Tulipano, L., 1996. Regional flow intruding sea in the carbonate aquifers of Apulia (Southern Italy). Proceedings of the 14th Salt Water Intrusion Meeting, Malrno¨, Sweden, pp. 230–240. Gomis-Yagu¨es, V., Boluda-Botella, N., Ruvia-Bevia´, F., 1997. Column displacement experiments to validate hydrogeochemical models of seawater intrusion. Journal of Contaminant Hydrology 29, 81–91. Gomis-Yagu¨es, V., Boluda-Botella, N., Ruvia-Bevia´, F., 1999. Gypsum precipitation dissolution as an explanation of sulphate concentration during seawater intrusion. Journal of Hydrology 228, 48–55. Grassi, D., Tadolini, T., Tulipano, L., 1975. Influenza delle caratteristiche morfologico - strutturali e paleogeografiche sull’idrologia della zona situata a nord di Otranto (penisola salentina). Atti del 3  Convegno Internazionale Sulle Acque Sotterranee, Palermo, 1–5 settembre 1975. Henry, H.R., 1964. Effects of dispersion on salt encroachment in coastal aquifers: in Sea Water in Coastal Aquifers, U.S. Geological Survey Water Supply Paper 1613-C, C71–C84. Huyakorn, P.S., Andersen, P.F., Mercer, J.W., White Jr., H.O., 1987. Saltwater intrusion in aquifers: developments and testing of a three-dimensional finite element model. Water Resources Research 23 (2), 293–312. Iribar, V., Carrera, J., Custodio, E., Medina, A., 1997. Inverse modelling of seawater intrusion in the Llobregat delta deep aquifer. Journal of Hydrology 198, 226–244. Kazemi, H., Gilman, J.R., 1993. Multiphase flow in fractured petroleum reservoirs. In: Bear, J., Tsang, Chin-Fu, De Marsily, G. (Eds.), Flow and Contaminant Transport in Fractured Rock, vols. 6–19. Academic Press, Inc., London, pp. 270–272. Kipp, K.L., 1986. HST3D, a computer code for simulation of heat and solute transport in three-dimensional groundwater flow systems. U.S.G.S. Water Resources Investigations Report 86-4095. Konikow, L.F., Goode, D.J., Hornberg, G.Z., 1996. A three dimensional method of characteristics solute transport model (MOC3D). U.S.G.S. Water Investigations Report 96-4267. Leone, A., Menenti, M., 2003. http://www.inea.it/otris/salinita/leone. htm, July 9. Li, L., Barry, D.A., Pattiaratchi, C.B., Masselink, G., 2002. BeachWin: modelling groundwater effects on swash sediment transport and beach profile changes. Environmental Modelling & Software 17, 313–320.

C. Masciopinto / Environmental Modelling & Software 21 (2006) 85–97 Magri, G., Troisi, S., 1969. Sull’influenza delle fluttuazioni di specchi d’acqua sui livelli delle falde costiere. Applicazioni allo studio della circolazione idrica sotterranea nella penisola salentina. Geologia Applicata e Idrogeologia IV, 25–40. Masciopinto, C., Palmisano, N., Tangorra, F., Vurro, M., 1991. A decision support system for artificial recharge plant. Water Science and Technology 24 (9), 331–342. Masciopinto, C., Passerella, G., Vurro, M., Castellano, L., 1994. Numerical simulations for the evaluation of the free surface history in porous media. Comparison between two different approaches. Advance in Engineering Software 21, 149–157. Masciopinto, C., Di Fazio, A., Benedini, M., Troisi, S., 1997. Evaluation of aquifer parameters by means of field tests and using different conceptual models. Proceedings of the XXVII IAHR Congress, San Francisco, California, USA, August 1997, pp. 493–498. Masciopinto, C., Barbiero, G., Benedini, M., September 1999. A large scale study for drinking water requirements in the Po basin (Italy). Water International 24 (3), 211–220. Masciopinto, C., 1999. Particles transport in a single fracture under variable flow regimes. Advances in Engineering Software 35 (5), 327–337. Masciopinto, C., Benedini, M., Troisi, S., Straface, S., 2000. Conceptual models and field test results in porous and fractured media. In: Katsifarakis, K.L. (Ed.), Groundwater Pollution Control. WIT Press, Boston, pp. 245–270. Masciopinto, C., Di Fazio, A., Troisi, S., 2001. Relation between experimental data and forecasting models for fractured aquifers. Proceedings of the Second International Symposium on Protection of Groundwater from Pollution and Seawater Intrusion, September 27–October 1, Bari (I). Futur Grafica Italia, Bari, pp. 327–340. Masciopinto, C., Carrieri, C., 2002. Assessment of water quality after 10 years of reclaimed water injection: the Nardo` fractured aquifer (Southern Italy). Ground Water Monitoring & Remediation 22 (1), 88–97. Masciopinto, C., La Mantia, R., Jatta, E., Iatta, R., 2003. Efficienza della filtrazione naturale per approvvigionamento idrico da falde fratturate del Salento. L’ACQUA 3, 23–32. Melloul, A.J., Goldenberg, L.C., 1997. Monitoring of seawater intrusion in coastal aquifer: basic and local concerns. Journal of Environmental Management 51, 73–86. Moreno, L., Tsang, Y.W., Tsang, C.F., Hale, V., Neretnieks, I., 1988. Flow and tracer transport in a single fracture: a stochastic model and its relation to some field observations. Water Resources Research 24, 2033–2048. Naji, A., Chend, A.H.-D., Ouazar, D., 1999. BEM solution of stochastic seawater intrusion problem. Engineering Analysis with Boundary Elements 23, 529–537.

97

Oude Essink, G.H.P., Boekelman, R.H., 2000. Saltwater intrusion in coastal aquifers. In: Katsifarakis, K.L. (Ed.), Groundwater Pollution Control. WIT Press, Boston, pp. 145–201. Oude Essink, G.H.P., 2001. Improving freshwater supply – problems and solutions. Ocean & Coastal Management 44, 429–449. Paniconi, C., Klhaifi, I., Giacomelli, A., Tarhouni, J., 2001. A modelling study of seawater intrusion in the Korba coastal plain, Tunisia. Physics and Chemistry of the Earth (B) Hydrogeology, Oceans and Atmosphere 26 (4), 345–351. Park, C.-H., Aral, M.M., 2004. A multi-objective optimisation of pumping rates and well placement in coastal aquifers. Journal of Hydrology 290, 80–99. Piano Regionale di Risanamento delle Acque (PRA), 1983. Bollettini della Regione Puglia, Bari. Polemio, M., Gallicchio, G., 2003. Proceedings of the IV Meeting on Crystallization Technology for Prevention of Salt Water Intrusion, September 26–29, 2002 – Scanzano Joinico (Matera, Italy). GEO-GRAPH S.n.c. Publishing Company, Segrate, Milano, Italy, pp. 75. Putti, M., Paniconi, C., 1995. Picard and Newton linearization for the couplet model of saltwater intrusion in aquifers. Advances in Water Resources 18 (3), 159–170. Ruskauff, G.J., 1998. Stochastic MODFLOW for Monte Carlo Simulation. Guide to using in GWVISTAS v. 4. Environmental Simulations Team: J.O. Rumbaugh and D.B. Rumbaugh, INC., Shrewsbury, UK, pp. 42–46. Sauter, F.J., Leijnse, A., Beusen, A.H.W., 1993. METROPOL’s User’s Guide, Report Number 725205.003. National Institute of Public Health Protection, Bilthoven, The Netherlands. Standard Methods for the Examination of Water and Wastewater, 19, 1995. American Public Health Association, 1015 Fifteenth Street, NW Washington, DC, 20005, 9-99, 9-102. Su, N., Liu, F., Anh, V., 2003. Tides as phase-modulated waves inducing periodic groundwater flow in coastal aquifers overlaying a sloping impervious base. Environmental Modelling & Software 18, 937–942. Troisi, S., Vurro, M., Di Fazio, A., 1985. Studio metodologico di parametri idrodisperivi mediante misure in situ. Quaderni IRSA 69, Rome, Italy. Uricchio, V., Masciopinto, C., 1999. Riconoscimento semiautomatico di linee preferenziali di scorrimento d’acqua da immagini satellitari. Acque Sotterranee 3, 63–71. Voss, C.I., Boldt, D., Shapiro, A.M., 1997. A graphicaluser interface for the U.S. Geological Survey’s SUTRA code using Argus ONE. U.S.G.S. Open-file Report 97-421, Reston, Virginia. Ward, D.S., 1991. Data Input for SWIFT/386, version 2.50, Geotrans Technical Report, Sterling, Virginia.