Numerical modeling of flow and transport in the Bari industrial area by means of rough walled parallel plate and random walk models

Modeling fluid flow and solute transport dynamics in fractured karst aquifers is one of the most challenging tasks in hydrogeology.

The present study investigates the hotspots of groundwater contamination in the industrial area of Modugno (Bari – southern Italy), where the limestone aquifer has a fractured and karstic nature.

A rough walled parallel plate model coupled with a geostatistical analysis to infer the values of the equivalent aperture has been implemented and calibrated on the basis of piezometric data. Using the random walk theory, the steady-state distribution of hypothetical contamination with the source at the hotspot has been carried out, reproducing a pollution scenario which is compatible with the observed one. From an analysis of the flow and transport pattern it is possible to infer that the anticline affecting the Calcare di Bari formation in direction ENE–WSW influences the direction of flow as well as the propagation of the contaminant.

The results also show that the presence of nonlinear flow influences advection, in that it leads to a delay in solute transport with respect to the linear flow assumption. This is due to the non-constant distribution of solutes according to different pathways for fractured media which is related to the flow rate.

1 Introduction

The characterization and the description of phenomena that involve fractured aquifers, especially when considered in relation to water resource exploitation, are important issues because fractured aquifers serve as the primary source of drinking water for many areas of the world. In fractured rock aquifers, groundwater is stored in the fractures, joints, bedding planes, and cavities of the rock mass. The ability of a fracture to transmit water as well as contaminants depends primarily on the size of the opening, or the fracture aperture.

The parallel plate model is widely used to simulate flow in a fracture due to its simplicity in idealizing a fracture. Many researchers (Baker, 1955; Huitt, 1956; Snow, 1965, 1970; Gale, 1977) have used flow between smooth parallel plates as a model for flow in fractures. The solution to the Navier–Stokes equation for flow between parallel plates, known as plane Poiseuille flow, has been known to fluid mechanicians since the nineteenth century.

Witherspoon et al. (1980) and Elliott and Brown (1988) suggested that a factor should be introduced into the parallel plate theory to take account of the effects of joint surface properties.

Zimmermann and Bodvarsson (1996) discussed the problem of fluid flow through a rock fracture within the context of fluid mechanics. The derivation of the “cubic law” was given as the solution to the Navier–Stokes equations for flow between smooth, parallel plates. They analyzed the various geometric and kinematic conditions that are necessary in order for the Navier–Stokes equations to be replaced by the more tractable lubrication or HeleShaw equations and reviewed various analytical and numerical results pertaining to the problem of relating the effective hydraulic aperture to the statistics of the aperture distribution.

Some researchers proposed a variable aperture model, stating that it is better adapted to describing flow and transport channeling effects than a parallel plate model (Neretnieks et al., 1982; Bourke, 1987; Pyrak-Nolte, 1988; Tsang and Tsang, 1989; Tsang et al., 2001) where fracture apertures can be described by normal (Lee et al., 2003), lognormal (e.g., Keller, 1998; Keller et al., 1999), or gamma distributions (Tsang and Tsang, 1987), or a self-affine scale invariance (Plouraboué et al., 1995).

Neuzil and Tracy (1981) presented a model for flow in fractures where the flow is envisioned as occurring in a set of parallel plate openings with different apertures whose distribution was lognormal and used a modified Poiseuille equation.

They showed that the flow conformed to the cubic law and also that the maximum flow occurs through the largest apertures, thereby emphasizing that flow occurs through preferred paths. Thus, in their analysis, the flow depended on the tail of the frequency distribution.

Tsang and Tsang (1987) proposed a theoretical approach to interpret flow in a tightly fractured medium in terms of flow through a system of statistically equivalent one-dimensional channels of variable aperture. The channels were statistically equivalent in the sense that the apertures along each flow channel are generated from the same aperture density distribution and spatial correlation length.

Oron and Berkowitz (1998) have examined the validity of applying the “local cubic law” (LCL) to flow in a fracture which is bounded by impermeable rock surfaces. A two-dimensional order-of-magnitude analysis of the Navier–Stokes equations yields three conditions for the applicability of LCL flow, as a leading-order approximation in a local fracture segment with parallel or nonparallel walls. These conditions demonstrate that the “cubic law” is valid provided that aperture is measured not on a point-by-point basis, but rather as an average over a certain length. Experimental work by Plouraboué et al. (2000) in self-affine rough fractures with various translations of the opposing fracture surfaces indicated that heterogeneity in the flow field caused deviations from the parallel plate model for fracture flow. Some researchers often find it convenient to represent aperture fields in terms of equivalent aperture in the parallel plate model (Zheng et al., 2008).

Brush and Thomson (2003) developed three-dimensional flow models to simulate fluid flow through various random synthetic rough walled fractures created by combining random fields of aperture and the mean wall topography or midsurface, which quantifies undulation about the fracture plane.

The total flow rates from three-dimensional Stokes simulations were within 10 % of LCL simulations with geometric corrections for all synthetic fractures. Differences between the Navier–Stokes (NS) and Stokes simulations clearly demonstrated that inertial forces can significantly influence the internal flow field within a fracture and the total flow rate across a fracture.

Klimczak et al. (2010) carried out flow simulations through fracture networks using the discrete fracture network model (DFN) where flow was modeled through fracture networks with the same spatial distribution of fractures for correlated and uncorrelated fracture length-to-aperture relationships. Results indicate that flow rates are significantly higher for correlated DFNs. Furthermore, the length-to-aperture relations lead to power-law distributions of network hydraulic conductivity which greatly influence equivalent permeability tensor values. These results confirm the importance of the correlated square root relationship of displacement to length scaling for total flow through natural opening-mode fractures and, hence, emphasize the role of these correlations for flow modeling.

Wang et al. (2015) developed and tested a modified LCL (MLCL) taking into account local tortuosity and roughness, and works across a low range of local Reynolds numbers. The MLCL was based on (1) modifying the aperture field by orienting it with the flow direction and (2) correcting for local roughness changes associated with local flow expansion/contraction. In order to test the MLCL, they compared it with direct numerical simulations with the Navier–Stokes equations using real and synthetic three-dimensional rough walled fractures, previously corrected forms of the LCL, and experimental flow tests. The MCL proved to be more accurate than previous modifications of the LCL.

The continuous time random walk (CTRW) approach provides a versatile framework for modeling (non-Fickian) solute transport in fractured media.

Berkowitz et al. (2001) examined a set of analytical solutions based on the CTRW approach to analyze breakthrough data from tracer tests to account for non-Fickian (or scale-dependent) dispersion behavior that cannot be properly quantified by using the advection–dispersion equation.

Cortis and Birkholzer (2008) developed a macroscopic model based on the CTRW framework, to characterize the interaction between the fractured and porous rock domains by using a probability distribution function of residence times. They presented a parametric study of how CTRW parameters evolve, describing transport as a function of the hydraulic conductivity ratio between fractured and porous domains.

Srinivasan et al. (2010) presented a particle-based algorithm that treats a particle trajectory as a subordinated stochastic process that is described by a set of Langevin equations, which represent a CTRW. They used convolution-based particle tracking (CBPT) to increase the computational efficiency and accuracy of these particle-based simulations. The combined CTRW–CBPT approach allows us to convert any particle tracking legacy code into a simulator capable of handling non-Fickian transport.

Dentz et al. (2015) developed a general CTRW approach for transport under radial flow conditions starting from the random walk equations for the quantification of non-local solute transport induced by heterogeneous flow distributions and by mobile–immobile mass transfer processes. They observed power-law tails of the solute breakthrough for broad distributions of particle transit times and particle trapping times. The combined model displayed an intermediate regime, in which the solute breakthrough is dominated by the particle transit times in the mobile zones, and a late time regime that is governed by the distribution of particle trapping times in immobile zones.

The present study is aimed at analyzing the scenario of groundwater contamination of the industrial area of Modugno (Bari – southern Italy) where the limestone aquifer has a fractured and karstic nature.

Previous studies carried out in the same aquifer have applied different conceptual models to model fluid flow and contaminant transport.

Cherubini (2008) applied the discrete feature approach (Diersch, 2002) where the three-dimensional geometry of the subsurface domain describing the matrix structure was combined by interconnected two-dimensional and one-dimensional discrete feature elements in two dimensions in order to simulate, respectively, fractures and karstic cavities in the Bari limestone aquifer. The fracture distribution was inferred from a nonparametric geostatistical analysis (indicator Kriging) of fracture frequency data which had been derived by Rock Quality Designation (RQD) (Priest and Hudson, 1976) data of the contaminated area of the ex Gasometer.

Cherubini et al. (2008) compared the flow modeling results of the previous work with those from a new hydrogeological reconstruction of the heterogeneities in the same aquifer by means of multiple realizations conditioned to borehole data (RQD population), in order to obtain a three-dimensional distribution of fracture frequency, cavities, and terra rossa lenses.

Cherubini and Pastore (2010) applied the nested sequential indicator simulation algorithm to represent the geological architecture of the Bari limestone aquifer which provided reliable prediction of fluid flow. According to phenol transport, the presence of preferential pathways was detected.

Cherubini et al. (2013)a realized a three-dimensional flow model of the Bari limestone aquifer supported by a detailed local-scale geologic model realized by means of sequential indicator simulation (SIS) of lithofacies unit sequences. In this study, a lumped parameter approach was used and calibrated on the groundwater discharge and global hydraulic gradient where fluid flow in fractures was represented by the cubic law, and the Darcy–Weisbach equation was used to estimate the resistance term in the karst network.

Masciopinto et al. (2010) adopted a conceptual model consisting of a three-dimensional parallel set of horizontal planar fractures in between rock layers, each fracture having a variable aperture generated by a stationary random field conditioned to the data derived from pumping-tracer tests. The particle tracking solution was combined with the PHREEQC-2 results to study two-dimensional laminar/non-laminar flow and reactive transport with biodegradation in each fracture of the conceptual model.

Masciopinto and Palmiotta (2013) derived new equations of the fracture aperture as functions of a tortuosity factor to simulate fluid flow and pollutant transport in fractured aquifers. MODFLOW/MT3DMS water velocity predictions were compared with those obtained using a specific software application which solves flow and transport problems in a three-dimensional set of parallel fissures. The results of a pumping/tracer test carried out in a fractured limestone aquifer in Bari (southern Italy) have been used to calibrate advective–dispersive tracer fluxes given by the applied models. Successful simulations of flow and transport in the fractured limestone aquifer were achieved by accommodating the new tortuosity factor in models whose importance lies in the possibility of switching from a discrete to a continuum model by taking into account the effective tracer velocity during flow and transport simulations in fractures.

Masciopinto and Visino (2017) carried out filtration tests on a set of 16 parallel limestone slabs with thicknesses of about 1 cm where rough surfaces and variable fracture apertures had been artificially created. The experimental filtration results suggest that model simulations of perturbed virus transport in fractured soils need to also consider pulse-like sources and sinks of viruses. This behavior cannot be simulated using conventional model equations without including a new kinetic model approach.

The present work focuses on the investigation of the hotspots of aquifer contamination in order to infer the location of the sources.

A rough walled parallel plate model has been implemented and calibrated on the basis of piezometric data and has coupled a geostatistical analysis to infer the values of the equivalent aperture.

The current study introduces a novel approach for simulating flow and transport in fractured aquifers by means of combining a rough walled parallel plate model for the flow simulation coupled with inverse modeling and geostatistical analysis to infer the values of the equivalent aperture together with the random walk theory to reproduce the scenario of contamination.

2 Geological and hydrogeological framework

It is well known that hydraulic properties and consequently fluid circulation and contaminant propagation in carbonate rocks are strongly influenced by the degree of rock fracturing and, in general, the presence of mechanical discontinuities, like faults, joints, or other tectonic elements such as syncline or anticline axes (Caine et al., 1996; Caine and Forster, 1999; Antonellini et al., 2014; Billi et al., 2003). Also, the deformation mechanisms are mainly controlled by the physico-chemical properties of rocks, which are, in turn, the result of different composition, depositional setting, and diagenetic evolution (Zhang and Spiers, 2005; Rustichelli et al., 2012).

From the geological point of view, the investigated area is located in the Murge Plateau corresponding to a broad antiformal structure oriented WNW–ESE and represents the bulging foreland of the Pliocene–Pleistocene Southern Apennines orogenic belt (Pieri et al., 1997; Doglioni, 1994; Forster and Evans, 1991; Korneva et al., 2014; Parise and Pascali, 2003).

The stratigraphy of the Murge area consists of a Variscan crystalline basement topped by 6–7 km thick Mesozoic sedimentary cover (represented by the Calcare di Bari formation) followed by relatively thin and discontinuous Cenozoic and Quaternary deposits (Calcareniti di Gravina formation). Figure 1 shows the simplified geological map of the area of Bari.

2.1 Calcare di Bari formation (Cretaceous)

The Calcare di Bari succession consists of biopeloidal and peloidal wackestones/packstones alternating with stromatolithic bindstones with frequent intercalations of dolomitic limestones and grey dolostones. The formation shows a thickness of about 470 m. Most of the Calcare di Bari formation shows facies features related to peritidal environments; only the upper part suggests a relatively more distal and deeper environment belonging to an external platform setting (Pieri et al., 2012; Fig. 1).

This succession appears stratified and fissured and, where it does not show tectonic discontinuities, it is characterized by a subhorizontal or slightly inclined lying position. This formation is subjected to the complex and relevant karstic phenomena that locally lead to the formation of cavities of different shapes and sizes, partially or completely filled by “terra rossa” deposits. The degree of fracturing degree affecting the Calcare di Bari formation is quite variable and mainly depends on the geological and structural (faults, anticline axis, etc.) evolution of the area, including faulting and folding. Also, the distribution of the local measurement of the RQD index is confirmed by the variability of the electrical resistivity along geoelectrical profiles (with lengths from 500 to 1000 m) and from the propagations of the P and S waves (seismic measurements; length of about 1000 m).

On the basis of borehole and in situ surveys, carried out by private companies, the following was observed.

• The fracturing degree of the Calcare di Bari formation is quite variable and it is expressed by the RQD values that vary between 16 % and 25 % (maximum borehole depths: 30 m). Based on the classification system of Deere and Deere (1988) the rock mass is of “very poor rock quality” (RQD < 25 %).

• Medium values of the Rock Mass Rating (RMR) of about 36, indicate, following the Bieniawski classification (1989), very poor rock mass (class IV).

• In addition, profiles of the electrical resistivity (depth < 30 m) allow us to emphasize the presence of very variable electric resistivity values with variations between 100 (low fracture carbonates rocks) and 1700 Ohm m for very fractured formations, with local values on the order of 3–4000 Ohm m in the case of underground cavities.

• Similarly, the velocities of seismic waves P and S have average values on the order of 1300 and 800 m s⁻¹, respectively, in highly fractured limestones, and 2300 (P) and 1400 (S) m s⁻¹ for compact formations.

2.2 Calcareniti di Gravina formation (lower Pleistocene)

This unit unconformably lies on the Calcare di Bari formation. Its thickness varies from a few meters to 20 m and its depositional environments are related to offshore settings. The lower boundary is transgressive and is locally marked by reddish residual deposits and/or by brackish silty deposits passing upward to shallow water calcarenite rich in bioclasts.

As regards the structural features of these deposits, it is possible to observe that the anticline affecting the Cretaceous succession of the Calcare di Bari formation with an ENE–WSW axial direction (Fig. 1) causes a partial diversion of the water courses, whose path seems to be also influenced by some NE–SW fault (NE of Modugno). The former phenomenon is due to the antithetically dipping flanks of the gentle fold, while the latter effect is likely a consequence of the denser fracturing along the shear zone and hence the increased erodibility of the local outcropping limestone enhancing the water flow concentration.

In general, the limestone bedrock hosts a wide and thick aquifer due to a diffuse rock fracturing and the karstic phenomena.

Moreover, the irregular spatial distribution of the fractures and karstic channels makes the Bari aquifer very anisotropic. The average hydraulic conductivity of this aquifer is generally estimated to be 10⁻³ to 10⁻⁴ m s⁻¹.

The groundwater flows toward the sea, under a low gradient, in different subparallel fractured layers separated by compact (i.e., not fractured) rock blocks.

In proximity to the coast, the carbonate (Mesozoic) stratum contains freshwater flowing in phreatic conditions and floating on underlying saltwater of continental intrusion. The location of the transition zone between freshwater and salt water has thickness and position variable and changes over time depending on the distribution of the hydrostatic pressures of the system.

Figure 1 Simplified geological map of the area of Bari: (1) Calcare di Bari formation (Cretaceous); (2) Calcareniti di Gravina formation (Lower Pleistocene); (3) hydrographic network; (4) anticlinal axis; (5) syncline axis; (6) fault (uncertain); (7): escarpment.

3 Hydrologic and hydrogeologic water budget

The effective infiltration has been estimated by means of the hydrologic and hydrogeologic water budget of the subtended basin. Climatic data registered in the thermopluviometric stations present in the area have been elaborated and the average rainfall module and the monthly evapotranspiration have been calculated for the 3 decades 1974–2005.

Twelve climatic stations have been considered (Bari – hydrographic station, Bari – observatory station, Bitonto, Grumo Appula, Adelfia, Casamassima, Mercadante, Ruvo di Puglia, Corato, Altamura, Santeramo, Gioia del Colle) and for each station the monthly rainfall and evapotranspiration map has been realized by means of the inverse distance weighting algorithm. The latter has been estimated by means of the Thornthwaite method, applying a crop coefficient of 0.40.

The hydrologic and hydrogeologic basins have been defined on the basis of literature data and the regional thematic cartography.

The lithotypes in the study area are principally limestones and calcarenites with secondary permeability, characterized by a high transmissivity. The zones in proximity to tectonic structures create preferential flow paths, but at the same time generate a dismemberment of the aquifer that could not be able to feed the flow downstream. Because of that it proves to be difficult to carry out a zonation of recharge areas, and therefore a constant runoff coefficient of 0.10 has been considered for the whole basin. In Fig. 2 the map of the (a) annual precipitation and (b) estimated evapotranspiration evaluated for the hydrological basin of the study area is shown.

Figure 2 Map of (a) annual precipitation and (b) estimated evapotranspiration evaluated for the hydrological basin of the study area.

4 Well-performance tests: step-drawdown tests

Ninety-eight long-term step-drawdown hydraulic tests have been analyzed in the study area.

A step-drawdown test is a single-well test in which the well is pumped at a low constant discharge rate until the drawdown within the well stabilizes.

Step-drawdown tests can be used to evaluate the characteristics of the well and its immediate environment. Unlike the aquifer test, it is not designed to produce reliable information concerning the aquifer, even though it is possible to estimate the transmissivity of the immediate surroundings of the catchment. This test determines the critical flow rate of the well, as well as the various head losses and drawdowns as functions of pumping rates and times. Finally, it is designed to estimate the well efficiency, to set an exploitation pumping rate, and to specify the depth of installation of the pump.

The total drawdown at a pumping well is given by

$$\begin{array}{}\text{(1)}& s=\left(A_{\mathrm{1}}+A_{\mathrm{2}}\right)\cdot Q+B\cdot Q^{\mathrm{2}},\end{array}$$

where s (L) represents the registered drawdown, Q (L³ T⁻¹) the pumped flow rate, A₁ (T L⁻²) the linear aquifer loss coefficient, and A₂ (T L⁻²), e, and B (T² L⁻⁵) are, respectively, the linear and nonlinear well-loss coefficients.

This equation can be made explicit in terms of aquifer transmissivity T (L² T⁻¹), the transmissivity of the damage zone T_SKIN (L² T⁻¹), and the nonlinear term β (T² L⁻⁴) (Cherubini and Pastore, 2011):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }s=\left[{\displaystyle \frac{\mathrm{1}}{T\mathrm{2}\mathit{\pi }}}\ln \left({\displaystyle \frac{R}{r_{\mathrm{w}}}}\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\left({\displaystyle \frac{\mathrm{1}}{T_{\mathrm{SKIN}}}}-{\displaystyle \frac{\mathrm{1}}{T}}\right)\ln \left({\displaystyle \frac{r_{\mathrm{SKIN}}}{r_{\mathrm{w}}}}\right)\right]Q\\ \text{(2)}& {\displaystyle }& {\displaystyle }+\left[{\displaystyle \frac{\mathit{\beta }}{\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}}}\left({\displaystyle \frac{\mathrm{1}}{r_{\mathrm{w}}}}-{\displaystyle \frac{\mathrm{1}}{R}}\right)\right]Q^{\mathrm{2}},\end{array}$$

where r_w (L) represents the well radius, r_SKIN (L) the radius of the damage zone, and R (L) the radius of influence of the well.

The total drawdown is formed from three components: the hydraulic component of the aquifer assuming a valid Thiem function, a skin function presented by Cooley and Cunningham (1979) assuming that the transmissivity and the radius of the damage zone are, respectively, equal to $T_{\mathrm{SKIN}}=T/\mathrm{2}$ e r_SKIN=2r_w, and a contribution related to nonlinear losses introduced by Wu (2001).

The radius of influence of the well is obtained by means of the Sichart equation:

$$\begin{array}{}\text{(3)}& R=\mathrm{3000}\cdot s\cdot \sqrt{K}.\end{array}$$

In Fig. 3 is reported the statistical distribution of the estimated transmissivity values along the study area.

Figure 3 Statistical distribution of log₁₀ (T).

5 Linear model of regionalization of transmissivity

The geostatistical analysis has been carried out on the log₁₀ transmissivity values using the S-GemS open-source code (Remy, 2004).

The experimental variogram, which provides a description of how the data are related (correlated) with distance, has been calculated (Fig. 4). Because the Kriging algorithm requires a positive definite model of spatial variability, the experimental variogram cannot be used directly. Instead, a model must be fitted to the data to approximately describe the spatial continuity of the data. An exponential model has been used to fit the experimental variogram described by the function:

$$\begin{array}{}\text{(4)}& \mathit{\gamma }\left(h\right)=C\left[\mathrm{1}-\exp \left(-{\displaystyle \frac{h}{a}}\right)\right],\end{array}$$

where C represents the variance (sill), h (L) the lag, and a (L) the correlation length (range). In our case C assumes a value of 1.2 and a a value of 10000 m.

Figure 4 Omnidirectional experimental variogram fitted with an exponential model; sill = 1.2, range = 10000 m.

Figure 5 shows the ordinary Kriging interpolation of log₁₀(T).

Figure 5 Ordinary Kriging interpolation of log₁₀(T).

6 Analysis of piezometric data

Figure 6 shows the spatial distribution of hydraulic heads on the basis of the 2012 measurement campaign. A global trend in the direction of groundwater flow from SW to NE is evident. A relevant aspect is the presence of high hydraulic head values in proximity to ASI and Bosch wells.

Possible explanations for the increase in hydraulic gradient are (1) lower transmissivity as shown through the step drawdown tests; (2) the transition from a more permeable outcropping lithotype to a less permeable one resulting in a decrease in the effective infiltration; (3) the presence of sinkholes and fissures at the surface giving rise to a point source recharge; and (4) hydraulic disconnection due to lower interconnectivity of the fracture system.

The aquifer transmissivity in that zone is on the order of 10⁻⁵ m² s⁻¹. The trend observed in the hydraulic gradient confirms the increase in the aquifer transmissivity from upstream to downstream: in fact, the tests carried out in proximity to the coast have returned a transmissivity value of 10⁻ to 10⁻³ m² s⁻¹.

Figure 6 Measured piezometric heads (m, slm) from the February 2012 monitoring campaign.

7 Analysis of the scenario of contamination for the study area

The various monitoring campaigns carried out have shown a contamination by chlorinated aliphatic hydrocarbons which, unlike petroleum products, are denser than water and can exist as dense non-aqueous phase liquids (DNAPLs).

The presence of two hotspot areas has been detected, located upstream of the groundwater flow, coherently with the state of contamination detected downstream.

Figure 7 shows the location of the detected contamination (µg L⁻¹).

The pollution indicator has been chosen on the basis of the toxicologic and cancirogenic parameters, the solubility, the sorption coefficient, and the maximum detected contaminant concentration. On the basis of the results of this screening, the tetrachloroethylene (PCE) has the highest concentration, as well as low values of reference dose factors (RDFs) and slope factors (SFs).

Figure 7 Location of the detected contamination by PCE (µg L⁻¹).

8 Parallel rough walled fracture model

The simplest model of flow through rock fractures is the parallel plate model (Huitt, 1956; Snow, 1965) which conceptualizes the fractured medium as made by a set of smooth parallel plates with the same hydraulic aperture b (L) that are separated by a uniform distance. This is actually the only geometrical fracture model for which an exact calculation of the hydraulic conductivity is possible.

Natural fractures present rough walls and complex geometries. Nonlinear flow may occur through rough walled rock fractures; as a consequence, the inertial effect dominates the flow dynamics, giving rise to a deviation from Darcy's law. Fluid flow through a set of natural fracture planes can be expressed using the Darcy–Weisbach equation:

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{d}h}{\mathrm{d}x}}=-{\displaystyle \frac{f}{D}}{\displaystyle \frac{v^{\mathrm{2}}}{\mathrm{2}g}},\end{array}$$

where D (L) represents the hydraulic diameter (2b for the parallel plate model), f is the Darcy–Weisbach coefficient, h (L) is the hydraulic head, x (L) is the distance, and v (L T⁻¹) is the average velocity in fracture calculated as

$$\begin{array}{}\text{(6)}& v={\displaystyle \frac{q}{n_{\mathrm{f}}b}},\end{array}$$

where q(L² T⁻¹) is the volumetric flow rate per unit length of fracture and n_f (–) is the number of fractures.

The Darcy–Weisbach equation can be rewritten in terms of volumetric flow per unit length:

$$\begin{array}{}\text{(7)}& \stackrel{\mathrm{‾}}{q}=-\left[n_{\mathrm{f}}b{\displaystyle \frac{\sqrt{\frac{\mathrm{4}b}{f}g}}{\sqrt{\left|\mathrm{\nabla }h\right|}}}\right]\mathrm{\nabla }h.\end{array}$$

The term in square brackets represents the equivalent hydraulic transmissivity T_eq(f,∇h) of the n_f rough walled fractures.

Figure 8 Map of log10 of aquifer transmissivity determined by means of the inverse modeling algorithm.

Figure 9 Map of simulated hydraulic heads.

Figure 10 Graph of the calibration.

The Darcy–Weisbach coefficient or friction factor depends on the flow regime. In the case of a smooth walled fracture and linear flow regime, f is equal to

$$\begin{array}{}\text{(8)}& f={\displaystyle \frac{\mathrm{96}}{\mathit{Re}}},\end{array}$$

where Re represents the Reynolds number:

$$\begin{array}{}\text{(9)}& \mathit{Re}={\displaystyle \frac{\mathit{\rho }vD}{\mathit{\mu }}}.\end{array}$$

Substituting Eq. (8) into Eq. (7), the cubic law (Witherspoon et al., 1980) where q is proportional to the cubic power of the fracture aperture is obtained:

$$\begin{array}{}\text{(10)}& q=n_{\mathrm{f}}{\displaystyle \frac{\mathit{\rho }g}{\mathit{\mu }}}{\displaystyle \frac{b^{\mathrm{3}}}{\mathrm{12}}}.\end{array}$$

The cubic law is not always adequate to represent the flow process in natural fractures, so a deviation from linearity can be observed.

The friction factor depends on the flow regime described by the Reynolds number and can be represented by the following relationship found by Nazridoust at al. (2006):

$$\begin{array}{}\text{(11)}& f={\displaystyle \frac{\mathrm{123}}{\mathit{Re}}}\left(\mathrm{1}+\mathrm{0.12}{\mathit{Re}}^{\mathrm{0.687}}\right).\end{array}$$

9 Inverse flow modeling

Inverse modeling is a technique used to estimate unknown model parameters using as input data punctual values of the state variables (hydraulic head flow). Generally, in real problems the number of parameters to estimate (n) is higher than the number of measured values (m). For example, this is the case with mapping hydraulic transmissivity values varying continuously in space.

For underdetermined inverse problems of this kind the objective function (L) can be written in this way:

$$\begin{array}{}\text{(12)}& L\left(\mathit{y}\mathbf{,}\mathit{s}\right)=L_{\mathrm{fitness}}\left(\mathit{y}\mathbf{,}\mathit{s}\right)+L_{\mathrm{penalty}}\left(\mathit{s}\right),\end{array}$$

where s represents the vector of measured values of state variables (e.g. hydraulic transmissivity), and y represents the vector of parameter values.

The fitness function responds to maximum likelihood criteria between the observed and simulated values and can be written as

$$\begin{array}{}\text{(13)}& L_{\mathrm{fitness}}\left(\mathit{y}\mathbf{,}\mathit{s}\right)={\left(\mathit{y}\mathbf{-}\mathit{h}\mathbf{\left(}\mathit{s}\mathbf{\right)}\right)}^T{\mathbf{R}}^{-\mathrm{1}}\left(\mathit{y}\mathbf{-}\mathit{h}\mathbf{\left(}\mathit{s}\mathbf{\right)}\right),\end{array}$$

where h represents the model that, starting from the parameter vector, estimates the state variable, and R is the measurement error covariance matrix. Generally this function can be reduced to the square root of the sum of the squared difference between the measured and simulated RMSE:

$$\begin{array}{}\text{(14)}& L_{\mathrm{fitness}}={\displaystyle \frac{{∥\mathit{y}\mathbf{-}\mathit{h}\left(\mathit{s}\right)∥}^{\mathrm{2}}}{\mathrm{\Delta }{H}^{\mathrm{2}}}},\end{array}$$

where ΔH represents a parameter of accuracy of observed data.

The penalty function is used to discriminate the solutions with values of the fitness function comparable by means of geostatistical criteria (Kitanidis, 1995):

$$\begin{array}{}\text{(15)}& L_{\mathrm{penalty}}={\left(\mathit{s}\mathbf{-}\mathit{X}\mathit{\beta }\right)}^T{\mathbf{Q}}^{-\mathrm{1}}\left(\mathit{s}\mathbf{-}\mathit{X}\mathit{\beta }\right).\end{array}$$

where Q represents the spatial covariance matrix, X is a unit vector, and β is the mean of the values of the parameters. The penalty function can be rewritten by eliminating β:

$$\begin{array}{}\text{(16)}& L_{\mathrm{penalty}}={\mathit{s}}^{\mathit{T}}\mathit{G}\mathit{s}\mathit{G}\equiv {\mathbf{Q}}^{-\mathbf{1}}-{\displaystyle \frac{{\mathbf{Q}}^{-\mathbf{1}}\mathit{X}{\mathit{X}}^{\mathit{T}}{\mathbf{Q}}^{-\mathrm{1}}}{{\mathit{X}}^T{\mathbf{Q}}^{-\mathrm{1}}\mathit{X}}}.\end{array}$$

The common assumption is that the spatial distribution of the parameters follows the geostatistical distribution defined by the variogram. Under this hypothesis the covariance matrix present in the penalty function can be defined as

$$\begin{array}{}\text{(17)}& Q_{ij}=\mathrm{2}\mathit{\gamma }\left(\left|x_i-x_j\right|\right)i,j=\mathrm{1},\mathrm{\dots },n.\end{array}$$

10 Solute transport modeling

Solute transport in fracture neglecting the effect of matrix diffusion and the chemical reactions can be described by the following advection–dispersion equation:

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{\partial c}{\partial t}}+\stackrel{\mathrm{‾}}{v}\cdot \mathrm{\nabla }c=\mathrm{\nabla }\cdot \left(\mathbf{D}\mathrm{\nabla }c\right),\end{array}$$

where c (M L⁻³) is the concentration of a solute and D (L² T⁻¹) is the symmetric dispersion tensor with the following components:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }Dxx=\left({\mathit{\alpha }}_{\mathrm{L}}{v}_x^{\mathrm{2}}+{\mathit{\alpha }}_{\mathrm{T}}{v}_y^{\mathrm{2}}\right)/\left|,v\right|\\ {\displaystyle }& {\displaystyle }Dyy=\left({\mathit{\alpha }}_{\mathrm{T}}{v}_x^{\mathrm{2}}+{\mathit{\alpha }}_{\mathrm{L}}{v}_y^{\mathrm{2}}\right)/\left|v\right|,\\ \text{(19)}& {\displaystyle }& {\displaystyle }Dxy=({\mathit{\alpha }}_{\mathrm{L}}-{\mathit{\alpha }}_{\mathrm{T}})v_x{v}_y/\left|v\right|,\end{array}$$

where α_L (L) and α_T (L) are the longitudinal and transverse dispersivities, respectively.

For pure advective transport the particle moves along the flow lines. In order to represent dispersion phenomena, the random walk method adds a random displacement to each particle, independently of the other particles, in addition to advective displacement.

For a given time step Δt, considering the tensorial nature of the dispersion and the spatially variable velocity field, each particle moves according to

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{x}_p\left(t+\mathrm{\Delta }t\right)=x_p\left(t\right)+v_x^{\prime }\mathrm{\Delta }t+Z_{\mathrm{1}}\sqrt{\mathrm{2}{D}_L\mathrm{\Delta }t}{\displaystyle \frac{v_x}{\left|v\right|}}\\ {\displaystyle }& {\displaystyle }-Z_{\mathrm{2}}\sqrt{\mathrm{2}{D}_T\mathrm{\Delta }t}{\displaystyle \frac{v_y}{\left|v\right|}},\\ {\displaystyle }& {\displaystyle }{y}_p\left(t+\mathrm{\Delta }t\right)=y_p\left(t\right)+v_y^{\prime }\mathrm{\Delta }t+Z_{\mathrm{1}}\sqrt{\mathrm{2}{D}_L\mathrm{\Delta }t}{\displaystyle \frac{v_y}{\left|v\right|}}\\ \text{(20)}& {\displaystyle }& {\displaystyle }+Z_{\mathrm{2}}\sqrt{\mathrm{2}{D}_T\mathrm{\Delta }t}{\displaystyle \frac{v_x}{\left|v\right|}},\end{array}$$

with

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{v}_x^{\prime }=v_x+{\displaystyle \frac{\partial D_{xx}}{\partial x}}+{\displaystyle \frac{\partial D_{xy}}{\partial y}},\\ {\displaystyle }& {\displaystyle }{v}_y^{\prime }=v_y+{\displaystyle \frac{\partial D_{xy}}{\partial x}}+{\displaystyle \frac{\partial D_{yy}}{\partial y}},\\ {\displaystyle }& {\displaystyle }{D}_L={\mathit{\alpha }}_L\left|v\right|,\\ \text{(21)}& {\displaystyle }& {\displaystyle }{D}_T={\mathit{\alpha }}_{\mathrm{T}}\left|v\right|,\end{array}$$

where Z₁ and Z₂ are two normally distributed random variables. For steady-state flow and for a source constant intensity, the assumption that the particles N released in time interval (t₁, t₁+Δt) follow exactly the same random trajectories of the particles N released during the previous interval (t₁, t₁−Δt) is possible. Under this assumption only N particles are needed to simulate the location of the particles at a previous time step.

11 Results and discussion

11.1 Flow modeling

The MODFLOW numerical code coupled with the inverse model approach presented in the previous section has been used to model groundwater flow.

The numerical simulations have been carried out on a two-dimensional domain of 968.7 km². The domain has been discretized by means of a structured grid of 100 m in size.

In correspondence to the coastline, a first type of boundary condition has been imposed (h=0 m), and along the detected watershed a second type of boundary condition (q=0 m² s⁻¹); the recharge from upstream is simulated by means of a first type of boundary condition where the hydraulic heads are equal to the detected regional values h=32–41 m (Piano di Tutela delle Acque Regione Puglia, Tav. 6.2 http://old.regione.puglia.it/index.php?page=documenti&id=29&opz=getdoc, last access: February 2018).

A second type of boundary condition on the whole simulation domain has been imposed that concerns the mean effective infiltration calculated from the hydrologic budget q=0.037 m³ d⁻¹.

The algorithm of inverse modeling has been applied to carry out the estimation of the spatial distribution of the equivalent transmissivity (Fig. 8) on the basis of the observed hydraulic head (vector y), the regionalization model (matrix Q) described by the variogram of the logarithm of the hydraulic transmissivity determined in the previous section.

The inverse model algorithm follows those steps. (1) Start from a conditional simulation of the log of T_eq determined by means of the hydraulic tests conducted in the area. (2) A set of pilot points are chosen in the area using a regularly spaced criterion and the value of T_eq has been determined for each pilot point (vector s). (3) By means of the ordinary Kriging interpolation of the pilot points the map of T_eq is obtained and represents the input datum of the flow numerical model. (4) The hydraulic head has been determined using the flow numerical model (vector h) and the values of the objective function have been determined using Eq. (12). (5) The values of T_eq are updated for each pilot point.

Using the Levenberg–Marquardt algorithm, the values of T_eq for each pilot point are updated as long as the objective function is minimized.

Figure 9 shows the results obtained from the flow model in steady-state condition, calibrated with the measurement campaign of February 2012 (Table 1).

Table 2 shows the data of model calibration and Fig. 10 shows the graph of the calibration.

The outcomes of the calibration are satisfactory. The comparison between the simulated and observed data has given a mean absolute residual equal to 0.57 m, an RMSE equal to 4.57 m, and a correlation coefficient r² equal to 0.997. In the following figures and tables are shown the results for the flow model.

Table 1 Comparison between the observed and simulated hydraulic heads with related residuals, relative to the measurement campaign of February 2012.

Table 2 Data of model calibration.

The simulated hydraulic head distribution together with the equivalent transmissivity map put into evidence how the anticline affecting the Calcare di Bari formation in direction ENE–WSW influences the flow directions. Furthermore, they highlight how the hydraulic circulation is more active along the coast coherently with a higher degree of fracturing and karst phenomena.

Once the equivalent hydraulic transmissivity map has been obtained, and assuming a value of the number of sets of fractures n_f, the spatial distribution of the mean equivalent aperture and the velocity field can be obtained.

Assuming the cubic law to be valid, the mean equivalent aperture can be obtained as

$$\begin{array}{}\text{(22)}& b_{\mathrm{eq}}=\sqrt[\mathrm{3}]{\mathrm{12}{\displaystyle \frac{T_{\mathrm{eq}}}{n_{\mathrm{f}}}}{\displaystyle \frac{\mathit{\mu }}{\mathit{\rho }g}}}.\end{array}$$

The velocity field results:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{v}_x=-{\displaystyle \frac{\mathit{\rho }g}{\mathit{\mu }}}{\displaystyle \frac{b_{\mathrm{eq}}^{\mathrm{2}}}{\mathrm{12}}}{h}_x,\\ \text{(23)}& {\displaystyle }& {\displaystyle }{v}_y=-{\displaystyle \frac{\mathit{\rho }g}{\mathit{\mu }}}{\displaystyle \frac{b_{\mathrm{eq}}^{\mathrm{2}}}{\mathrm{12}}}{h}_y,\end{array}$$

whereas assuming the Darcy–Weisbach equation to be valid, the mean equivalent aperture and the flow field can be obtained by means of the following iterative steps starting from the values of b_eq, v_x, and v_y previously evaluated:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{Re}}^k={\displaystyle \frac{\mathit{\rho }\left|v^k\right|\mathrm{2}{b}_{\mathrm{eq}}^k}{\mathit{\mu }}},\\ {\displaystyle }& {\displaystyle }{f}^{k+\mathrm{1}}={\displaystyle \frac{\mathrm{123}}{{\mathit{Re}}^k}}\left(\mathrm{1}+\mathrm{0.12}\left({\mathit{Re}}^{k\mathrm{0.687}}\right)\right),\\ {\displaystyle }& {\displaystyle }{b}_{\mathrm{eq}}^{k+\mathrm{1}}=\sqrt[\mathrm{3}]{{\displaystyle \frac{T_{\mathrm{eq}}}{n_{\mathrm{f}}^{\mathrm{2}}}}{\displaystyle \frac{f^{k+\mathrm{1}}}{\mathrm{4}g}}\left|\mathrm{\nabla }h\right|},\\ {\displaystyle }& {\displaystyle }{v}_x^{k+\mathrm{1}}=-{\displaystyle \frac{\sqrt{\frac{\mathrm{4}{b}_{\mathrm{eq}}^{k+\mathrm{1}}}{f^{k+\mathrm{1}}}g}}{\sqrt{\left|\mathrm{\nabla }h\right|}}}{h}_x,\\ \text{(24)}& {\displaystyle }& {\displaystyle }{v}_y^{k+\mathrm{1}}=-{\displaystyle \frac{\sqrt{\frac{\mathrm{4}{b}_{\mathrm{eq}}^{k+\mathrm{1}}}{f^{k+\mathrm{1}}}g}}{\sqrt{\left|\mathrm{\nabla }h\right|}}}{h}_y.\end{array}$$

Figure 11 shows the relative percentage of error in the flow velocity magnitude for different numbers of fractures.

Figure 11 Relative percentage of error in the flow velocity magnitude for different numbers of fractures: (a) n_f= 4; (b) n_f= 28.

Figure 12 Steady-state distribution of hypothetical contamination using the random walk model with the source contamination localized in correspondence to the hotspot of the contamination considering a number of fracture of n_f= 20 and a longitudinal and transversal dispersion coefficient equal to α_L= 70 m and α_T= 7 m.

Figure 13 Breakthrough curves of hypothetical continuous contamination released in correspondence to the hotspot, determined for linear and nonlinear flow models, evaluated at the downstream boundary.

Figure 14 Mean travel time t_m at varying numbers of fractures n_f for linear and nonlinear models.

As the number of fractures increases, the velocity magnitude decreases; therefore, the friction factor reaches a value of 96∕Re. Anyway, the percentage of error in the flow velocity magnitude seems not to be negligible. In fact, for n_f=28 a minimum value of 8 % is obtained, reaching a value of 28 % for n_f=4. These results show that under natural hydraulic gradient conditions in fractured limestone the nonlinearity of the flow cannot be negligible. It is clear that under a forced hydraulic gradient due to anthropic stresses the equivalent transmissivity decreases dramatically, with a value less than 40 % of Darcian-like hydraulic transmissivity (Cherubini et al., 2012).

11.2 Analysis of the scenario of contamination

A particle tracking transport method has been applied for the simulation of contaminant transport. A punctual source contamination has been imposed in correspondence to the detected hotspot equal to a PCE concentration of 1283 µg L⁻¹. This localization is coherent with the soil contamination detected in the area.

In order to solve the advective transport equation a numerical Lagrangian particle-based random walk method is implemented. For each time step a constant number of 500 particles have been released into the domain from the source. According to the Rauch et al. (2005) assumption reported in the previous section, only 500 particles are needed to simulate the steady-state distribution of the hypothetical contamination. Even if the source of contamination has been considered punctual, the obtained simulation scenario proves to be compatible with the observed one, and therefore it is possible to assume that the sources of contamination are located in correspondence to the detected hotspot (Fig. 12).

Figure 13 shows the breakthrough curves of hypothetical continuous contamination released in correspondence to the hotspot, determined for linear and nonlinear flow models, evaluated at the downstream boundary for n_f = 20.

Figure 14 shows the mean travel time at a varying number of fractures for the linear and nonlinear models. With an increasing number of fractures, the travel time increases in a linear way, because the cross-sectional area increases as well. The figures highlight that travel time for the nonlinear model is higher than the linear assumption. In a particular way, the percentages of error are in the range of 6.22 %–5.34 % passing from n_f= 4 (Re = 0.02–10.60) to n_f=28 (Re  = 0.002–1.51). This is coherent with what was detected by Cherubini et al. (2012, 2013b, 2014), who carried out hydraulic and tracer tests on an artificially created fractured rock sample and found a pronounced mobile–immobile zone interaction leading to a non-equilibrium behavior of solute transport.

The existence of a non-Darcian flow regime was shown to influence the velocity field by giving rise to a delay in solute migration with respect to the values that could be obtained under the assumption of a linear flow field. Furthermore, the presence of inertial effects was shown to enhance non-equilibrium behavior. In a particular manner, they found that the percentage of error in the travel time with respect to the linear flow assumption varied in the range of 5.90 %–40.75 %, corresponding to a range of Re of 29.48–52.16. These results highlight the fact that as the scale of observation increases, the error in the mean travel time with respect to the linear flow model becomes more relevant. In fact, at field scale also for Re, just above the unit (n_f = 28), the error is equal to 5.34 %, comparable with the error of 5.90 % found at laboratory scale for Re equal to 29.48. This means that under anthropic stresses multiple pumpings or injections give rise to a higher flow velocity and then higher Re, leading to a dramatic delay on contaminant transport. Therefore, nonlinear flow must be considered in order to have a more accurate estimation of the breakthrough curve and mean travel time of contaminated scenarios.

In fracture networks, the presence of nonlinear flow plays an important role in the distribution of the solutes according to the different pathways. In fact, the energy spent to cross the path should be proportional to the resistance to flow associated with the single pathway, which in a nonlinear flow regime is not constant but depends on the flow rate. This means that by changing the boundary conditions, the resistance to flow varies and consequently the distribution of solute in the main and secondary pathways also changes, giving rise to a different behavior of solute transport (Cherubini et al., 2014).

This concept has to be taken into account in case of cleanup of the aquifer using for example the Pump & Treat system. The multiple pumping and reinjection of the treated groundwater give rise to a higher flow velocity in the aquifer, resulting in a much greater hydraulic gradient. In this case nonlinear flow behavior has to be taken into account in order to obtain more accurate cleanup strategies.

12 Conclusions

The present study is aimed at analyzing the scenario of groundwater contamination (by investigating the hotspots) of the industrial area of Modugno (Bari – southern Italy) where the limestone aquifer has a fractured and karstic nature.

The presence of hotspot areas has been detected, located upstream of the groundwater flow, coherently with the state of contamination detected downstream and the soil contamination.

A rough walled parallel plate model has been implemented and calibrated on the basis of piezometric data and has coupled a geostatistical analysis to infer the values of the equivalent aperture. Using the random walk theory, the steady-state distribution of hypothetical contamination with the source contamination at the hotspot has been carried out.

The flow and transport models have reproduced the flow pattern well and have given a pollution scenario that is compatible with the observed one.

From an analysis of the flow and transport pattern it is possible to infer that the anticline affecting the Calcare di Bari formation in direction ENE–WSW influences the direction of flow as well as the propagation of the contaminant.

The results also show that the presence of nonlinear flow influences advection, in that it leads to a delay in solute transport with respect to the linear flow assumption. Moreover, the distribution of solute according to different pathways is not constant, but is related to the flow rate.

This is due to the non-proportionality between the energy spent to cross the path and the resistance to flow for fractured media, which affects the distribution of the solutes according to the different pathways.

The obtained results represent the fundamental basis for a detailed study of the contaminant propagation in correspondence to the hotspot area in order to find the best cleanup strategies and optimize any anthropic intervention on the industrial site.

Future developments of the current study will be to implement a transient model and to include the density-dependent flow in the simulations.