**Modeling the Long-Term Fate of Agricultural Nitrate in Groundwater in the San Joaquin Valley, California**

Francis H. Chapelle, Bruce G. Campbell, Mark A. Widdowson and Mathew K. Landon

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/53652

## **1. Introduction**

Nitrate contamination of groundwater systems used for human water supplies is a major environmental problem in many parts of the world. Fertilizers containing a variety of re‐ duced nitrogen compounds are commonly added to soils to increase agricultural yields. But the amount of nitrogen added during fertilization typically exceeds the amount of nitrogen taken up by crops. Oxidation of reduced nitrogen compounds present in residual fertilizers can produce substantial amounts of nitrate which can be transported to the underlying wa‐ ter table. Because nitrate concentrations exceeding 10 mg/L in drinking water can have a va‐ riety of deleterious effects for humans, agriculturally derived nitrate contamination of groundwater can be a serious public health issue.

The Central Valley aquifer of California accounts for 13 percent of all the groundwater with‐ drawals in the United States [1]. The Central Valley, which includes the San Joaquin Valley, is one of the most productive agricultural areas in the world and much of this groundwater is used for crop irrigation. However, rapid urbanization has led to increasing groundwater withdrawals for municipal public water supplies. That, in turn, has led to concern about how contaminants associated with agricultural practices will affect the chemical quality of groundwater in the San Joaquin Valley [2]. Crop fertilization with various forms of nitrogencontaining compounds can greatly increase agricultural yields. However, leaching of nitrate from soils due to irrigation has led to substantial nitrate contamination of shallow ground‐ water [3]. That shallow nitrate-contaminated groundwater has been moving deeper into the

© 2013 Chapelle et al.; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 Chapelle et al.; licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Central Valley aquifer since the 1960s [3]. Denitrification can be an important process limit‐ ing the mobility of nitrate in groundwater systems [4]. However, substantial denitrification requires adequate sources of electron donors in order to drive the process. In many cases, dissolved organic carbon (DOC) and particulate organic carbon (POC) are the primary elec‐ tron donors driving active denitrification in groundwater. The purpose of this chapter is to use a numerical mass balance modeling approach to quantitatively compare sources of elec‐ tron donors (DOC, POC) and electron acceptors (dissolved oxygen, nitrate, and ferric iron) in order to assess the potential for denitrification to attenuate nitrate migration in the Cen‐ tral Valley aquifer.

**2. Methods**

**2.1. Study area**

The San Joaquin Valley occupies the southern two-thirds of the Central Valley of California (Figure 1), a large northwest-trending, asymmetric structural trough filled with marine and continental sediments up to 10 km thick [18]. East of the valley the Sierra Nevada rise to an elevation of more than 4,200 m. West of the valley are the Coast Ranges, a series of parallel ridges of moderate elevation. Streams in the northern part of the San Joaquin Valley drain through the San Joaquin River northward to the San Francisco Bay. During predevelopment, groundwater generally moved toward the center of the valley where it discharged to the San Joaquin River. However, extensive development of groundwater for agriculture and public

Modeling the Long-Term Fate of Agricultural Nitrate in Groundwater in the San Joaquin Valley, California

http://dx.doi.org/10.5772/53652

153

The hydrologic system in the Modesto area (Figure 1) is complex, in part because of the het‐ erogeneous nature of the hydrogeological setting. The primary aquifers in the study area are a complex sequence of overlapping alluvial fan deposits that have been eroded from the Si‐ erra Nevada (Figure 2). These alluvial fan deposits consist of coarse-grained sands and grav‐ els with discontinuous clayey silts and clays [19]. A relevant characteristic of these sediments for the present study is that they contain very low amounts of organic carbon, typically in the range of 0.01 to 0.1 weight percent organic carbon (Figure 2). The low sedi‐ ment organic carbon content reflects the generally arid climate in the recent geologic past,

The relatively low amounts of available organic carbon have resulted in groundwater, that is predominantly oxic [20] (Figure 1). Most of the individual analyses (Figure 1) of ground‐ water are either oxic or mixed oxic and anoxic. The specific criteria used to distinguish re‐ dox conditions are described in [20]. In general, concentrations of dissolved oxygen tend to decrease and concentrations of dissolved iron tend to increase as groundwater approaches

The numerical model in this study is based on a regional model of groundwater conditions in the Central-Eastern San Joaquin Valley [19, 21]. The U.S. Geological Survey (USGS) threedimensional finite-difference code MODFLOW-2000 [22] was used to simulate groundwater flow and water-level distributions across the study area (Fig. 1). This regional model was constructed using a three-dimensional grid consisting of 153 rows and 137 columns and 16 layers. The 400 m by 400 m cells were uniform in dimension, and model layers varied from 0.5 m to 16 m above layer 8 and from 20 to 74 m below layer 8. The model area extends from the Coast Ranges to the Sierra Nevada foothills, although the area west of the San Joaquin River was not simulated. The external boundaries of the regional model were a no-flow boundary on the northeastern boundary and general head boundaries on the other three sides. Hydraulic conductivities were estimated based on sediment texture documented in drilling logs. Complete details of the model, including input parameters, calibration, flow

and the overland transport of alluvial fan deposits prior to final deposition.

water supply has substantially altered the flow system.

the discharge area near the San Joaquin River.

budget and travel times are described in [21].

**2.2. MODFLOW model of the San Joaquin Valley aquifer**

#### **1.1. Interactions of dissolved organic carbon, oxygen, and nitrate**

There are at least three distinct compartments present in groundwater systems that store natural organic carbon capable of reacting with dissolved oxygen, nitrate, and other electron acceptors. Dissolved organic carbon (DOC) is present at varying concentrations in all groundwater systems [5,6,7], and this dissolved compartment can store substantial amounts of organic carbon. In addition to DOC, many groundwater systems contain particulate or‐ ganic carbon (POC) in various stages of diagenesis [5,6,8,9,10]. Microbial degradation of POC can be an additional source of DOC to soil interstitial water [10] and groundwater [11]. Finally, silicate, iron oxyhydroxide, and other minerals present in aquifer solids have a ca‐ pacity to adsorb DOC [12], removing it from the dissolved phase [13,14,9,15]. These adsorp‐ tion processes are partially reversible, so that desorption of organic carbon from aquifer materials is also a potential source of DOC [16]. A mass balance model of organic carbon dynamics in groundwater systems, therefore, will need to account for each of these carbonstoring compartments and their interactions.

In contrast to the complexity inherent in the multiple sources, sinks, and composition of DOC, atmospheric oxygen carried through the unsaturated zone by infiltrating precipitation is the sole source of dissolved oxygen (DO) to groundwater systems which lack active pho‐ tosynthesis. In addition, DO's relatively limited solubility in fresh water (10.1 mg/L at 15 °C) provides a convenient upper limit to concentrations of DO that can be delivered to the water table. These characteristics will be useful in constraining a quantitative mass balance be‐ tween DOC and DO.

The interaction of DO with the three compartments of organic carbon present in groundwa‐ ter systems determines the transformation or lack of transformation of nitrate. The usual ecological succession of electron acceptor utilization in groundwater systems (oxygen>ni‐ trate>Fe(III)> sulfate> carbon dioxide [17] implies that once concentrations of DO drop be‐ low approximately 0.5 mg/L, reduction of nitrate will occur and may coincide with any of the succeeding predominant terminal electron-accepting processes. Constructing a mass bal‐ ance between the sources and bioavailability of DOC, DO, and nitrate, therefore, is central to assessing the fate and transport of nitrate.

## **2. Methods**

Central Valley aquifer since the 1960s [3]. Denitrification can be an important process limit‐ ing the mobility of nitrate in groundwater systems [4]. However, substantial denitrification requires adequate sources of electron donors in order to drive the process. In many cases, dissolved organic carbon (DOC) and particulate organic carbon (POC) are the primary elec‐ tron donors driving active denitrification in groundwater. The purpose of this chapter is to use a numerical mass balance modeling approach to quantitatively compare sources of elec‐ tron donors (DOC, POC) and electron acceptors (dissolved oxygen, nitrate, and ferric iron) in order to assess the potential for denitrification to attenuate nitrate migration in the Cen‐

There are at least three distinct compartments present in groundwater systems that store natural organic carbon capable of reacting with dissolved oxygen, nitrate, and other electron acceptors. Dissolved organic carbon (DOC) is present at varying concentrations in all groundwater systems [5,6,7], and this dissolved compartment can store substantial amounts of organic carbon. In addition to DOC, many groundwater systems contain particulate or‐ ganic carbon (POC) in various stages of diagenesis [5,6,8,9,10]. Microbial degradation of POC can be an additional source of DOC to soil interstitial water [10] and groundwater [11]. Finally, silicate, iron oxyhydroxide, and other minerals present in aquifer solids have a ca‐ pacity to adsorb DOC [12], removing it from the dissolved phase [13,14,9,15]. These adsorp‐ tion processes are partially reversible, so that desorption of organic carbon from aquifer materials is also a potential source of DOC [16]. A mass balance model of organic carbon dynamics in groundwater systems, therefore, will need to account for each of these carbon-

In contrast to the complexity inherent in the multiple sources, sinks, and composition of DOC, atmospheric oxygen carried through the unsaturated zone by infiltrating precipitation is the sole source of dissolved oxygen (DO) to groundwater systems which lack active pho‐ tosynthesis. In addition, DO's relatively limited solubility in fresh water (10.1 mg/L at 15 °C) provides a convenient upper limit to concentrations of DO that can be delivered to the water table. These characteristics will be useful in constraining a quantitative mass balance be‐

The interaction of DO with the three compartments of organic carbon present in groundwa‐ ter systems determines the transformation or lack of transformation of nitrate. The usual ecological succession of electron acceptor utilization in groundwater systems (oxygen>ni‐ trate>Fe(III)> sulfate> carbon dioxide [17] implies that once concentrations of DO drop be‐ low approximately 0.5 mg/L, reduction of nitrate will occur and may coincide with any of the succeeding predominant terminal electron-accepting processes. Constructing a mass bal‐ ance between the sources and bioavailability of DOC, DO, and nitrate, therefore, is central to

**1.1. Interactions of dissolved organic carbon, oxygen, and nitrate**

152 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

storing compartments and their interactions.

assessing the fate and transport of nitrate.

tral Valley aquifer.

tween DOC and DO.

## **2.1. Study area**

The San Joaquin Valley occupies the southern two-thirds of the Central Valley of California (Figure 1), a large northwest-trending, asymmetric structural trough filled with marine and continental sediments up to 10 km thick [18]. East of the valley the Sierra Nevada rise to an elevation of more than 4,200 m. West of the valley are the Coast Ranges, a series of parallel ridges of moderate elevation. Streams in the northern part of the San Joaquin Valley drain through the San Joaquin River northward to the San Francisco Bay. During predevelopment, groundwater generally moved toward the center of the valley where it discharged to the San Joaquin River. However, extensive development of groundwater for agriculture and public water supply has substantially altered the flow system.

The hydrologic system in the Modesto area (Figure 1) is complex, in part because of the het‐ erogeneous nature of the hydrogeological setting. The primary aquifers in the study area are a complex sequence of overlapping alluvial fan deposits that have been eroded from the Si‐ erra Nevada (Figure 2). These alluvial fan deposits consist of coarse-grained sands and grav‐ els with discontinuous clayey silts and clays [19]. A relevant characteristic of these sediments for the present study is that they contain very low amounts of organic carbon, typically in the range of 0.01 to 0.1 weight percent organic carbon (Figure 2). The low sedi‐ ment organic carbon content reflects the generally arid climate in the recent geologic past, and the overland transport of alluvial fan deposits prior to final deposition.

The relatively low amounts of available organic carbon have resulted in groundwater, that is predominantly oxic [20] (Figure 1). Most of the individual analyses (Figure 1) of ground‐ water are either oxic or mixed oxic and anoxic. The specific criteria used to distinguish re‐ dox conditions are described in [20]. In general, concentrations of dissolved oxygen tend to decrease and concentrations of dissolved iron tend to increase as groundwater approaches the discharge area near the San Joaquin River.

## **2.2. MODFLOW model of the San Joaquin Valley aquifer**

The numerical model in this study is based on a regional model of groundwater conditions in the Central-Eastern San Joaquin Valley [19, 21]. The U.S. Geological Survey (USGS) threedimensional finite-difference code MODFLOW-2000 [22] was used to simulate groundwater flow and water-level distributions across the study area (Fig. 1). This regional model was constructed using a three-dimensional grid consisting of 153 rows and 137 columns and 16 layers. The 400 m by 400 m cells were uniform in dimension, and model layers varied from 0.5 m to 16 m above layer 8 and from 20 to 74 m below layer 8. The model area extends from the Coast Ranges to the Sierra Nevada foothills, although the area west of the San Joaquin River was not simulated. The external boundaries of the regional model were a no-flow boundary on the northeastern boundary and general head boundaries on the other three sides. Hydraulic conductivities were estimated based on sediment texture documented in drilling logs. Complete details of the model, including input parameters, calibration, flow budget and travel times are described in [21].

**Figure 1.** Location of study area, orientation of the SEAM 3D model, and redox conditions observed in the Central Valley aquifer from wells located within 5 km of the line of section through the regional model.

**Figure 2.** Cross section showing the lithology, borehole resistivity logs, and Formations (modified from [19]), and or‐

Modeling the Long-Term Fate of Agricultural Nitrate in Groundwater in the San Joaquin Valley, California

http://dx.doi.org/10.5772/53652

155

ganic carbon content of sediments in the San Joaquin aquifer.

Modeling the Long-Term Fate of Agricultural Nitrate in Groundwater in the San Joaquin Valley, California http://dx.doi.org/10.5772/53652 155

**Figure 2.** Cross section showing the lithology, borehole resistivity logs, and Formations (modified from [19]), and or‐ ganic carbon content of sediments in the San Joaquin aquifer.

**Figure 1.** Location of study area, orientation of the SEAM 3D model, and redox conditions observed in the Central

Valley aquifer from wells located within 5 km of the line of section through the regional model.

154 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

#### **2.3. SEAM-3D model of the San Joaquin Valley aquifer**

The Sequential Electron Acceptor Model in three dimensions (SEAM-3D) was used to con‐ struct a quantitative mass balance between electron donors and acceptors in this study [24]. Because of the computational complexity of this mass balance, it was not feasible to model the entire area of the aquifer (Fig. 1). Rather, a cross-sectional model approximately 25 km in length (Fig. 1) was constructed from the regional model of [21] from the Tuolumne River to the San Joaquin River (Fig. 1). Only the top 11 of the 16 layers of [21] were included in the SEAM cross-sectional simulations; the bottom 5 layers were below a regional confining unit (Corcoran clay) and were likely to contain water that is too old for evaluating changes in re‐ dox conditions as a result of anthropogenic processes. Specified heads cells were used at the eastern and western boundaries of the model to approximately reproduce the configuration of the flow system prior to development. The USGS program MODPATH was used to simu‐ late groundwater travel times. MODPATH is a particle tracking post-processing model that computes three dimensional flow paths using output from groundwater flow simulations based on MODFLOW [23]. The simulated head distribution and approximate times-of-travel for the cross-sectional model are shown in Figure 3. In general, the simulated flow system reflects recharge near the Sierra Nevada in the vicinity of the Tuolumne River and discharge at the San Joaquin River. The undulations of the flowpaths shown in Figure 3 reflect litho‐ logical heterogeneities that cause variations in the distribution of hydraulic conductivity. The triangles shown on each of the flowpaths delineated in Figure 3 show the lateral dis‐ tance traveled by groundwater in 10 years. The shallowest flowpaths have travel times on the order of 50 years and the deepest on the order of hundreds of years.

Equations for the mass balance of bioavailable organic carbon and electron acceptors (EAs) assume that DOC serves as the primary electron donor (carbon/energy source) for a hetero‐ trophic microbial population in the aquifer system. Physical and biogeochemical processes incorporated in equations of transport include advection, dispersion, microbially-mediated biotransformation, rate-limited sorption and desorption. For example, the mass balance

Modeling the Long-Term Fate of Agricultural Nitrate in Groundwater in the San Joaquin Valley, California

sin , *DOC s bio DOC DOC*

¶

r

is the average pore water velocity [L T-1]; D*ij* is the tensor for

 ¶

http://dx.doi.org/10.5772/53652

*bio* is a biodegradation sink term de‐

a

b

is distance along the respective Cartesian coor‐

(1)

157

*DOC*is

(2)

*Fe*, is ex‐

equation for bioavailable DOC is given as

¶

¶

*i ij*

or sink flux [M L-3]; θ is aquifer porosity [-];*xij*

*i ij*

 ¶

 ¶

*i ij*

 ¶

 ¶

 ¶

subsurface material [M L-3];*vi*

( )

¶

¶

¶

where *EO*<sup>2</sup>

tively;*EO*<sup>2</sup>

pressed as

¶

( )

\* and *EN <sup>O</sup>*<sup>3</sup>

respectively; and *R*sin*<sup>k</sup>* ,*O*<sup>2</sup>

 ¶

( ) \*

the hydrodynamic dispersion coefficient [L2 T-1]; *R*sin*<sup>k</sup>* ,*DOC*

senting fluid sources (positive) and sinks (negative) [T-1].

¶

æ ö

*i DOC ij DOC k DOC b*


¶q

*C q C C v C <sup>D</sup> C R x xx t t*

where *CDOC* is the concentration of bioavailable DOC in the aqueous phase [M L-3]; *C*¯

the concentration of bioavailable DOC in the solid phase [M M-3]; *ρb*is the bulk density of the

pendent on the mode of respiration [M L-3 T-1]; *CDOC* \* is the DOC concentration of the source

dinate axis [L]; *t* is time [T]; and *q*s is the volume flow rate per unit volume of aquifer repre‐

2 2

*O O s bio*

\*

spectively. This treatment assumes the effects of sorption on nitrate transport are small. The

consumption of the bioavailable Fe(III) concentration [M M-3] in the solid phase, *E*¯

sink, *Bio Fe Fe dE <sup>R</sup>*

sin ,

3 3

*NO s bio NO*

and *EN <sup>O</sup>*<sup>3</sup> are the aqueous phase concentrations [M L-3] of DO and nitrate, respec‐

sin ,

\* are the DO and nitrate concentrations [M L-3] of source or sink fluxes,

 ¶

> ¶

> > ¶

*bio* are the EA biodegradation sink terms [M L-3 T-1], re‐

*dt* - = (3)

¶

Mass balance equations of the aqueous phase EAs (DO and nitrate, respectively) are

\*

2 2 2

¶q

*E E <sup>q</sup> vE D E R x xx t*

*i NO ij NO k NO*

*i O ij O kO*

¶

æ ö - + +- = ç ÷ ç ÷ è ø

*bio* and *R*sin*<sup>k</sup>* ,*NO*<sup>3</sup>

¶

æ ö - + +- = ç ÷ ç ÷ è ø

3 3 3

¶q

*E E <sup>q</sup> vE D E R x xx t*

**Figure 3.** Model-derived flowpaths and times of travel in the Central Valley aquifer. See Fig. 1 for the location of the model cross section.

Equations for the mass balance of bioavailable organic carbon and electron acceptors (EAs) assume that DOC serves as the primary electron donor (carbon/energy source) for a hetero‐ trophic microbial population in the aquifer system. Physical and biogeochemical processes incorporated in equations of transport include advection, dispersion, microbially-mediated biotransformation, rate-limited sorption and desorption. For example, the mass balance equation for bioavailable DOC is given as

**2.3. SEAM-3D model of the San Joaquin Valley aquifer**

156 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

the order of 50 years and the deepest on the order of hundreds of years.

Distance (km) <sup>25</sup> <sup>0</sup>

**Figure 3.** Model-derived flowpaths and times of travel in the Central Valley aquifer. See Fig. 1 for the location of the

Head (m)

model cross section.

Elevation (m)

45


The Sequential Electron Acceptor Model in three dimensions (SEAM-3D) was used to con‐ struct a quantitative mass balance between electron donors and acceptors in this study [24]. Because of the computational complexity of this mass balance, it was not feasible to model the entire area of the aquifer (Fig. 1). Rather, a cross-sectional model approximately 25 km in length (Fig. 1) was constructed from the regional model of [21] from the Tuolumne River to the San Joaquin River (Fig. 1). Only the top 11 of the 16 layers of [21] were included in the SEAM cross-sectional simulations; the bottom 5 layers were below a regional confining unit (Corcoran clay) and were likely to contain water that is too old for evaluating changes in re‐ dox conditions as a result of anthropogenic processes. Specified heads cells were used at the eastern and western boundaries of the model to approximately reproduce the configuration of the flow system prior to development. The USGS program MODPATH was used to simu‐ late groundwater travel times. MODPATH is a particle tracking post-processing model that computes three dimensional flow paths using output from groundwater flow simulations based on MODFLOW [23]. The simulated head distribution and approximate times-of-travel for the cross-sectional model are shown in Figure 3. In general, the simulated flow system reflects recharge near the Sierra Nevada in the vicinity of the Tuolumne River and discharge at the San Joaquin River. The undulations of the flowpaths shown in Figure 3 reflect litho‐ logical heterogeneities that cause variations in the distribution of hydraulic conductivity. The triangles shown on each of the flowpaths delineated in Figure 3 show the lateral dis‐ tance traveled by groundwater in 10 years. The shallowest flowpaths have travel times on

$$-\frac{\partial}{\partial \mathbf{x}\_i} \left( \mathbf{v}\_i \mathbf{C}\_{\mathrm{DOC}} \right) + \frac{\partial}{\partial \mathbf{x}\_i} \left( D\_{ij} \frac{\partial \mathbf{C}\_{\mathrm{DOC}}}{\partial \mathbf{x}\_j} \right) + \frac{q\_s}{\theta} \mathbf{C}\_{\mathrm{DOC}}^\* - \mathbf{R}\_{\mathrm{sink}, \mathrm{DOC}}^{\mathrm{bio}} - \rho\_b \frac{\partial \overline{\mathbf{C}}\_{\mathrm{DOC}}}{\partial t} = \frac{\partial \mathbf{C}\_{\mathrm{DOC}}}{\partial t} \tag{1}$$

where *CDOC* is the concentration of bioavailable DOC in the aqueous phase [M L-3]; *C*¯ *DOC*is the concentration of bioavailable DOC in the solid phase [M M-3]; *ρb*is the bulk density of the subsurface material [M L-3];*vi* is the average pore water velocity [L T-1]; D*ij* is the tensor for the hydrodynamic dispersion coefficient [L2 T-1]; *R*sin*<sup>k</sup>* ,*DOC bio* is a biodegradation sink term de‐ pendent on the mode of respiration [M L-3 T-1]; *CDOC* \* is the DOC concentration of the source or sink flux [M L-3]; θ is aquifer porosity [-];*xij* is distance along the respective Cartesian coor‐ dinate axis [L]; *t* is time [T]; and *q*s is the volume flow rate per unit volume of aquifer repre‐ senting fluid sources (positive) and sinks (negative) [T-1].

Mass balance equations of the aqueous phase EAs (DO and nitrate, respectively) are

$$\begin{aligned} &-\frac{\partial}{\partial \mathbf{x}\_{i}} \left( v\_{i} \mathbf{E}\_{\mathrm{O\_{2}}} \right) + \frac{\partial}{\partial \mathbf{x}\_{i}} \left( D\_{\mathrm{ij}} \frac{\partial \mathbf{E}\_{\mathrm{O\_{2}}}}{\partial \mathbf{x}\_{j}} \right) + \frac{q\_{s}}{\theta} \mathbf{E}\_{\mathrm{O\_{2}}}^{\*} - R\_{\mathrm{sink}, \mathrm{O\_{2}}}^{\mathrm{bio}} = \frac{\partial \mathbf{E}\_{\mathrm{O\_{2}}}}{\partial t} \\ &- \frac{\partial}{\partial \mathbf{x}\_{i}} \left( v\_{i} \mathbf{E}\_{\mathrm{NO\_{3}}} \right) + \frac{\partial}{\partial \mathbf{x}\_{i}} \left( D\_{\mathrm{ij}} \frac{\partial \mathbf{E}\_{\mathrm{NO\_{3}}}}{\partial \mathbf{x}\_{j}} \right) + \frac{q\_{s}}{\theta} \mathbf{E}\_{\mathrm{NO\_{3}}}^{\*} - R\_{\mathrm{sink}, \mathrm{NO\_{3}}}^{\mathrm{bio}} = \frac{\partial \mathbf{E}\_{\mathrm{NO\_{3}}}}{\partial t} \mathbf{t} \end{aligned} \tag{2}$$

where *EO*<sup>2</sup> and *EN <sup>O</sup>*<sup>3</sup> are the aqueous phase concentrations [M L-3] of DO and nitrate, respec‐ tively;*EO*<sup>2</sup> \* and *EN <sup>O</sup>*<sup>3</sup> \* are the DO and nitrate concentrations [M L-3] of source or sink fluxes, respectively; and *R*sin*<sup>k</sup>* ,*O*<sup>2</sup> *bio* and *R*sin*<sup>k</sup>* ,*NO*<sup>3</sup> *bio* are the EA biodegradation sink terms [M L-3 T-1], re‐ spectively. This treatment assumes the effects of sorption on nitrate transport are small. The consumption of the bioavailable Fe(III) concentration [M M-3] in the solid phase, *E*¯ *Fe*, is ex‐ pressed as

$$-R\_{\rm sink,Fe}^{B\dot{o}o} = \frac{d\overline{E}\_{Fe}}{dt} \tag{3}$$

Biodegradation of DOC is a function of EA availability and is described using modified Monod kinetics [24]. In summary, the overall approach is to write an equation of mass bal‐ ance for each individual solute and solid-phase constituent considered and then to solve these equations simultaneously in order to compute a true mass balance as a function of time and space.

**Parameter Initial Estimate Predevelopment**

**O2-DOC Vmax (d-1)** 0.0002 0.00068 0.0032 **NO3-DOC Vmax(d-1)** 0.0001 0.000029 0.00001 **Fe-DOC Vmax (d-1)** 0.00002 0.00002 0.000013

**Sediment organic carbon (wt%)** 0.02 0.02 0.02

**DOC solubility (mg/L)** 1 1 1

The sediments used in this study to quantify concentrations of particulate and adsorbed or‐ ganic carbon were collected by the USGS using a wire-line core barrel driven into the bot‐ tom of a borehole drilled with mud rotary coring methods. Cores were collected from two bore holes named MRWA and MREA located west and east of the city of Modesto respec‐ tively (Figure 2). The cores were collected on site, logged, and stored in core boxes prior to transporting them to a storage facility. The cores were subsampled for analysis of organic carbon approximately 4 years after collection. In the storage facility, the cores were broken in half and the center of the core subsampled for analysis in order to minimize the effects of residual drilling mud. Concentrations of sediment organic carbon (Fig. 2) were analyzed with an organic carbon analyzer (Costech Analytical Technologies, Inc., Valencia, California)

The soil science literature has extensively considered the dynamics of DOC formation and transport in soils [10], and this literature is a useful starting point for considering DOC dy‐ namics in groundwater systems. The standard conceptual model of DOC dynamics in soils [10] considers that the total amount of carbon present at any given time and place reflects DOC and POC delivery from surface sources, production of DOC from bioavailable POC, biodegradation of DOC and POC, and the adsorption/desorption of DOC on soil particles.

**O2 half saturation constant (g m-3)**

**NO3 half saturation constant (g m-3)**

**Organic carbon dissolution rate (d-1)**

**Table 1.** Kinetic model parameters used in the SEAM-3D model.

**2.5. Sediment organic carbon measurements**

with a detection limit 0.001 w% organic carbon.

**2.6. Modeling organic carbon dynamics in groundwater systems**

DOC in recharge = 1 mg/L O2/DOC Vmax = 0.0002 d-1

**observations (PEST)**

Modeling the Long-Term Fate of Agricultural Nitrate in Groundwater in the San Joaquin Valley, California

1.0 1.85 3.0

1 1 1

0.01 0.01 0.0041

**Flowpath observations (PEST)**

159

http://dx.doi.org/10.5772/53652

## **2.4. Parameter estimation**

The SEAM-3D model was initially parameterized to reproduce the approximate steady-state distribution of dissolved oxygen prior to large-scale agricultural development. A relatively high (6 mg/L) concentration of dissolved oxygen was assumed to enter the aquifer at the eastern boundary, simulating recharge from the Sierra Nevada foothills (Figure 4). In addi‐ tion, concentrations of sediment organic carbon were fixed at 0.02 weight percent (wt%) throughout the model domain, concentrations of DOC entering the model with recharge were assumed to be 1 mg/L, the half-saturation constant (Ks) and maximum oxygen utiliza‐ tion rate (Vmax) were initially set at 5.0 g m-3 and 0.002 d-1, respectively. The steady-state dis‐ tribution of dissolved oxygen given these assumptions is shown in Figure 4. Simulated dissolved oxygen concentrations ranged from 6.0 to about 1.0 mg/L and generally decreased with depth.

The next step in constraining the model was to use parameter estimation (PEST) to further refine the monod kinetic parameters [25]. First, water-quality data estimated from predevel‐ opment observations of dissolved oxygen and nitrate were used to refine the Vmax values (Table 1). In a second step, water-quality observations measured along the flowpath of the cross sectional model [20] (Figure 1) were used to constrain the Vmax values (Table 1). In gen‐ eral, the PEST approach tended to lower Vmax values for the O2-DOC, NO3-DOC, and Fe-DOC reactions. The final Vmax values as constrained by the flowpath water-quality data are listed in Table 1.

**Figure 4.** Model-derived concentrations of dissolved oxygen in the Central Valley aquifer indicated by the initial pa‐ rameter estimates. See Figure 1 for the location of the model cross section; Sediment organic carbon = 0.02 wt%


**Table 1.** Kinetic model parameters used in the SEAM-3D model.

#### **2.5. Sediment organic carbon measurements**

Biodegradation of DOC is a function of EA availability and is described using modified Monod kinetics [24]. In summary, the overall approach is to write an equation of mass bal‐ ance for each individual solute and solid-phase constituent considered and then to solve these equations simultaneously in order to compute a true mass balance as a function of

158 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

The SEAM-3D model was initially parameterized to reproduce the approximate steady-state distribution of dissolved oxygen prior to large-scale agricultural development. A relatively high (6 mg/L) concentration of dissolved oxygen was assumed to enter the aquifer at the eastern boundary, simulating recharge from the Sierra Nevada foothills (Figure 4). In addi‐ tion, concentrations of sediment organic carbon were fixed at 0.02 weight percent (wt%) throughout the model domain, concentrations of DOC entering the model with recharge were assumed to be 1 mg/L, the half-saturation constant (Ks) and maximum oxygen utiliza‐ tion rate (Vmax) were initially set at 5.0 g m-3 and 0.002 d-1, respectively. The steady-state dis‐ tribution of dissolved oxygen given these assumptions is shown in Figure 4. Simulated dissolved oxygen concentrations ranged from 6.0 to about 1.0 mg/L and generally decreased

The next step in constraining the model was to use parameter estimation (PEST) to further refine the monod kinetic parameters [25]. First, water-quality data estimated from predevel‐ opment observations of dissolved oxygen and nitrate were used to refine the Vmax values (Table 1). In a second step, water-quality observations measured along the flowpath of the cross sectional model [20] (Figure 1) were used to constrain the Vmax values (Table 1). In gen‐ eral, the PEST approach tended to lower Vmax values for the O2-DOC, NO3-DOC, and Fe-DOC reactions. The final Vmax values as constrained by the flowpath water-quality data are

Distance (km) <sup>25</sup> <sup>0</sup>

**Figure 4.** Model-derived concentrations of dissolved oxygen in the Central Valley aquifer indicated by the initial pa‐ rameter estimates. See Figure 1 for the location of the model cross section; Sediment organic carbon = 0.02 wt%

Elevation (m)

45


time and space.

with depth.

listed in Table 1.

Dissolved Oxygen (mg/L)

**2.4. Parameter estimation**

The sediments used in this study to quantify concentrations of particulate and adsorbed or‐ ganic carbon were collected by the USGS using a wire-line core barrel driven into the bot‐ tom of a borehole drilled with mud rotary coring methods. Cores were collected from two bore holes named MRWA and MREA located west and east of the city of Modesto respec‐ tively (Figure 2). The cores were collected on site, logged, and stored in core boxes prior to transporting them to a storage facility. The cores were subsampled for analysis of organic carbon approximately 4 years after collection. In the storage facility, the cores were broken in half and the center of the core subsampled for analysis in order to minimize the effects of residual drilling mud. Concentrations of sediment organic carbon (Fig. 2) were analyzed with an organic carbon analyzer (Costech Analytical Technologies, Inc., Valencia, California) with a detection limit 0.001 w% organic carbon.

#### **2.6. Modeling organic carbon dynamics in groundwater systems**

The soil science literature has extensively considered the dynamics of DOC formation and transport in soils [10], and this literature is a useful starting point for considering DOC dy‐ namics in groundwater systems. The standard conceptual model of DOC dynamics in soils [10] considers that the total amount of carbon present at any given time and place reflects DOC and POC delivery from surface sources, production of DOC from bioavailable POC, biodegradation of DOC and POC, and the adsorption/desorption of DOC on soil particles. This conceptual model, in turn, can be used to build a quantitative mass-balance model of organic carbon dynamics in groundwater systems.

DOC mobilized from surface soils is one source of bioavailable DOC to groundwater, but DOC can also be generated from POC present in aquifer material as well [11]. Of all the par‐ ticulate organic carbon (POCtot) present in an aquifer, only a fraction is available to support microbial metabolism:

$$\text{POC}\_{bio} = f \text{POC}\_{tot} \tag{4}$$

Many of the silt and clay sediments showed visible remains of roots indicating that they rep‐ resented fossilized soils, which developed in between sedimentation events. The modern surface soils contained between 0.2 and 0.5 wt% organic carbon. The alluvial fan sediments, however, showed very low concentrations of organic carbon, typically between 0.01 and 0.05 wt%. Modern fluvial sediments typically contain well in excess of 0.1 wt% organic car‐ bon. Thus, the San Joaquin alluvial fan deposits contain unusually low amounts of organic

Modeling the Long-Term Fate of Agricultural Nitrate in Groundwater in the San Joaquin Valley, California

The simulated transport of nitrate through the Central Valley aquifer, using the flowpath PEST-estimated kinetic parameters (Table 1) is shown in Figure 5. This simulation assumes that recharge coming in from the agricultural areas east of the San Joaquin River contains 20 mg/L of nitrate. This assumption is consistent with ambient conditions that existed approxi‐ mately in the year 2000 [3]. The starting point for the simulations that follow, therefore, may

25 0

Distance (km) <sup>25</sup> <sup>0</sup>

**Figure 5.** Simulated transport of nitrate in the San Joaquin aquifer for 100 years. See Fig. 1 for the location of the

After ten years of simulation, nitrate at concentrations of about 10 mg/L have moved into the shallow portions of the aquifer, which is consistent with observed nitrate contamination in this system [3]. After 50 years of simulation, the shallowest portion of the aquifer shows nitrate concentrations ranging from 15 to 20 mg/L, and nitrate has reached the discharge area near the San Joaquin River. After 100 years of simulation, nitrate concentrations in por‐

 10 years

http://dx.doi.org/10.5772/53652

161

 50 years

 100 years

carbon. That characteristic, in turn, will affect the fate and transport of nitrate.

be thought of as beginning in 2000.

45


Nitrate (mg/L)

model cross section.

Elevation (m)

where *f* is the fraction of bioavailable POC. Microbial metabolism of POC*bio* can generate ad‐ ditional DOC, and this process can be conceived of as a mass-transfer from the particulate to the dissolved phase according to the equation:

$$\rho\_b \frac{\partial \overline{\mathbb{C}}\_{\text{DOC}}}{\partial t} = -k\_{\text{DOC}} \left( \mathbb{C}\_{\text{DOC}}^{eq} - \mathbb{C}\_{\text{DOC}} \right) \tag{5}$$

Where *ρb* is the bulk density of the subsurface material [M L-3], *kDOC* is the rate constant for DOC production [T-1]; *CDOC eq* is the aqueous concentration of DOC in equilibrium with POC at any point in the system [M L-3]. The sorption-desorption of DOC onto aquifer materials such as silicate minerals, ferric oxyhydroxides, and POC, can be approximated using a sim‐ ple linear isotherm:

$$
\overline{\overline{C}}\_{D\text{OC}} = K\_d \overline{C}\_{D\text{OC}}^{eq} \tag{6}
$$

where *<sup>C</sup>*¯ *DOC* is the concentration of DOC adsorbed to aquifer material; and *Kd* is the distri‐ bution coefficient [L3 M-1]. Once the fraction of bioavailable POCtot is specified, the numerical model used in this paper uses equations 2-7 to iteratively calculate concentrations of bioa‐ vailable DOC and POC in the system as a function of time and space. This bioavailable DOC then drives the sequential utilization of electron acceptors (EAs) such as dissolved oxygen (DO), nitrate, and Fe(III). Note that this approach yields a closed mass balance, as envi‐ sioned by the model of [10], for the total amount of organic carbon stored in all three com‐ partments as a function of time. This approach builds on previous numerical methods used to simulate of redox processes in groundwater systems [26].

#### **3. Results and discussion**

The sediments analyzed for particulate organic carbon in this study (Fig. 2) were predomi‐ nantly coarse to fine-grained poorly sorted sands with inter-bedded lenses of silt and clay. Many of the silt and clay sediments showed visible remains of roots indicating that they rep‐ resented fossilized soils, which developed in between sedimentation events. The modern surface soils contained between 0.2 and 0.5 wt% organic carbon. The alluvial fan sediments, however, showed very low concentrations of organic carbon, typically between 0.01 and 0.05 wt%. Modern fluvial sediments typically contain well in excess of 0.1 wt% organic car‐ bon. Thus, the San Joaquin alluvial fan deposits contain unusually low amounts of organic carbon. That characteristic, in turn, will affect the fate and transport of nitrate.

This conceptual model, in turn, can be used to build a quantitative mass-balance model of

DOC mobilized from surface soils is one source of bioavailable DOC to groundwater, but DOC can also be generated from POC present in aquifer material as well [11]. Of all the par‐ ticulate organic carbon (POCtot) present in an aquifer, only a fraction is available to support

where *f* is the fraction of bioavailable POC. Microbial metabolism of POC*bio* can generate ad‐ ditional DOC, and this process can be conceived of as a mass-transfer from the particulate to

> ( ) *DOC eq b DOC DOC DOC <sup>C</sup> kC C*

Where *ρb* is the bulk density of the subsurface material [M L-3], *kDOC* is the rate constant for DOC production [T-1]; *CDOC eq* is the aqueous concentration of DOC in equilibrium with POC at any point in the system [M L-3]. The sorption-desorption of DOC onto aquifer materials such as silicate minerals, ferric oxyhydroxides, and POC, can be approximated using a sim‐

*DOC* is the concentration of DOC adsorbed to aquifer material; and *Kd* is the distri‐

model used in this paper uses equations 2-7 to iteratively calculate concentrations of bioa‐ vailable DOC and POC in the system as a function of time and space. This bioavailable DOC then drives the sequential utilization of electron acceptors (EAs) such as dissolved oxygen (DO), nitrate, and Fe(III). Note that this approach yields a closed mass balance, as envi‐ sioned by the model of [10], for the total amount of organic carbon stored in all three com‐ partments as a function of time. This approach builds on previous numerical methods used

The sediments analyzed for particulate organic carbon in this study (Fig. 2) were predomi‐ nantly coarse to fine-grained poorly sorted sands with inter-bedded lenses of silt and clay.

M-1]. Once the fraction of bioavailable POCtot is specified, the numerical

*bio tot POC fPOC* = (4)

=- - (5)

*eq C KC DOC d DOC* <sup>=</sup> (6)

organic carbon dynamics in groundwater systems.

160 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

the dissolved phase according to the equation:

*t*

¶

¶

to simulate of redox processes in groundwater systems [26].

r

microbial metabolism:

ple linear isotherm:

bution coefficient [L3

**3. Results and discussion**

where *<sup>C</sup>*¯

The simulated transport of nitrate through the Central Valley aquifer, using the flowpath PEST-estimated kinetic parameters (Table 1) is shown in Figure 5. This simulation assumes that recharge coming in from the agricultural areas east of the San Joaquin River contains 20 mg/L of nitrate. This assumption is consistent with ambient conditions that existed approxi‐ mately in the year 2000 [3]. The starting point for the simulations that follow, therefore, may be thought of as beginning in 2000.

**Figure 5.** Simulated transport of nitrate in the San Joaquin aquifer for 100 years. See Fig. 1 for the location of the model cross section.

After ten years of simulation, nitrate at concentrations of about 10 mg/L have moved into the shallow portions of the aquifer, which is consistent with observed nitrate contamination in this system [3]. After 50 years of simulation, the shallowest portion of the aquifer shows nitrate concentrations ranging from 15 to 20 mg/L, and nitrate has reached the discharge area near the San Joaquin River. After 100 years of simulation, nitrate concentrations in por‐ tions of the aquifer near the discharge area have risen above 10 mg/L. The results of these simulations predict that, given the kinetic parameters estimated by PEST, nitrate will pene‐ trate deeper into the aquifer for the foreseeable future. The results also suggest that concen‐ trations of nitrate in the deeper portions of the Central Valley aquifer may increase above the 10 mg/L maximum concentration limit (MCL) established by the U.S. Environmental Protection Agency.

Dissolved Oxygen

Modeling the Long-Term Fate of Agricultural Nitrate in Groundwater in the San Joaquin Valley, California

NO3 rate x 20, DOC input x 10

http://dx.doi.org/10.5772/53652

163

time (yrs) 0 20 40 60 80 100 120

Nitrate

NO3 rate x 20, DOC input x 10

NO3 rate x 20, DOC input x 10

time (yrs) 0 20 40 60 80 100 120

Dissolved Organic Carbon

time (yrs) 0 20 40 60 80 100 120

**Figure 7.** Simulated concentration changes of dissolved oxygen, nitrate, and DOC at the San Joaquin River discharge area assuming the nitrate Vmax is increased by a factor of 20 and DOC concentrations in recharge are increased to 10

Because the kinetic parameters used in the model are subject to uncertainty, the next step in the analysis assessed the sensitivity of the results to changes in those parameters. For this step, the focus was on the cells at the discharge area of the San Joaquin river at the western terminus of the cross-sectional model in layer 6. The first of these simulations assumed that the Vmax for nitrate reduction coupled to DOC oxidation was increased by a factor of 20. This scenario reflects the possibility that the PEST-derived Vmax values could be underestimated.

Dissolved Oxygen (mg/L)

Nitrate (mg/L)

DOC (mg/L)

mg/L.

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

**Figure 6.** Simulated changes in the concentrations of dissolved oxygen, nitrate, and dissolved iron at the San Joaquin River discharge area assuming the nitrate Vmax is increased by a factor of 20.

Modeling the Long-Term Fate of Agricultural Nitrate in Groundwater in the San Joaquin Valley, California http://dx.doi.org/10.5772/53652 163

Dissolved Oxygen

tions of the aquifer near the discharge area have risen above 10 mg/L. The results of these simulations predict that, given the kinetic parameters estimated by PEST, nitrate will pene‐ trate deeper into the aquifer for the foreseeable future. The results also suggest that concen‐ trations of nitrate in the deeper portions of the Central Valley aquifer may increase above the 10 mg/L maximum concentration limit (MCL) established by the U.S. Environmental

Dissolved Oxygen

NO3 rate x 20

NO3 rate x 20

time (yrs) 0 20 40 60 80 100 120

Nitrate

time (yrs) 0 20 40 60 80 100 120

Dissolved Iron

time (yrs) 0 20 40 60 80 100 120

**Figure 6.** Simulated changes in the concentrations of dissolved oxygen, nitrate, and dissolved iron at the San Joaquin

Protection Agency.

Dissolved

Nitrate (mg/L)

Dissolved Iron (mg/L) 0.00 0.01 0.02 0.03 0.04 0.05

River discharge area assuming the nitrate Vmax is increased by a factor of 20.

NO3 rate x 20

162 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

Oxygen (mg/L)

**Figure 7.** Simulated concentration changes of dissolved oxygen, nitrate, and DOC at the San Joaquin River discharge area assuming the nitrate Vmax is increased by a factor of 20 and DOC concentrations in recharge are increased to 10 mg/L.

Because the kinetic parameters used in the model are subject to uncertainty, the next step in the analysis assessed the sensitivity of the results to changes in those parameters. For this step, the focus was on the cells at the discharge area of the San Joaquin river at the western terminus of the cross-sectional model in layer 6. The first of these simulations assumed that the Vmax for nitrate reduction coupled to DOC oxidation was increased by a factor of 20. This scenario reflects the possibility that the PEST-derived Vmax values could be underestimated. The results are summarized in Figure 6. The decrease in DO concentrations reflects the in‐ crease in the the O2-DOC Vmax from 0.0002 to 0.0032 d-1 indicated by the flowpath PEST. As was the case in the simulation shown in Figure 5, nitrate is predicted to begin arriving at the discharge area after about 40 years of simulation. The maximum nitrate concentrations (~10 mg/L) predicted in the simulation of Figure 6, however, were about 20% lower than those shown in the simulation of Figure 5. So, an increase in the NO3-DOC Vmax by a factor of 20 does increase the attenuation of nitrate. However, nitrate concentrations are still predicted to increase substantially at the discharge area over time.

hypothesis-testing tool to evaluate different possible future scenarios of nitrate transport in the Central Valley aquifer. Model simulations using the PEST-derived model parameters in‐ dicate that the amount of dissolved and particulate organic carbon available in the aquifer is inadequate to consume the amount of DO that typically recharges the aquifer, and that NO3 derived from agricultural activities will be drawn deeper into the aquifer in the foreseeable future. Model simulations that increase the assumed rate of nitrate reaction with DOC by a factor of 20 decrease simulated concentrations of nitrate near the discharge area of the aqui‐

Modeling the Long-Term Fate of Agricultural Nitrate in Groundwater in the San Joaquin Valley, California

Nitrate concentrations have increased substantially in shallow groundwater of the San Joa‐ quin Valley in recent years. This phenomenon has led many public supply well operators to tap progressively deeper groundwater in order to avoid nitrate contamination. That prac‐ tice, in turn, has led to concern that elevated nitrate concentrations will continue to be drawn deeper into this groundwater system over time. The current study used a mass bal‐ ance modeling approach to assess the possible transport and attenuation of nitrate in this system. The results indicate that the amount of available organic carbon in this system, ei‐ ther DOC or particulate organic carbon in aquifer solids, is inadequate to fully attenuate the nitrate that is now entering the shallow portion of the aquifer. This finding, in turn, suggests that nitrate contamination in the Central Valley aquifer will continue to move deeper into

the system, and may eventually reach the discharge area near the San Joaquin River.

and does not constitute endorsement by the U.S. Government.

, Bruce G. Campbell1

\*Address all correspondence to: chapelle@usgs.gov

1 U.S. Geological Survey, Columbia, SC, USA

2 Virginia Tech University, Blacksburg, VA, USA

This research was supported by the National Water-Quality Assessment (NAWQA) pro‐ gram of the U.S. Geological Survey. Use of trade names is for identification purposes only

, Mark A. Widdowson2

[1] Maupin, M.A., and N.L. Barber. 2005. Estimated withdrawals prom principal aqui‐

fers in the United States, 2000. U.S. Geological Survey Circular 1270. 52 p.

and Mathew K. Landon1

http://dx.doi.org/10.5772/53652

165

**Acknowledgements**

**Author details**

Francis H. Chapelle1

**References**

fer by about 20 percent, but simulated concentrations are still substantial.

Note the simulated behavior of nitrate and dissolved iron indicated in Figure 6. As dis‐ solved oxygen concentrations decrease (due to the PEST-indicated increase in O2-DOC Vmax), iron concentrations at the discharge area begin to increase. This is consistent with observed detections of dissolved iron near the discharge area (Figure 1). However, as nitrate en‐ croaches on the discharge area, nitrate metabolism begins to replace iron metabolism and iron concentrations begin to decrease. This reflects the design of SEAM-3D which uses the most efficient available electron acceptor, in this case nitrate, preferentially to Fe(III).

In addition to increasing nitrate concentrations in recharge water, agricultural practices such as fertilization, tilling, and irrigation have the capacity to increase concentrations of DOC as well. The next simulation, therefore, explored the sensitivity of the model to increasing con‐ centrations of DOC in recharge water from 1 to 10 mg/L. The higher NO3-DOC Vmax used in Figure 6 was also used. The results are shown in Figure 7 and indicate additional attenua‐ tion of nitrate, with maximum concentrations decreasing from 10 mg/L (Figure 6) to about 6 mg/L. DOC concentrations at the discharge area also increase from 1 to 4 mg/L. This, in turn, indicates that increasing concentrations of DOC entering the aquifer with recharge may in‐ deed increase nitrate attenuation. However, the model results still indicate substantial ni‐ trate concentrations reaching the discharge area.

#### **4. Summary and conclusions**

Here, a numerical mass-balance modeling approach was used to simulate the long-term fate and transport of agriculturally-derived nitrate in the aquifer system underlying the San Joa‐ quin Valley in California. The SEAM-3D code (Sequential Electron Acceptor Model in three dimensions) used in this study couples the oxidation of dissolved organic carbon (DOC) to the reduction of dissolved oxygen (DO), nitrate (NO3), and ferric iron (Fe(III)) using monod kinetics and including inhibition functions to force the utilization of electron acceptors in the order DO> NO3>Fe(III). A cross-sectional model 25 kilometers in length was constructed by taking the hydraulic conductivity distribution from a calibrated regional model and provid‐ ing boundary conditions that approximate the historical steady-state distribution of hy‐ draulic head. The model was initially parameterized by matching model-derived concentrations of DOC and DO to historically measured concentrations. The parameteriza‐ tion was then refined using parameter estimation (PEST) on measured point concentrations of DOC, DO, NO3, and dissolved iron (Fe(II)). The parameterized model was then used as a hypothesis-testing tool to evaluate different possible future scenarios of nitrate transport in the Central Valley aquifer. Model simulations using the PEST-derived model parameters in‐ dicate that the amount of dissolved and particulate organic carbon available in the aquifer is inadequate to consume the amount of DO that typically recharges the aquifer, and that NO3 derived from agricultural activities will be drawn deeper into the aquifer in the foreseeable future. Model simulations that increase the assumed rate of nitrate reaction with DOC by a factor of 20 decrease simulated concentrations of nitrate near the discharge area of the aqui‐ fer by about 20 percent, but simulated concentrations are still substantial.

Nitrate concentrations have increased substantially in shallow groundwater of the San Joa‐ quin Valley in recent years. This phenomenon has led many public supply well operators to tap progressively deeper groundwater in order to avoid nitrate contamination. That prac‐ tice, in turn, has led to concern that elevated nitrate concentrations will continue to be drawn deeper into this groundwater system over time. The current study used a mass bal‐ ance modeling approach to assess the possible transport and attenuation of nitrate in this system. The results indicate that the amount of available organic carbon in this system, ei‐ ther DOC or particulate organic carbon in aquifer solids, is inadequate to fully attenuate the nitrate that is now entering the shallow portion of the aquifer. This finding, in turn, suggests that nitrate contamination in the Central Valley aquifer will continue to move deeper into the system, and may eventually reach the discharge area near the San Joaquin River.

## **Acknowledgements**

The results are summarized in Figure 6. The decrease in DO concentrations reflects the in‐ crease in the the O2-DOC Vmax from 0.0002 to 0.0032 d-1 indicated by the flowpath PEST. As was the case in the simulation shown in Figure 5, nitrate is predicted to begin arriving at the discharge area after about 40 years of simulation. The maximum nitrate concentrations (~10 mg/L) predicted in the simulation of Figure 6, however, were about 20% lower than those shown in the simulation of Figure 5. So, an increase in the NO3-DOC Vmax by a factor of 20 does increase the attenuation of nitrate. However, nitrate concentrations are still predicted

Note the simulated behavior of nitrate and dissolved iron indicated in Figure 6. As dis‐ solved oxygen concentrations decrease (due to the PEST-indicated increase in O2-DOC Vmax), iron concentrations at the discharge area begin to increase. This is consistent with observed detections of dissolved iron near the discharge area (Figure 1). However, as nitrate en‐ croaches on the discharge area, nitrate metabolism begins to replace iron metabolism and iron concentrations begin to decrease. This reflects the design of SEAM-3D which uses the

In addition to increasing nitrate concentrations in recharge water, agricultural practices such as fertilization, tilling, and irrigation have the capacity to increase concentrations of DOC as well. The next simulation, therefore, explored the sensitivity of the model to increasing con‐ centrations of DOC in recharge water from 1 to 10 mg/L. The higher NO3-DOC Vmax used in Figure 6 was also used. The results are shown in Figure 7 and indicate additional attenua‐ tion of nitrate, with maximum concentrations decreasing from 10 mg/L (Figure 6) to about 6 mg/L. DOC concentrations at the discharge area also increase from 1 to 4 mg/L. This, in turn, indicates that increasing concentrations of DOC entering the aquifer with recharge may in‐ deed increase nitrate attenuation. However, the model results still indicate substantial ni‐

Here, a numerical mass-balance modeling approach was used to simulate the long-term fate and transport of agriculturally-derived nitrate in the aquifer system underlying the San Joa‐ quin Valley in California. The SEAM-3D code (Sequential Electron Acceptor Model in three dimensions) used in this study couples the oxidation of dissolved organic carbon (DOC) to the reduction of dissolved oxygen (DO), nitrate (NO3), and ferric iron (Fe(III)) using monod kinetics and including inhibition functions to force the utilization of electron acceptors in the order DO> NO3>Fe(III). A cross-sectional model 25 kilometers in length was constructed by taking the hydraulic conductivity distribution from a calibrated regional model and provid‐ ing boundary conditions that approximate the historical steady-state distribution of hy‐ draulic head. The model was initially parameterized by matching model-derived concentrations of DOC and DO to historically measured concentrations. The parameteriza‐ tion was then refined using parameter estimation (PEST) on measured point concentrations of DOC, DO, NO3, and dissolved iron (Fe(II)). The parameterized model was then used as a

most efficient available electron acceptor, in this case nitrate, preferentially to Fe(III).

to increase substantially at the discharge area over time.

164 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

trate concentrations reaching the discharge area.

**4. Summary and conclusions**

This research was supported by the National Water-Quality Assessment (NAWQA) pro‐ gram of the U.S. Geological Survey. Use of trade names is for identification purposes only and does not constitute endorsement by the U.S. Government.

## **Author details**

Francis H. Chapelle1 , Bruce G. Campbell1 , Mark A. Widdowson2 and Mathew K. Landon1

\*Address all correspondence to: chapelle@usgs.gov

1 U.S. Geological Survey, Columbia, SC, USA

2 Virginia Tech University, Blacksburg, VA, USA

## **References**

[1] Maupin, M.A., and N.L. Barber. 2005. Estimated withdrawals prom principal aqui‐ fers in the United States, 2000. U.S. Geological Survey Circular 1270. 52 p.

[2] Belitz, Kenneth, Dubrovsky, N.M., Burow, Karen, Jurgens, Bryant, and Johnson, Ty‐ ler, 2003, Framework for a ground-water quality monitoring and assessment pro‐ gram for California: U.S. Geological Survey Water-Resources Investigations Report 03-4166, 78.

[16] Gu, B., J. Schmitt, Z. Chen, L. Liang, and J.F. McCarthy. 1995. Adsorption and de‐ sorption of different organic matter fractions on iron oxide. Geochimica et Cosmochi‐

Modeling the Long-Term Fate of Agricultural Nitrate in Groundwater in the San Joaquin Valley, California

http://dx.doi.org/10.5772/53652

167

[17] McMahon P.B. and F.H. Chapelle. 2008. Redox processes and the water quality of se‐

[18] Page, R.W. 1986. Geology of the fresh ground-water basin of the Central Valley, Cali‐ fornia, with texture maps and sections. U.S. Geological Survey Professional Paper

[19] Burow, K.R., Shelton, J.L., Hevesi, J.A., and Weissmann, G.S., 2004, Hydrogeologic Characterization of the Modesto Area, San Joaquin Valley, California: U.S. Geological

[20] Landon, M.K., Green, C.T., Belitz, K., Singleton, M.J., and Esser, B.K., 2011, Relations of hydrogeologic factors, groundwater reduction-oxidation conditions, and temporal and spatial distributions of nitrate, Central-Eastside San Joaquin Valley, California:

[21] Phillips, S.P. C.T. Green, K.R. Burow, J.L. Shelton, and D.L. Rewis. 2007. Simulation of multiscale ground-water flow in part of the northeastern San Joaquin Valley, Cali‐ fornia.. U.S. Geological Survey Scientific Investigations Report 2007-5009, 43 p.

[22] Harbaugh, A.W. E.R. Banta, M.C. Hill, and M.G. McDonald. 2000. MODFLOW-2000, the U.S. Geological Survey modular ground-water model-User guide to modulariza‐ tion concepts and the ground-water flow process. U.S. Geological Survey Open-File

[23] Pollock, D.W., 1994, User's guide for MODPATH/MODPATH-PLOT, version 3: A particle-tracking post-processing package for MODFLOW, the U.S. Geological Sur‐ vey finite-difference ground-water flow model: U.S. Geological Survey Open-File Re‐

[24] Waddill, D.W. and M.A. Widdowson. 1998. A three-dimensional model for subsur‐ face transport and biodegradation. ASCE Journal of Environmental Engineering,

[25] Doherty, John, 2005, PEST, model independent parameter estimation users manual,

[26] Feinstein, D.T. and M.A. Thomas. 2008. Hypothetical modeling of redox conditions within a complex ground-water flow field in a glacial setting. U.S. Geological Survey

5th edition: Watermark Numerical Computing, 336 p.

Scientific Investigations Report 2008-5066, 28 pp.

lected principal aquifer systems. Groundwater 46(2):259-285.

Survey Scientific Investigations Report 2004-5232, 54 p.

Hydrogeology Journal, 19: 1203-1224

mica Acta 59(2):219-229.

1401-C, 54 p.

Report 00-92, 121 p.

port 94–464, 249 p.

124(4), 336-344.


[16] Gu, B., J. Schmitt, Z. Chen, L. Liang, and J.F. McCarthy. 1995. Adsorption and de‐ sorption of different organic matter fractions on iron oxide. Geochimica et Cosmochi‐ mica Acta 59(2):219-229.

[2] Belitz, Kenneth, Dubrovsky, N.M., Burow, Karen, Jurgens, Bryant, and Johnson, Ty‐ ler, 2003, Framework for a ground-water quality monitoring and assessment pro‐ gram for California: U.S. Geological Survey Water-Resources Investigations Report

[3] Jurgens, B.C., K.R. Burow, B.A. Dalgish, and J.L. Shelton. 2008. Hydrogeology, water chemistry, and factors affecting the transport of contaminants in the zone of contri‐ bution of a public-supply well in Modesto, Eastern San Joaquin Valley, California.

[4] Korom, S.F., 1992. Natural denitrification in the saturated zone: a review. Water Re‐

[5] Leenheer, J.A. 1974. Occurrence of dissolved organic carbon in selected groundwater samples in the United States. U.S. Geological Survey Journal of Research 2: 361-369.

[6] Thurman, E.M. 1985. Organic Geochemistry of Natural Waters. Dordrecht/Boston/

[7] Aiken, G.R. 1989. Organic matter in groundwater. U.S. Geological Survey Open File

[8] McMahon, P.B., D.F. Williams, and J.T. Morris. 1990. Production and carbon isotopic composition of bacterial CO2 in deep Coastal Plain sediments of South Carolina.

[9] Lilienfein, J., R.G. Qualls, S.M. Uselman, and S.D. Bridgham. 2004. Adsorption of dis‐ solved organic carbon and nitrogen in soils of a weathering chronosequence. Soil Sci‐

[10] Kalbiz, K., S. Solinger, J.H. Park, B. Michalzik, B. Matzner. Controls on the dynamics of dissolved organic matter in soils: A Review. Soil Science 165(4):277-304.

[11] McMahon, P.B. and F.H. Chapelle. 1991. Microbial organic-acid production in aqui‐

[12] Kahle, M., M. Kleber, and R. Jahn. 2003. Retention of dissolved organic matter by phyllosilcate and soil clay fractions in relation to mineral properties. Organic Geo‐

[13] Davis, J.A. 1982. Adsorption of natural dissolved organic matter at the oxide/water

[14] Findlay, S. and W.V. Sobczak. 1996. Variability in removal of dissolved organic car‐ bon in hyporheic sediments. Journal of the North American Benthological Society 15

[15] Jardine, P.M. M.A. Mayes, P.J. Mulholland, P.J. Hanson, J.R. Tarver, R.J. Luxmoore, J.F. McCarthy, and G.V. Wilson. 2006. Vadose zone flow and transport of dissolved organic carbon at multiple scales in humid regimes. Vadose Zone Journal 5:140-152.

tard sediments and its role in aquifer geochemistry. Nature 349:233-235.

interface. Geochimica et Cosmochimica Acta 46:2381-2393.

U.S. Geological Survey Scientific Investigations Report 2008-5156, 78 pp.

Lancaster: Martinus Nijhoff/DR W. Junk Publishers, 497 pp.

03-4166, 78.

Report 02-89, 7pp.

ence 68:292-305.

chemistry 35:269-276.

no. 1:35-41.

Ground Water 28(5): 693-702.

sources Research, 28(6): 1657-1668.

166 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability


**Chapter 7**

**Groundwater and Contaminant Hydrology**

Groundwater and Contaminant Hydrology has a range of research relating to the transport and fate of contaminants in soils and groundwater. The scope of the center includes: 1) the development of new sampling and site characterization techniques; and 2) other improved

Contaminant hydrology is the study of processes that affect both ground and surface water pollution. It draws on the principles of hydrology and chemistry. Contaminant hydrology and water quality research seeks to understand the role of soil properties and hydrologic processes on ground and surface water pollution and develop strategies to mitigate their impacts. Research is done at all scales varying from soil pore to basin scale and covers both traditional and emerging contaminants. Groundwater and contaminant hydrology studies include fate and transport of jet fuel leakages from oil depots, producing water injection in shallow wells from the oil and gas exploration field concession areas, veterinary pharmaceuticals from landapplied manure, pathogen losses from manure application, fate and transport of disposal wastes in unlined evaporation ponds from pharmaceutical industries, impacts of tile drainage on sediment and nutrient pollution on Rivers, sediment-turbidity relationships, water quality

Some research institutes address national and international needs for subsurface contaminant characterization and remediation across a spectrum of approaches - laboratory experiments, field tests, and theoretical and numerical groundwater flow and transport investigations. Some of the developing countries most critical subsurface contamination issues, including the chemical evolution of highly alkaline radioactive waste in storage tanks; reduction, reoxidation, and diffusion of uranium forms in sediments; hydraulic properties of unsaturated gravels; and the natural production of transport-enhancing mobile nanoparticles in the

> © 2013 Ahmad et al.; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2013 Ahmad et al.; licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Zulfiqar Ahmad, Arshad Ashraf, Gulraiz Akhter and

Additional information is available at the end of the chapter

Iftikhar Ahmad

**1. Introduction**

http://dx.doi.org/10.5772/54732

groundwater remediation techniques.

modeling, and TMDL and paired watershed studies.

## **Groundwater and Contaminant Hydrology**

Zulfiqar Ahmad, Arshad Ashraf, Gulraiz Akhter and Iftikhar Ahmad

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/54732

## **1. Introduction**

Groundwater and Contaminant Hydrology has a range of research relating to the transport and fate of contaminants in soils and groundwater. The scope of the center includes: 1) the development of new sampling and site characterization techniques; and 2) other improved groundwater remediation techniques.

Contaminant hydrology is the study of processes that affect both ground and surface water pollution. It draws on the principles of hydrology and chemistry. Contaminant hydrology and water quality research seeks to understand the role of soil properties and hydrologic processes on ground and surface water pollution and develop strategies to mitigate their impacts. Research is done at all scales varying from soil pore to basin scale and covers both traditional and emerging contaminants. Groundwater and contaminant hydrology studies include fate and transport of jet fuel leakages from oil depots, producing water injection in shallow wells from the oil and gas exploration field concession areas, veterinary pharmaceuticals from landapplied manure, pathogen losses from manure application, fate and transport of disposal wastes in unlined evaporation ponds from pharmaceutical industries, impacts of tile drainage on sediment and nutrient pollution on Rivers, sediment-turbidity relationships, water quality modeling, and TMDL and paired watershed studies.

Some research institutes address national and international needs for subsurface contaminant characterization and remediation across a spectrum of approaches - laboratory experiments, field tests, and theoretical and numerical groundwater flow and transport investigations. Some of the developing countries most critical subsurface contamination issues, including the chemical evolution of highly alkaline radioactive waste in storage tanks; reduction, reoxidation, and diffusion of uranium forms in sediments; hydraulic properties of unsaturated gravels; and the natural production of transport-enhancing mobile nanoparticles in the

© 2013 Ahmad et al.; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 Ahmad et al.; licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

subsurface. Inverse modeling of reactive transport and joint hydrologic and geophysical inversion are investigated to develop new tools and approaches for estimating field-scale reactive transport parameters and characterizing contamination sites.

**•** Non aqueous phase liquids (NAPL)

spills of fuels like gasoline or jet fuel.

the related issues, implications and concerns.

**2. Case study — I**

**2.1. Introduction**

required.

**1.2. Groundwater flow and contaminant transport modeling**

Organic liquids that have densities greater than water are referred to as DNAPL (dense nonaqueous phase liquids). Nonaqueous phase liquids that have densities less than water are called LNAPLs (light nonaqueous phase liquids). Contamination by LNAPL typically involves

Groundwater and Contaminant Hydrology http://dx.doi.org/10.5772/54732 171

The preliminary steps in modeling groundwater flow and contaminant transport include development of a conceptual model, selection of a computer code, and developing model design [3]. Defining a numerical groundwater flow model is based on parameters like (a) sources and sinks of water in the field system; (b) the available data on geohydrologic system; (c) the system geometry i.e. types and extent of model layers; (d) the spatial and temporal structure of the hydraulic properties; and (e) boundary condition. The widely used MOD‐ FLOW [4] and MT3D solute transport [5] numerical codes use finite differences schemes and are considered very reliable. MODFLOW is a three-dimensional modular finite-difference model of U.S. Geological Survey widely used for the description and prediction of the behavior of groundwater system. The program uses variable grid spacing in x and y directions. Parameter estimation can be approached to find the set of parameter values that provides the best fit of model results to field observations,. At first stage, the model computes drawdown, direction of flow and hydraulic heads on each nodal point using a finite difference grid system. Using the steady-state hydraulic heads calculated by the model as the initial condition, the MT3D model is run to simulate contaminant transport in a groundwater system. Once a model is calibrated, it can be used to make predictions for management or other purposes [6].

Two case studies including i) simulated transport of jet fuel leaking into groundwater, Sindh Pakistan; and ii) deep-seated disposal of hydrocarbon exploration produced water using threedimensional contaminant transport model, Sindh Pakistan have been discussed to highlight

Groundwater is a major source for domestic and industrial uses in many urban settlements of the world. Effluents from industrial areas as well as accidental spills and leaks from surface and underground storage tanks are the main sources of natural groundwater contamination. When such contamination is detected, it becomes essential to estimate the spatial extent of contamination. Conventionally, determination of the extent of contamination is undertaken by taking many samples within time and budgetary constraints from several points, which in general requires the installation of several observation wells. As a result, the cost of such operations can be very high, especially when measurements with higher resolution are

Contaminants can migrate directly into groundwater from below-ground sources (e.g. storage tanks, pipelines) that lie within the saturated zone. Additionally contaminants can enter the groundwater system from the surface by vertical leakage through the seals around well casings, through wells abandoned without proper procedures, or as a result of contaminant disposal of improperly constructed wells [1].

#### **1.1. Governing processes of contaminant transport**

Generally three processes can be distinguished which govern the transport of contaminants in groundwater: advection, dispersion and retardation. Dispersion and density/viscosity differ‐ ences may accelerate contaminant movement, while retardation processes can slow the rate of movement. Some contamination problems involve two or more fluids. Examples include air, water and organic liquids in the unsaturated zone, or organic liquids and water in an aquifer. Tracers are useful for characterizing water flow in the saturated and unsaturated zone.

**•** Advection

The term advection refers to the movement caused by the flow of groundwater. Groundwater flow or advection is calculated based on Darcy's law. Particle tracking can be used to calculate advective transport paths [2]. Particle tracking is a numerical method by placing a particle into the flow field and numerically integrating the flow path.

**•** Dispersion

Dispersive spreading within and transverse to the main flow direction causes a gradual dilution of the contaminant plume. The dispersive spreading of a contaminant plume is due to aquifer heterogeneities. Dispersion on the macroscopic scale is caused by variations in hydraulic conductivity and porosity. Solute transport can be influenced by preferential flowpaths, arising from variations of hydraulic conductivity, at a decimetre scale.

**•** Retardation

Two major mechanisms that retard contaminant movement are sorption and biodegradation. If the sorptive process is rapid compared with the flow velocity, the solute will reach an equilibrium condition with the sorbed phase and the process can be described by an equili‐ brium sorption isotherm. The linear sorption isotherm can be described by the equation:

$$\mathbf{C}^\* = \mathbf{K}\_d \mathbf{C} \tag{1}$$

Where *C\**= mass of solute sorbed per dry unit weight of solid (mg/kg), *C*= concentration of solute in solution in equilibrium with the mass of solute sorbed onto the solid (mg/l) and *Kd* = distribution coefficient (L/kg)

**•** Non aqueous phase liquids (NAPL)

subsurface. Inverse modeling of reactive transport and joint hydrologic and geophysical inversion are investigated to develop new tools and approaches for estimating field-scale

Contaminants can migrate directly into groundwater from below-ground sources (e.g. storage tanks, pipelines) that lie within the saturated zone. Additionally contaminants can enter the groundwater system from the surface by vertical leakage through the seals around well casings, through wells abandoned without proper procedures, or as a result of contaminant

Generally three processes can be distinguished which govern the transport of contaminants in groundwater: advection, dispersion and retardation. Dispersion and density/viscosity differ‐ ences may accelerate contaminant movement, while retardation processes can slow the rate of movement. Some contamination problems involve two or more fluids. Examples include air, water and organic liquids in the unsaturated zone, or organic liquids and water in an aquifer.

The term advection refers to the movement caused by the flow of groundwater. Groundwater flow or advection is calculated based on Darcy's law. Particle tracking can be used to calculate advective transport paths [2]. Particle tracking is a numerical method by placing a particle into

Dispersive spreading within and transverse to the main flow direction causes a gradual dilution of the contaminant plume. The dispersive spreading of a contaminant plume is due to aquifer heterogeneities. Dispersion on the macroscopic scale is caused by variations in hydraulic conductivity and porosity. Solute transport can be influenced by preferential flow-

Two major mechanisms that retard contaminant movement are sorption and biodegradation. If the sorptive process is rapid compared with the flow velocity, the solute will reach an equilibrium condition with the sorbed phase and the process can be described by an equili‐ brium sorption isotherm. The linear sorption isotherm can be described by the equation:

Where *C\**= mass of solute sorbed per dry unit weight of solid (mg/kg), *C*= concentration of solute in solution in equilibrium with the mass of solute sorbed onto the solid (mg/l) and *Kd* =

\* *C KCd* <sup>=</sup> (1)

paths, arising from variations of hydraulic conductivity, at a decimetre scale.

Tracers are useful for characterizing water flow in the saturated and unsaturated zone.

reactive transport parameters and characterizing contamination sites.

170 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

disposal of improperly constructed wells [1].

**•** Advection

**•** Dispersion

**•** Retardation

distribution coefficient (L/kg)

**1.1. Governing processes of contaminant transport**

the flow field and numerically integrating the flow path.

Organic liquids that have densities greater than water are referred to as DNAPL (dense nonaqueous phase liquids). Nonaqueous phase liquids that have densities less than water are called LNAPLs (light nonaqueous phase liquids). Contamination by LNAPL typically involves spills of fuels like gasoline or jet fuel.

## **1.2. Groundwater flow and contaminant transport modeling**

The preliminary steps in modeling groundwater flow and contaminant transport include development of a conceptual model, selection of a computer code, and developing model design [3]. Defining a numerical groundwater flow model is based on parameters like (a) sources and sinks of water in the field system; (b) the available data on geohydrologic system; (c) the system geometry i.e. types and extent of model layers; (d) the spatial and temporal structure of the hydraulic properties; and (e) boundary condition. The widely used MOD‐ FLOW [4] and MT3D solute transport [5] numerical codes use finite differences schemes and are considered very reliable. MODFLOW is a three-dimensional modular finite-difference model of U.S. Geological Survey widely used for the description and prediction of the behavior of groundwater system. The program uses variable grid spacing in x and y directions. Parameter estimation can be approached to find the set of parameter values that provides the best fit of model results to field observations,. At first stage, the model computes drawdown, direction of flow and hydraulic heads on each nodal point using a finite difference grid system. Using the steady-state hydraulic heads calculated by the model as the initial condition, the MT3D model is run to simulate contaminant transport in a groundwater system. Once a model is calibrated, it can be used to make predictions for management or other purposes [6].

Two case studies including i) simulated transport of jet fuel leaking into groundwater, Sindh Pakistan; and ii) deep-seated disposal of hydrocarbon exploration produced water using threedimensional contaminant transport model, Sindh Pakistan have been discussed to highlight the related issues, implications and concerns.

## **2. Case study — I**

## **2.1. Introduction**

Groundwater is a major source for domestic and industrial uses in many urban settlements of the world. Effluents from industrial areas as well as accidental spills and leaks from surface and underground storage tanks are the main sources of natural groundwater contamination. When such contamination is detected, it becomes essential to estimate the spatial extent of contamination. Conventionally, determination of the extent of contamination is undertaken by taking many samples within time and budgetary constraints from several points, which in general requires the installation of several observation wells. As a result, the cost of such operations can be very high, especially when measurements with higher resolution are required.

A three-dimensional model of the contaminant transport was developed to predict the fate of jet fuel, which leaked from above surface storage tanks in urban site of Karachi, Pakistan. Since the tanks were situated in a sandy layer, the dissolved product entered the groundwater system and started spreading beyond the site. The modelling process comprised of steadystate simulation of the groundwater system, transient simulation of the groundwater system in the period from January 1986 through December 2015, and calibration of jet fuel that was performed in context of different parameters in groundwater system. The fuel was simulated using a modular three-dimensional finite-difference groundwater model (PMWIN) ModFlow and solute transport model (MT3D) in the 1986-2001 periods under a hypothetical scenario. After a realistic distribution of piezometric heads within the aquifer system, calibration was achieved and matched to known conditions; the solute transport component was therefore coupled to the flow. Jet fuel concentration contour maps show the expanding plume over a given time, which become almost prominent in the preceding years.

The vadose zone is contaminated with up to 1300 ppm of total petroleum hydrocarbons (TPH) within the storage site. In the previous study by MMP [8], soil samples were collected from different locations with 3 feet (1 m) below land surface (bls) and analyzed to estimate the concentration of hydrocarbon compound. Soil samples associated with the storage area have indicated higher TPH concentrations. The contaminant plume follows the hydraulic gradient

Groundwater and Contaminant Hydrology http://dx.doi.org/10.5772/54732 173

Kim and Corapcioglu [7] developed two-dimensional model to describe areal spreading and migration of light nonaqueous-phase liquids (LNAPLs) introduced into the subsurface by spills or leaks from underground storage tanks. The nonaqueous-phase liquids (NAPL) transport model was coupled with two-dimensional contaminant transport models to predict contamination of soil gas and groundwater resulting from a LNAPL migrating on the water table. Simulations were performed using the finite-difference method to study LNAPL migration and groundwater contamination. The model was applied to subsurface contami‐ nation by jet-fuel. Results indicated that LNAPL migration was affected mostly by volatiliza‐ tion. Further, the spreading and movement of the dissolved plume was affected by the geology of the area and the free-product plume. Most of the spilled mass remained as a free LNAPL phase 20 years after the spill. The migration of LNAPL for such a long period resulted in the

El-Kadi [10] investigated the US Navy's bulk fuel storage facility at Red Hill located in the island of Oahu. The facility consisted of 20 buried steel tanks with a capacity of about 12.5 million gallons each. The tanks contain jet-fuel and diesel fuel marine. The bottoms of the tanks are situated about 80 feet above the basal water table. The geology of the area is primarily basaltic lava flows. Investigations found evidence of releases from several tanks. Two borings were drilled to identify and monitor potential migration of contamination to the potable water source. A numerical model of the regional hydro‐ geology at the Red Hill Fuel Storage Facility (RHFSF) was developed to simulate the fate and transport of potential contamination from the jet-fuel tanks and the effect on the

Periago et al. [11] investigated infiltration into soil of contaminants present in cattle slurry. Column experiments were performed in order to characterize the release of contaminants at the slurry-soil interface after surface application of slurry with subsequent rainfall or irrigation. The shape of the release curves suggests that the release of substances from slurry can be modeled by a single-parameter release function. They compared prediction of solute transport (a) with input defined by the release function and (b) assuming rectangular-pulse input.

Eric et al. [12] developed a parameter identification (PI) procedure and implemented with the United States Geological Survey's Method of Characteristics (USGS-MOC) model. The test results showed that the proposed algorithm could identify transmissivity and dispersivity accurately under ideal situations. Because of the improved efficiency in model calibration,

contamination of both groundwater and a large volume of soil.

saltwater/freshwater transition zone of various pumping scenarios.

extended application to field conditions was effective.

to the southwest.

*2.1.2. Literature review*

Two-dimensional (2-D) solute transport models can be used to predict the effects of transverse dispersion of the contaminant plume (spreading). Additionally, 2-D models are appropriate where the contaminant source may lie within or near the radius of influence of a continuously pumping well. While three-dimensional (3-D) numerical models should only be used if extensive data are available regarding vertical and horizontal heterogeneity, and spatial variability in contaminant concentrations. A localized contaminant transport model for groundwater is developed to gain insight into the dynamics of the leakage of jet fuel from above-ground storage tanks in the metropolitan area of Karachi, Pakistan. Jet fuel consists of refined, kerosene-type hydrocarbons, which are mixtures of benzene, toluene, ethyl benzene, and isomers of xylene [7]. Hypothetical monitoring wells were established to estimate the concentrations of jet fuel over a stipulated time period as a result of continuous seepage from the storage depot. Although, no specific data on the history of seepage were available, in view of the results inferred from an electrical resistivity sounding survey (ERSS) regarding the nature of the subsurface lithologies coupled with the findings of previous investigators from Mott MacDonald Pakistan (MMP) [8], it is envisaged that seepage from the storage tanks occurred for more than a decade. ERSS is used to obtain the subsurface resistivity values that are assigned to different geological material. MMP [8] conducted the study on soil and groundwater assessment of environmental damage due to oil pollution and remedial measures were suggested for depots / installations / airfields of Shell Pakistan, scattered throughout the country. Appendix A provides the composition of jet-fuel (Table 1).

## *2.1.1. Site description*

The project site is located between longitudes 67º 07´ 20" and 67º 10´ 30" and latitudes 24º 52´ 20" and 24º 54´ 20". The land surface elevation ranges from approximately 15 to 33 meters above mean sea level (masl). In the far south, the Malir River drains into the Arabian Sea (Figure 1). The typical lithology of the site is silty to sandy clay from 0 to 15 ft (4.6 m) bls, gravelly sand from 15 to 43 ft (4.6 to 13 m) bls, and clayey to silty sand from 43 to 80 ft (13 to 24 m) bls. The region is arid with an average annual rainfall of about 200 mm (7.9 in). Out of this, only 10% [9] is considered to recharge the aquifer system (6.34 x 10-9 m/sec).

The vadose zone is contaminated with up to 1300 ppm of total petroleum hydrocarbons (TPH) within the storage site. In the previous study by MMP [8], soil samples were collected from different locations with 3 feet (1 m) below land surface (bls) and analyzed to estimate the concentration of hydrocarbon compound. Soil samples associated with the storage area have indicated higher TPH concentrations. The contaminant plume follows the hydraulic gradient to the southwest.

## *2.1.2. Literature review*

A three-dimensional model of the contaminant transport was developed to predict the fate of jet fuel, which leaked from above surface storage tanks in urban site of Karachi, Pakistan. Since the tanks were situated in a sandy layer, the dissolved product entered the groundwater system and started spreading beyond the site. The modelling process comprised of steadystate simulation of the groundwater system, transient simulation of the groundwater system in the period from January 1986 through December 2015, and calibration of jet fuel that was performed in context of different parameters in groundwater system. The fuel was simulated using a modular three-dimensional finite-difference groundwater model (PMWIN) ModFlow and solute transport model (MT3D) in the 1986-2001 periods under a hypothetical scenario. After a realistic distribution of piezometric heads within the aquifer system, calibration was achieved and matched to known conditions; the solute transport component was therefore coupled to the flow. Jet fuel concentration contour maps show the expanding plume over a

Two-dimensional (2-D) solute transport models can be used to predict the effects of transverse dispersion of the contaminant plume (spreading). Additionally, 2-D models are appropriate where the contaminant source may lie within or near the radius of influence of a continuously pumping well. While three-dimensional (3-D) numerical models should only be used if extensive data are available regarding vertical and horizontal heterogeneity, and spatial variability in contaminant concentrations. A localized contaminant transport model for groundwater is developed to gain insight into the dynamics of the leakage of jet fuel from above-ground storage tanks in the metropolitan area of Karachi, Pakistan. Jet fuel consists of refined, kerosene-type hydrocarbons, which are mixtures of benzene, toluene, ethyl benzene, and isomers of xylene [7]. Hypothetical monitoring wells were established to estimate the concentrations of jet fuel over a stipulated time period as a result of continuous seepage from the storage depot. Although, no specific data on the history of seepage were available, in view of the results inferred from an electrical resistivity sounding survey (ERSS) regarding the nature of the subsurface lithologies coupled with the findings of previous investigators from Mott MacDonald Pakistan (MMP) [8], it is envisaged that seepage from the storage tanks occurred for more than a decade. ERSS is used to obtain the subsurface resistivity values that are assigned to different geological material. MMP [8] conducted the study on soil and groundwater assessment of environmental damage due to oil pollution and remedial measures were suggested for depots / installations / airfields of Shell Pakistan, scattered throughout the

The project site is located between longitudes 67º 07´ 20" and 67º 10´ 30" and latitudes 24º 52´ 20" and 24º 54´ 20". The land surface elevation ranges from approximately 15 to 33 meters above mean sea level (masl). In the far south, the Malir River drains into the Arabian Sea (Figure 1). The typical lithology of the site is silty to sandy clay from 0 to 15 ft (4.6 m) bls, gravelly sand from 15 to 43 ft (4.6 to 13 m) bls, and clayey to silty sand from 43 to 80 ft (13 to 24 m) bls. The region is arid with an average annual rainfall of about 200 mm (7.9 in). Out of this, only 10%

given time, which become almost prominent in the preceding years.

172 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

country. Appendix A provides the composition of jet-fuel (Table 1).

[9] is considered to recharge the aquifer system (6.34 x 10-9 m/sec).

*2.1.1. Site description*

Kim and Corapcioglu [7] developed two-dimensional model to describe areal spreading and migration of light nonaqueous-phase liquids (LNAPLs) introduced into the subsurface by spills or leaks from underground storage tanks. The nonaqueous-phase liquids (NAPL) transport model was coupled with two-dimensional contaminant transport models to predict contamination of soil gas and groundwater resulting from a LNAPL migrating on the water table. Simulations were performed using the finite-difference method to study LNAPL migration and groundwater contamination. The model was applied to subsurface contami‐ nation by jet-fuel. Results indicated that LNAPL migration was affected mostly by volatiliza‐ tion. Further, the spreading and movement of the dissolved plume was affected by the geology of the area and the free-product plume. Most of the spilled mass remained as a free LNAPL phase 20 years after the spill. The migration of LNAPL for such a long period resulted in the contamination of both groundwater and a large volume of soil.

El-Kadi [10] investigated the US Navy's bulk fuel storage facility at Red Hill located in the island of Oahu. The facility consisted of 20 buried steel tanks with a capacity of about 12.5 million gallons each. The tanks contain jet-fuel and diesel fuel marine. The bottoms of the tanks are situated about 80 feet above the basal water table. The geology of the area is primarily basaltic lava flows. Investigations found evidence of releases from several tanks. Two borings were drilled to identify and monitor potential migration of contamination to the potable water source. A numerical model of the regional hydro‐ geology at the Red Hill Fuel Storage Facility (RHFSF) was developed to simulate the fate and transport of potential contamination from the jet-fuel tanks and the effect on the saltwater/freshwater transition zone of various pumping scenarios.

Periago et al. [11] investigated infiltration into soil of contaminants present in cattle slurry. Column experiments were performed in order to characterize the release of contaminants at the slurry-soil interface after surface application of slurry with subsequent rainfall or irrigation. The shape of the release curves suggests that the release of substances from slurry can be modeled by a single-parameter release function. They compared prediction of solute transport (a) with input defined by the release function and (b) assuming rectangular-pulse input.

Eric et al. [12] developed a parameter identification (PI) procedure and implemented with the United States Geological Survey's Method of Characteristics (USGS-MOC) model. The test results showed that the proposed algorithm could identify transmissivity and dispersivity accurately under ideal situations. Because of the improved efficiency in model calibration, extended application to field conditions was effective.

Values of hydraulic conductivities (K) for the layers are taken from the literature [14]. The K value for depth range 16-31 ft (4.9-9.4m) is taken as 1.7x10-6 ft/sec (5.2 x10-7 m/sec) and for

The transport and fate of hydrocarbons depend on multi physical and chemical process‐ es, including advection, dispersion, volatilization, dissolution, biodegradation, and sorp‐ tion. When a solute undergoes chemical reactions, its rate of movement may be substantially less than the average rate of groundwater flow. In this study, retardation of the movement of dissolved hydrocarbons is simulated as a sorption process, which in‐ cludes both adsorption and partitioning into soil organic matter or organic solvents. The MT3D software was used for simulation [5]. It uses a linear isotherm to simulate parti‐ tioning of a contaminant species between the porous media and the fluid phase due to sorption. This sorption process is approximated by the following equilibrium relationship

Where *S* is the concentration of the adsorbed phase (M/M), *C* is the concentration of the

organic materials are commonly calculated as the product of the fraction of organic carbon in the soil, *foc*, and the organic carbon partitioning coefficient, *Koc*, or *Kd* = *foc Koc. Koc* values are contaminant specific and reported in various sources [15-17]. The *foc* in the uncontaminated soil was estimated to range from 0.001 to 0.02 based on guidelines by [18]. Assuming the linear

Where rb is the bulk density of the porous material (M/L3) and neff is the effective porosity.

Processing ModFlow for Windows (PMWIN5), a modular 3-D finite-difference groundwater model, is used to configure the flow field [4]. The model consists of 41 columns and 39 rows in each layer (Figure 2). The size of cells is 410 ft x410 ft (125 m x 125 m) outside the fuel storage domain and 205 ft x 205 ft (62.5 m x 62.5 m) within the storage domain. Automatic calibration of the water table was made with algorithm - UCODE and a perfect match obtained with the known condition prior to developing the transport model [6]. Using the steady-state hydraulic heads calculated by PMWIN5 as the initial condition, the solute transport model MT3D was run to simulate the dispersion of the dissolved jet-fuel plume [5]. The parameters adjusted were the retardation factor R for each cell within the finite-difference grid, and the dispersion

dissolved phase (M/L3) and *Kd* is the sorption or distribution coefficient (L3

isotherm, the retardation factor (*R*) is expressed as follows:

For jet fuel, the distribution coefficient Kd is taken as 0.004415 ft3

*R* = 1 + (48 kg/ft3)/0.25 *x* 0.004415 = 1.848

*2.2.1. Numerical ground water flow modeling*

*<sup>d</sup> S KC* = (2)

( ) 1 / *R rn K b eff d* = + (3)

/M). *Kd* values for

Groundwater and Contaminant Hydrology http://dx.doi.org/10.5772/54732 175

/kg. With these values [7],

depth range 31-97 ft (9.4-30 m) as 1.5x10-5 ft/sec (4.6 x10-6 m/sec).

between the dissolved and adsorbed phases:

**Figure 1.** Location of the study area in Southern Pakistan

Jin et al. [13] investigated hydrocarbon plumes in groundwater through the installation of extensive monitoring wells. Electromagnetic induction survey was carried out as an alternative technique for mapping petroleum contaminants in the subsurface. The surveys were conduct‐ ed at a coal mining site near Gillette, Wyoming, using the EM34-XL ground conductivity meter. Data from this survey used to validate with the known concentrations of diesel compounds detected in groundwater. Groundwater data correlated perfectly with the electromagnetic survey data, which was used to generate a site model to identify subsurface diesel plumes. Results from this study indicated that this geophysical technique was an effective tool for assessing subsurface petroleum hydrocarbon sources and plumes at contaminated sites.

#### **2.2. Model conceptualization and simulation**

The groundwater flow system is treated as a two layer. The upper layer (predominantly silty sand) is bounded above by the water table and is 15 ft (4.6 m) thick, while the lower layer (predominantly clayey sand) is 65 ft (20 m) thick. These are unconfined and recharged from the surface by infiltrating rain, but only over permeable surfaces. A small stream runs along east of model area acts as a drain to groundwater, which flows from the northeast to the southwest. With the exception of the stream, all other boundaries are artificial, that is neither constant head, nor constant flow boundaries. The processes that control the groundwater flow are: (i) recharges from infiltrated rainfall; (ii) flow entering the model across the eastern boundary (also across the northern boundary); (iii) flow reaching the stream; (iv) flow leaving the model across the southern boundary; and (v) pumping from one well near tank no. 9. Values of hydraulic conductivities (K) for the layers are taken from the literature [14]. The K value for depth range 16-31 ft (4.9-9.4m) is taken as 1.7x10-6 ft/sec (5.2 x10-7 m/sec) and for depth range 31-97 ft (9.4-30 m) as 1.5x10-5 ft/sec (4.6 x10-6 m/sec).

The transport and fate of hydrocarbons depend on multi physical and chemical process‐ es, including advection, dispersion, volatilization, dissolution, biodegradation, and sorp‐ tion. When a solute undergoes chemical reactions, its rate of movement may be substantially less than the average rate of groundwater flow. In this study, retardation of the movement of dissolved hydrocarbons is simulated as a sorption process, which in‐ cludes both adsorption and partitioning into soil organic matter or organic solvents. The MT3D software was used for simulation [5]. It uses a linear isotherm to simulate parti‐ tioning of a contaminant species between the porous media and the fluid phase due to sorption. This sorption process is approximated by the following equilibrium relationship between the dissolved and adsorbed phases:

$$S = K\_d \mathbb{C} \tag{2}$$

Where *S* is the concentration of the adsorbed phase (M/M), *C* is the concentration of the dissolved phase (M/L3) and *Kd* is the sorption or distribution coefficient (L3 /M). *Kd* values for organic materials are commonly calculated as the product of the fraction of organic carbon in the soil, *foc*, and the organic carbon partitioning coefficient, *Koc*, or *Kd* = *foc Koc. Koc* values are contaminant specific and reported in various sources [15-17]. The *foc* in the uncontaminated soil was estimated to range from 0.001 to 0.02 based on guidelines by [18]. Assuming the linear isotherm, the retardation factor (*R*) is expressed as follows:

$$R = 1 + \left(r\_b / n\_{\rm eff} \right) K\_d \tag{3}$$

Where rb is the bulk density of the porous material (M/L3) and neff is the effective porosity. For jet fuel, the distribution coefficient Kd is taken as 0.004415 ft3 /kg. With these values [7], *R* = 1 + (48 kg/ft3)/0.25 *x* 0.004415 = 1.848

#### *2.2.1. Numerical ground water flow modeling*

**Figure 1.** Location of the study area in Southern Pakistan

174 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

**2.2. Model conceptualization and simulation**

Jin et al. [13] investigated hydrocarbon plumes in groundwater through the installation of extensive monitoring wells. Electromagnetic induction survey was carried out as an alternative technique for mapping petroleum contaminants in the subsurface. The surveys were conduct‐ ed at a coal mining site near Gillette, Wyoming, using the EM34-XL ground conductivity meter. Data from this survey used to validate with the known concentrations of diesel compounds detected in groundwater. Groundwater data correlated perfectly with the electromagnetic survey data, which was used to generate a site model to identify subsurface diesel plumes. Results from this study indicated that this geophysical technique was an effective tool for assessing subsurface petroleum hydrocarbon sources and plumes at contaminated sites.

The groundwater flow system is treated as a two layer. The upper layer (predominantly silty sand) is bounded above by the water table and is 15 ft (4.6 m) thick, while the lower layer (predominantly clayey sand) is 65 ft (20 m) thick. These are unconfined and recharged from the surface by infiltrating rain, but only over permeable surfaces. A small stream runs along east of model area acts as a drain to groundwater, which flows from the northeast to the southwest. With the exception of the stream, all other boundaries are artificial, that is neither constant head, nor constant flow boundaries. The processes that control the groundwater flow are: (i) recharges from infiltrated rainfall; (ii) flow entering the model across the eastern boundary (also across the northern boundary); (iii) flow reaching the stream; (iv) flow leaving the model across the southern boundary; and (v) pumping from one well near tank no. 9.

Processing ModFlow for Windows (PMWIN5), a modular 3-D finite-difference groundwater model, is used to configure the flow field [4]. The model consists of 41 columns and 39 rows in each layer (Figure 2). The size of cells is 410 ft x410 ft (125 m x 125 m) outside the fuel storage domain and 205 ft x 205 ft (62.5 m x 62.5 m) within the storage domain. Automatic calibration of the water table was made with algorithm - UCODE and a perfect match obtained with the known condition prior to developing the transport model [6]. Using the steady-state hydraulic heads calculated by PMWIN5 as the initial condition, the solute transport model MT3D was run to simulate the dispersion of the dissolved jet-fuel plume [5]. The parameters adjusted were the retardation factor R for each cell within the finite-difference grid, and the dispersion coefficient. Concentration-time curves have been calculated for ten monitoring wells. PMPATH [19] is used to retrieve the groundwater flow model and simulation result from PMWIN5. A semi-analytical particle-tracking scheme is used to calculate the groundwater flow paths, travel times, and time-related capture zones resulting from pumping a neighboring well at the storage facility [20]. As a preprocessor to modeling and creating input data files, the PMWIN5 utility package was used. Prior to initiating the modeling work, a groundwater information system was established with all data in binary and / or ASCII files that could be exported to other softwares.

EJHD Shell Depot

Tank# 8 Tank# 9

General flow direction

The model was calibrated for steady-state conditions. Since it was speculated that the seepage of the fuel might have started as early as 1986, simulation of the groundwater flow was begun in January 1971, the opening of storage facility. In the steady-state phase the only input comes from constant head boundaries (along the east and west of the model) and from infiltrated rainfall. All output goes into constant head boundaries (the stream and the southwest boun‐ dary of the model). The differences in the known and simulated heads were calibrated to less than 0.30 ft (0.091 m) by making slight adjustments in the *K* values of both the layers. Transient calibration of groundwater flow was accomplished using the time-variant hydraulic head values. Parameters such as recharge rates during each stress period, hydraulic heads in the stream and along the model boundaries, aquifer storage properties, pumping rates, and timedependent capture zone were adjusted during the calibration. To be objective and consistent, the recharge from infiltration was made equal to 10% of rainfall in each month. Effective porosity of the aquifer was varied between 15% and 25% until the value of 25% was determined

Tank# 10

Pum ping W ell

**Figure 3.** Layout of the storage tanks marked on a finite-difference grid

*2.2.3. Model calibration*

Tank# 7

Source of Seepage

Groundwater and Contaminant Hydrology http://dx.doi.org/10.5772/54732 177

Jet Fuel Storage Tanks (10 million liters)

Scale: 1 small square = 205 ft by 205 ft

## *2.2.2. Locations of hypothetical wells*

The dissolved phase jet-fuel plume was traced using a combination of ten hypothetical monitoring wells (Figure 2) known as MW-1 through MW-10. The wells served to identify lithology, observe water levels, and monitor concentrations of organic compounds. The wells extend to a depth of 80 ft (24 m). In addition, actual well was completed to a depth of 100 ft (30.5 m) near storage tank no.9 in case of emergency need. In the modeling study, this well was used to track the time-related capture zone. The general layout of the storage tanks over the finite-difference grid is shown in Figure 3. The location of the pumping well is marked as a small red square in Figure 2 and Figure 3.

**Figure 2.** Model design indicating finite difference grid and locations of hypothetical monitoring wells

Scale: 1 small square = 205 ft by 205 ft

**Figure 3.** Layout of the storage tanks marked on a finite-difference grid

## *2.2.3. Model calibration*

coefficient. Concentration-time curves have been calculated for ten monitoring wells. PMPATH [19] is used to retrieve the groundwater flow model and simulation result from PMWIN5. A semi-analytical particle-tracking scheme is used to calculate the groundwater flow paths, travel times, and time-related capture zones resulting from pumping a neighboring well at the storage facility [20]. As a preprocessor to modeling and creating input data files, the PMWIN5 utility package was used. Prior to initiating the modeling work, a groundwater information system was established with all data in binary and / or ASCII files that could be

176 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

The dissolved phase jet-fuel plume was traced using a combination of ten hypothetical monitoring wells (Figure 2) known as MW-1 through MW-10. The wells served to identify lithology, observe water levels, and monitor concentrations of organic compounds. The wells extend to a depth of 80 ft (24 m). In addition, actual well was completed to a depth of 100 ft (30.5 m) near storage tank no.9 in case of emergency need. In the modeling study, this well was used to track the time-related capture zone. The general layout of the storage tanks over the finite-difference grid is shown in Figure 3. The location of the pumping well is marked as

**Figure 2.** Model design indicating finite difference grid and locations of hypothetical monitoring wells

exported to other softwares.

*2.2.2. Locations of hypothetical wells*

a small red square in Figure 2 and Figure 3.

The model was calibrated for steady-state conditions. Since it was speculated that the seepage of the fuel might have started as early as 1986, simulation of the groundwater flow was begun in January 1971, the opening of storage facility. In the steady-state phase the only input comes from constant head boundaries (along the east and west of the model) and from infiltrated rainfall. All output goes into constant head boundaries (the stream and the southwest boun‐ dary of the model). The differences in the known and simulated heads were calibrated to less than 0.30 ft (0.091 m) by making slight adjustments in the *K* values of both the layers. Transient calibration of groundwater flow was accomplished using the time-variant hydraulic head values. Parameters such as recharge rates during each stress period, hydraulic heads in the stream and along the model boundaries, aquifer storage properties, pumping rates, and timedependent capture zone were adjusted during the calibration. To be objective and consistent, the recharge from infiltration was made equal to 10% of rainfall in each month. Effective porosity of the aquifer was varied between 15% and 25% until the value of 25% was determined to be the best predictor for the model. Constant pumping rates of 1500 US gallons/hr (1.58 x 10-3 m3 /sec) and 1000 US gallons/hr (1.05 x 10-3 m3 /sec) were used in layer 1 and layer 2 respectively. The period from 1986 through 2000 was divided into seven stress periods, each of 2 years in duration. From 2000 to 2001, one stress period was assigned. The length of a stress period was made equal to the number of days in that month.

#### **2.3. Results and discussion**

The effect of the pumping well is clearly visible as a cone of depression (Figure 4). The drawdown was determined to be 14.0 ft (4.3 m) near the storage facility.

**Figure 5.** Capture zone of the pumping well with arrows indicating flow directions

Groundwater and Contaminant Hydrology http://dx.doi.org/10.5772/54732 179

**Figure 6.** Days-capture zone calculated by PMPATH

**Figure 4.** Cone of Depression visible around pumping well developed in layer 1

The model assumes uniform recharge from infiltrated rainfall to every "recharging" cell. Although effective porosity, hydraulic conductivity, and recharge may vary in space and time, the model is expected to have produced a reasonable configuration of the groundwater flow pattern throughout the whole period of simulation. The time-related capture zones produced due to constant pumping are shown in Figure 5 and Figure 6. Water balances, which was calculated for each year of the simulated period, showed a perfect match between the input and output components.

**Figure 5.** Capture zone of the pumping well with arrows indicating flow directions

to be the best predictor for the model. Constant pumping rates of 1500 US gallons/hr (1.58 x

respectively. The period from 1986 through 2000 was divided into seven stress periods, each of 2 years in duration. From 2000 to 2001, one stress period was assigned. The length of a stress

The effect of the pumping well is clearly visible as a cone of depression (Figure 4). The

/sec) were used in layer 1 and layer 2

/sec) and 1000 US gallons/hr (1.05 x 10-3 m3

178 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

period was made equal to the number of days in that month.

**Figure 4.** Cone of Depression visible around pumping well developed in layer 1

and output components.

The model assumes uniform recharge from infiltrated rainfall to every "recharging" cell. Although effective porosity, hydraulic conductivity, and recharge may vary in space and time, the model is expected to have produced a reasonable configuration of the groundwater flow pattern throughout the whole period of simulation. The time-related capture zones produced due to constant pumping are shown in Figure 5 and Figure 6. Water balances, which was calculated for each year of the simulated period, showed a perfect match between the input

drawdown was determined to be 14.0 ft (4.3 m) near the storage facility.

10-3 m3

**2.3. Results and discussion**

**Figure 6.** Days-capture zone calculated by PMPATH

## *2.3.1. Calibration of plume dispersion*

For simulation of the movement of dissolved jet fuel, lateral hydraulic conductivity values equal to 0.864 ft/day (0.263 m/day) for layer 1 and 1.29 ft/day (0.393 m/day) for layer 2 were accepted, while the vertical hydraulic conductivity were taken as 0.0864 ft/day (0.0263 m/ day) and 0.129 ft/day (0.0393 m/day) for each layer, respectively [14]. The hydraulic gradient and flow-net were obtained by running the flow component of the model derived from wa‐ ter level information in the previous study.

Using steady-state hydraulic heads as initial conditions, the evolution of the plume was modeled over nine stress periods as a result of continuous seepage from cells (18,16,1; tank7),

 2 6.30 x 107 1986 – 88 2 12.60 x 107 1988 – 90 2 18.92 x 107 1990 – 92 2 25.23 x 107 1990 – 94 2 3.15 x 108 1994 – 96 2 3.78 x 108 1996 - 98 2 4.41 x 108 1998 –00 1 4.73 x 108 2000 – 01 14 9.14 x 108 2001 – 15

Parameters describing various processes are used after calibration with different combination

**Elapsed Time (sec)**

**Period**

Groundwater and Contaminant Hydrology http://dx.doi.org/10.5772/54732 181

(19,16,1; tank 10), (19,17,1; tank 8), and (20,17,1; tank9) as shown in Table 2.

**Time interval (years)**

**Table 2.** Stress period used in time-dependent solute transport modeling of jet fuel

**Parameters Value** Longitudinal Dispersion 10 ft Transverse Dispersion 3.3 ft

Effective Porosity (Layer-1) 0.25 Effective Porosity (Layer-2) 0.30

**Table 3.** Preliminary and final values of parameters used in modeling

Molecular Diffusion 0.001 ft2/day Distribution Coefficient 0.004415 ft3/kg Retardation Factor (*R*) 1.80 to 1.848 Decay Coefficient 1 x10-9 day--1

Hydraulic Conductivity *K* (Layer-1) 1x10-5 ft/sec (0.864 ft/day) Hydraulic Conductivity *K* (Layer-2) 1.49x10-5 ft/sec (1.29 ft/day)

**Stress Period**

*2.3.3. Calibration scenario*

of parameters (Table 3).

The United States Environmental Protection Agency (USEPA) and the Georgia Environmen‐ tal Protection Division (GAEPD) recommend that the value for longitudinal dispersion should be one-tenth of the distance from the place where a contaminant enters the ground‐ water system to the down-gradient receptor (a well, stream, or other point of compliance). The distance from the storage facility (tank no. 7) to the pumping well is approximately 100 ft (30.48 m). In all calibration runs, as recommended, the value for longitudinal dispersion was set at 10 ft (3 m). USEPA and GAEPD also recommend for a solute transport model that the value for transverse dispersion equal one-third for the longitudinal dispersion. For this model, transverse dispersion would equal 3.3 ft (1.0 m). In the simulation of the fate of jet fuel, the transverse dispersion coefficient was varied within a range of 2.0 ft to 3.3 ft (0.61 m to 1.0 m). In the model a value of 0.001 ft2 /day (1 x 10-5 cm2 /sec) was used for molecular dif‐ fusion. With a retardation factor of 1.80, dissolved jet fuel takes 1.33 years to travel a lateral distance of about 70 to 80 ft (21 to 24 m) in groundwater beneath tank no. 8. The best value of the microbial decay coefficient for jet fuel is estimated to be 1-10 / day with a microbial yield coefficient for oxygen of 0.52 [14].

#### *2.3.2. Strategy development for release of Jet fuel*

The previous integrity test run on the storage tanks containing 10 million liters of jet fuel indicated no loss. The date when the leak initially began is unknown, although inventory records indicated that the leak was not present before tank integrity testing. The product has been detected in several hypothetical-monitoring wells (notably in MW-2 and MW-4) and in many soil samples taken within several tens of feet of the tank. The initial concentration of jet fuel entering the system is not of prime concern for the modeling. The product of the in‐ flux (in *L3 /T*) and the concentration (in *M/L3* ) gives the total mass of jet fuel entering the sys‐ tem in a certain time interval. For the purpose of calibrating the jet fuel input, the initial concentration used, based upon field data [8], varied from 0.095 to 0.19 g/ft3 (0.0027 to 0.0054 g/m3 ). The initial mass of jet fuel, as simulated by the model, was equal to each of four cells "injecting" at a mass rate of 95 to 190 g/ft3 (2.7 to 5.4 g/m3 ) following the initial period of 15 years during which no groundwater contamination was assumed (Table 1).


**Table 1.** Strategy developed for the plume modeling scenarios

Using steady-state hydraulic heads as initial conditions, the evolution of the plume was modeled over nine stress periods as a result of continuous seepage from cells (18,16,1; tank7), (19,16,1; tank 10), (19,17,1; tank 8), and (20,17,1; tank9) as shown in Table 2.


**Table 2.** Stress period used in time-dependent solute transport modeling of jet fuel

#### *2.3.3. Calibration scenario*

*2.3.1. Calibration of plume dispersion*

ter level information in the previous study.

180 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

to 1.0 m). In the model a value of 0.001 ft2

yield coefficient for oxygen of 0.52 [14].

flux (in *L3*

g/m3

*2.3.2. Strategy development for release of Jet fuel*

"injecting" at a mass rate of 95 to 190 g/ft3

**Table 1.** Strategy developed for the plume modeling scenarios

*/T*) and the concentration (in *M/L3*

For simulation of the movement of dissolved jet fuel, lateral hydraulic conductivity values equal to 0.864 ft/day (0.263 m/day) for layer 1 and 1.29 ft/day (0.393 m/day) for layer 2 were accepted, while the vertical hydraulic conductivity were taken as 0.0864 ft/day (0.0263 m/ day) and 0.129 ft/day (0.0393 m/day) for each layer, respectively [14]. The hydraulic gradient and flow-net were obtained by running the flow component of the model derived from wa‐

The United States Environmental Protection Agency (USEPA) and the Georgia Environmen‐ tal Protection Division (GAEPD) recommend that the value for longitudinal dispersion should be one-tenth of the distance from the place where a contaminant enters the ground‐ water system to the down-gradient receptor (a well, stream, or other point of compliance). The distance from the storage facility (tank no. 7) to the pumping well is approximately 100 ft (30.48 m). In all calibration runs, as recommended, the value for longitudinal dispersion was set at 10 ft (3 m). USEPA and GAEPD also recommend for a solute transport model that the value for transverse dispersion equal one-third for the longitudinal dispersion. For this model, transverse dispersion would equal 3.3 ft (1.0 m). In the simulation of the fate of jet fuel, the transverse dispersion coefficient was varied within a range of 2.0 ft to 3.3 ft (0.61 m

/day (1 x 10-5 cm2

fusion. With a retardation factor of 1.80, dissolved jet fuel takes 1.33 years to travel a lateral distance of about 70 to 80 ft (21 to 24 m) in groundwater beneath tank no. 8. The best value of the microbial decay coefficient for jet fuel is estimated to be 1-10 / day with a microbial

The previous integrity test run on the storage tanks containing 10 million liters of jet fuel indicated no loss. The date when the leak initially began is unknown, although inventory records indicated that the leak was not present before tank integrity testing. The product has been detected in several hypothetical-monitoring wells (notably in MW-2 and MW-4) and in many soil samples taken within several tens of feet of the tank. The initial concentration of jet fuel entering the system is not of prime concern for the modeling. The product of the in‐

tem in a certain time interval. For the purpose of calibrating the jet fuel input, the initial

**Phase Stress period Condition**

). The initial mass of jet fuel, as simulated by the model, was equal to each of four cells

(2.7 to 5.4 g/m3

concentration used, based upon field data [8], varied from 0.095 to 0.19 g/ft3

years during which no groundwater contamination was assumed (Table 1).

Safe Period 15 years (1971 to 1986) No leakage

Hazardous Period 10 years (1986 to 1996) Low to moderate leakage Risk Assessment 5 years (1996-2001) Moderate leakage Future Prediction 14 years (2001 – 2015) Accretion in leakage

/sec) was used for molecular dif‐

) gives the total mass of jet fuel entering the sys‐

) following the initial period of 15

(0.0027 to 0.0054

Parameters describing various processes are used after calibration with different combination of parameters (Table 3).


**Table 3.** Preliminary and final values of parameters used in modeling

The release of the jet fuel is simulated in four cells, all along columns 18 to 20 from row 16 to row 17. The area of injection is equal to 42025 ft2 (3,904 m2 ). The concentration of jet fuel at the source (95 to 190 mg/ft3 [2.7 to 5.4 g/m3 ]) maintained constant throughout the designated "hazardous period" simulation period (1986-1996). The concentration was increased slightly (about 0.01 %) from 1996 through 2011 and further up to longer time duration of 4 years i.e., up to 2015. The plume simulations are shown in Figure 7.

The modeled concentration at MW-1, MW-2, and MW-6 was much higher than the concen‐ trations in the remaining wells. The maximum level of concentrations recorded in MW-6 was

Groundwater and Contaminant Hydrology http://dx.doi.org/10.5772/54732 183

The shape of the plume is elliptical, with the major axis in the direction of groundwater flow. This shape results from advection and longitudinal dispersion. The lateral spread of the plume results from transverse dispersion and molecular diffusion. Upgradient spread of the plume results from molecular diffusion [25-26]. The plume travels toward the stream, which is still far away in the west. By the end of 2015 the effect of the plume becomes evident and monitoring wells MW-4, MW-5, and MW-6 indicated increased concentration of jet fuel (Figure 10). The resultant plume appears to be spreading more in the elliptical path but in the direction of

. Figure 9 reflects the plume spreading of year 2001.

2990 μg/ft3

groundwater flow.

**Figure 9.** Extent of simulated plume in 2001

**Figure 7.** Simulated Jet-fuel plumes(1986 to1988 and 1988 to 1990)

The jet-fuel break-through curves for the hypothetical monitoring wells are shown in Figure 8. Conventionally, determination of the extent and level of contamination is undertaken by taking multiple measurements in wells [21-23]. However, higher spatial resolution generally requires installation of monitoring wells, which is costly [24].

**Figure 8.** Concentration versus time based on data from 10 monitoring wells

The modeled concentration at MW-1, MW-2, and MW-6 was much higher than the concen‐ trations in the remaining wells. The maximum level of concentrations recorded in MW-6 was 2990 μg/ft3 . Figure 9 reflects the plume spreading of year 2001.

The shape of the plume is elliptical, with the major axis in the direction of groundwater flow. This shape results from advection and longitudinal dispersion. The lateral spread of the plume results from transverse dispersion and molecular diffusion. Upgradient spread of the plume results from molecular diffusion [25-26]. The plume travels toward the stream, which is still far away in the west. By the end of 2015 the effect of the plume becomes evident and monitoring wells MW-4, MW-5, and MW-6 indicated increased concentration of jet fuel (Figure 10). The resultant plume appears to be spreading more in the elliptical path but in the direction of groundwater flow.

**Figure 9.** Extent of simulated plume in 2001

The release of the jet fuel is simulated in four cells, all along columns 18 to 20 from row 16 to

"hazardous period" simulation period (1986-1996). The concentration was increased slightly (about 0.01 %) from 1996 through 2011 and further up to longer time duration of 4 years i.e.,

The jet-fuel break-through curves for the hypothetical monitoring wells are shown in Figure 8. Conventionally, determination of the extent and level of contamination is undertaken by taking multiple measurements in wells [21-23]. However, higher spatial resolution generally

(3,904 m2

). The concentration of jet fuel at the

]) maintained constant throughout the designated

row 17. The area of injection is equal to 42025 ft2

[2.7 to 5.4 g/m3

182 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

up to 2015. The plume simulations are shown in Figure 7.

**Figure 7.** Simulated Jet-fuel plumes(1986 to1988 and 1988 to 1990)

requires installation of monitoring wells, which is costly [24].

**Figure 8.** Concentration versus time based on data from 10 monitoring wells

source (95 to 190 mg/ft3

biological origin, but now it is used to describe materials most of which contain carbon and hydrogen and which may contain oxygen, nitrogen, the halogens, and lesser amounts of other elements. The simplest of these are the hydrocarbons, molecules of hydrogen and carbon, many of which are the components of natural gas, petroleum, and coal. Petroleum, however, has a very large number of components ranging from methane to the high molecular weight materials asphalt and paraffin. Typical fractions into which crude oil is separated in an oil

refinery and some principal molecular species are shown in Table 4.

Gas Below 20o C Gas

**Boiling range Product of secondary**

Naphtha 20o – 175o C Naphtha gasoline C 11H24 C 18H38

Kerosene 175o – 400o C Kerosene diesel fuel C 11H24 C 18H38

**Table 4.** The Fractions and Representative Components obtained from Crude Oil

Residue above 400o C Asphalt Heavier hydrocarbons

It is important to maintain the existing quality of groundwater because once contamination occurs; it is sometime difficult or rather impossible to clean the aquifer. There is high proba‐ bility associated with certain landuses like agriculture, industrial/urban land and drainage wells for contaminating the groundwater. The early detection of such contamination can be executed through proper monitoring of the groundwater quality. Many oil & gas companies are disposing off their producing water into the deep Ranikot formation in Bhit oil-field area of southern Pakistan. The producing water contains Total dissolved solid (TDS) within range of 18,000 - 22,000 mg/l besides oil condensate. There were concerns that the producing water is affecting the fresh water aquifers belonging to the overlying Nari and Kirthar formations. This phenomenon has been studied by the utilization of groundwater contaminant transport model. Injection has been monitored at 2100 meters depth in the Pab sandstone formation. A three-dimensional contaminant transport model was developed to simulate and monitor the migration of disposal of hydrocarbon exploration produced water in Injection well at 2000 meters depth in the Upper Cretaceous Pab sandstone in the study area. Framework of regional

**treatment**

Liquefied Pet. Gas (LPG)

Lubricating oil C 15H32 C 40H82

**Typical molecular components**

Groundwater and Contaminant Hydrology http://dx.doi.org/10.5772/54732 185

CH4 methane, C2H6

ethane C 3H8 propane, C 4H10 butane

**Fraction from distillation**

**3. Case study — II**

**3.1. Introduction**

**Figure 10.** Extent of plume spreading in year 2015

#### **2.4. Conclusions and recommendation**

Based on the modeling, it is concluded that the jet-fuel plume has neither expanded nor moved considerably. It is less than 250 ft (76.2 m) beyond the storage tanks and is oriented northeast to southwest. The level of concentration found in the simulated monitoring wells is significant, but because groundwater is brackish and thus unlikely to be used, no harmful effects are expected. However, with the continuous process of leaking jet-fuel from the storage depot, the level of concentration is expected to increase over the period of time. Regionally, the jet-fuel expansion will have on prominent effect over a longer areal coverage and will be confined to a localized area.

An interdisciplinary investigation of the processes controlling the fate and transport hydro‐ carbons in the subsurface is needed. Concentrations should be performed in wells within and down-gradient of the plume, as field data would help develop a stronger argument for the fate of jet fuel in groundwater. Periodic observations need to be carried out in wells to have good control on the changes of groundwater chemistry. Defective storage depots need to be mended to stop the release of jet-fuel in future.

## **2.5. Appendix A**

## Composition of contaminant (Jet fuel)

Knowledge of the geochemistry of a contaminated aquifer is important to understand the chemical and biological processes controlling the migration of hydrocarbon contaminants in the subsurface. Originally, the jet fuel (kerosene oil) is the name assigned to a material with a biological origin, but now it is used to describe materials most of which contain carbon and hydrogen and which may contain oxygen, nitrogen, the halogens, and lesser amounts of other elements. The simplest of these are the hydrocarbons, molecules of hydrogen and carbon, many of which are the components of natural gas, petroleum, and coal. Petroleum, however, has a very large number of components ranging from methane to the high molecular weight materials asphalt and paraffin. Typical fractions into which crude oil is separated in an oil refinery and some principal molecular species are shown in Table 4.


**Table 4.** The Fractions and Representative Components obtained from Crude Oil

## **3. Case study — II**

#### **3.1. Introduction**

**Figure 10.** Extent of plume spreading in year 2015

184 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

**2.4. Conclusions and recommendation**

to stop the release of jet-fuel in future.

Composition of contaminant (Jet fuel)

a localized area.

**2.5. Appendix A**

Based on the modeling, it is concluded that the jet-fuel plume has neither expanded nor moved considerably. It is less than 250 ft (76.2 m) beyond the storage tanks and is oriented northeast to southwest. The level of concentration found in the simulated monitoring wells is significant, but because groundwater is brackish and thus unlikely to be used, no harmful effects are expected. However, with the continuous process of leaking jet-fuel from the storage depot, the level of concentration is expected to increase over the period of time. Regionally, the jet-fuel expansion will have on prominent effect over a longer areal coverage and will be confined to

An interdisciplinary investigation of the processes controlling the fate and transport hydro‐ carbons in the subsurface is needed. Concentrations should be performed in wells within and down-gradient of the plume, as field data would help develop a stronger argument for the fate of jet fuel in groundwater. Periodic observations need to be carried out in wells to have good control on the changes of groundwater chemistry. Defective storage depots need to be mended

Knowledge of the geochemistry of a contaminated aquifer is important to understand the chemical and biological processes controlling the migration of hydrocarbon contaminants in the subsurface. Originally, the jet fuel (kerosene oil) is the name assigned to a material with a

It is important to maintain the existing quality of groundwater because once contamination occurs; it is sometime difficult or rather impossible to clean the aquifer. There is high proba‐ bility associated with certain landuses like agriculture, industrial/urban land and drainage wells for contaminating the groundwater. The early detection of such contamination can be executed through proper monitoring of the groundwater quality. Many oil & gas companies are disposing off their producing water into the deep Ranikot formation in Bhit oil-field area of southern Pakistan. The producing water contains Total dissolved solid (TDS) within range of 18,000 - 22,000 mg/l besides oil condensate. There were concerns that the producing water is affecting the fresh water aquifers belonging to the overlying Nari and Kirthar formations. This phenomenon has been studied by the utilization of groundwater contaminant transport model. Injection has been monitored at 2100 meters depth in the Pab sandstone formation. A three-dimensional contaminant transport model was developed to simulate and monitor the migration of disposal of hydrocarbon exploration produced water in Injection well at 2000 meters depth in the Upper Cretaceous Pab sandstone in the study area. Framework of regional stratigraphic and structural geology, landform characteristics, climate and hydrogeological setup were used to model the subsurface aquifer. The shallow and deep-seated characteristics of geological formations were obtained from electrical resistivity sounding surveys, geophys‐ ical well-logging information and available drilling data. The modeling process comprised of steady-state and transient simulations of the prolific groundwater system and, predictive simulation of contaminants transport after 1-, 10- and 30-year of injection. The contaminant transport was evaluated from the bottom of the injection well and its short and, long-term effects were determined on aquifer system lying in varying hydrogeological and geological conditions.

#### *3.1.1. Description of study area*

The study area of Bhit oil field is located about 43 km south of the Manchhar Lake within longitudes 67 <sup>o</sup> 25' - 67 <sup>o</sup> 48' E and latitudes 26 <sup>o</sup> 01' - 26 <sup>o</sup> 30' N in Dadu district in Sindh province of Pakistan (Figure 11). Manchhar Lake, one of the largest lake in Pakistan and in Asia, is formed in a depression in the western side of the Indus River in Sindh province. The total catchment area of the lake is about 97,125 km2 [27]. The surface area of the lake fluctuates with the seasons from as little as 350 km² to as much as 520 km² [28]. The lake is fed by two canals, the Aral and the Danister emerging from the river Indus in their eastern side. Due to less rainfall and contamination of surface water, the Manchhar Lake contains brackish water. The elevation ranges between 45 m at the Manchhar Lake and 163 m towards Bhit study area. On regional scale the area is a part of Gigantic Indus river basin composed of alluvium transported by the river and its tributaries. The main surface water sources are Naig Nai stream originating from Bhit Mountain range in the west, Dhanar Dhoro stream passing close to the oil-field plan, besides other small intermittent streams which remain dry in most parts of the year. The discharge data of these streams were not available. The water supply for local communities is maintained by the springs originating from the nearby Limestone Mountains of the Kirthar Range. According to reconnaissance studies conducted in the area, fresh groundwater source is available at different locations in the alluvium. Presently, there is no significant groundwater development in the area. Only few water wells were constructed to fulfill the requirements of local communities. The major sources of potable water to humans and livestock and for irrigation are unconsolidated aquifers. A water supply scheme consisting of four tube wells at Jhangara village is providing water to the villages between Manchhar Lake and Jhangara. The rainfall is scanty. Average annual rainfall is about 200 mm. It is higher in summer months like July and August due to prevalence of monsoon conditions. The aquifer area is located in the alluvial deposits along the Naig Nai stream.

**Figure 11.** Location of study area in southern Pakistan

The data of subsoil properties, aquifer characteristics and existing groundwater conditions were collected through reconnaissance level field investigations including geophysical survey. The electrical resistivity soundings survey (ERSS) was conducted primarily to collect required input data for the modeling. The surface drainage and topographic information were extracted from the topo-sheets prepared by Survey of Pakistan on 1:50,000 scale. The climate data i.e. precipitation and temperature, of 1961-1976 period were acquired of nearby meteorological stations like Karachi, Nawabshah, Moenjodaro and Khuzdar from Meteorological Department and Water and Power Development Authority (WAPDA) Pakistan. The published literatures

Groundwater and Contaminant Hydrology http://dx.doi.org/10.5772/54732 187

The geophysical/drill logs of the injection well field suggest that the subsurface material is composed of layers of sandstone, limestone, dolomite and clay-stone of different formations. One of the interpretative seismic sections of the Bhit concession area (shown in Figure 12) indicates the deeper Chiltan limestone formation beyond the depth of injection well. The hard rock aquifers are mainly composed of partially fractured limestone and sandstone belonging to Nari and Kirthar formations. Limestone, which is the dominant formation, has solution channels due to water action having secondary permeability characteristics. Further, chances

of the region i.e. see [29-32] had been used to firm up the study results.

**3.2. Material and methods**

*3.2.1. Geophysical data analysis*

The water table is generally in phreatic to semi-unconfined conditions. The observation wells drilled in the area indicated water table depth of about 12 m and hydraulic head value of 148.8 meter above sea level (masl). The groundwater flow is generally from southwest to‐ wards northeast direction. The flow direction of groundwater is true replica of the flow di‐ rection of Naig Nai stream draining the area. The groundwater level is mainly influenced by seasonal floods, stream flows and tubewell (water well) discharge. The fluctuations are small in the deep water table in the piedmont plain.

**Figure 11.** Location of study area in southern Pakistan

## **3.2. Material and methods**

stratigraphic and structural geology, landform characteristics, climate and hydrogeological setup were used to model the subsurface aquifer. The shallow and deep-seated characteristics of geological formations were obtained from electrical resistivity sounding surveys, geophys‐ ical well-logging information and available drilling data. The modeling process comprised of steady-state and transient simulations of the prolific groundwater system and, predictive simulation of contaminants transport after 1-, 10- and 30-year of injection. The contaminant transport was evaluated from the bottom of the injection well and its short and, long-term effects were determined on aquifer system lying in varying hydrogeological and geological

The study area of Bhit oil field is located about 43 km south of the Manchhar Lake within

of Pakistan (Figure 11). Manchhar Lake, one of the largest lake in Pakistan and in Asia, is formed in a depression in the western side of the Indus River in Sindh province. The total catchment area of the lake is about 97,125 km2 [27]. The surface area of the lake fluctuates with the seasons from as little as 350 km² to as much as 520 km² [28]. The lake is fed by two canals, the Aral and the Danister emerging from the river Indus in their eastern side. Due to less rainfall and contamination of surface water, the Manchhar Lake contains brackish water. The elevation ranges between 45 m at the Manchhar Lake and 163 m towards Bhit study area. On regional scale the area is a part of Gigantic Indus river basin composed of alluvium transported by the river and its tributaries. The main surface water sources are Naig Nai stream originating from Bhit Mountain range in the west, Dhanar Dhoro stream passing close to the oil-field plan, besides other small intermittent streams which remain dry in most parts of the year. The discharge data of these streams were not available. The water supply for local communities is maintained by the springs originating from the nearby Limestone Mountains of the Kirthar Range. According to reconnaissance studies conducted in the area, fresh groundwater source is available at different locations in the alluvium. Presently, there is no significant groundwater development in the area. Only few water wells were constructed to fulfill the requirements of local communities. The major sources of potable water to humans and livestock and for irrigation are unconsolidated aquifers. A water supply scheme consisting of four tube wells at Jhangara village is providing water to the villages between Manchhar Lake and Jhangara. The rainfall is scanty. Average annual rainfall is about 200 mm. It is higher in summer months like July and August due to prevalence of monsoon conditions. The aquifer area is located in the

The water table is generally in phreatic to semi-unconfined conditions. The observation wells drilled in the area indicated water table depth of about 12 m and hydraulic head value of 148.8 meter above sea level (masl). The groundwater flow is generally from southwest to‐ wards northeast direction. The flow direction of groundwater is true replica of the flow di‐ rection of Naig Nai stream draining the area. The groundwater level is mainly influenced by seasonal floods, stream flows and tubewell (water well) discharge. The fluctuations are

01' - 26 <sup>o</sup>

30' N in Dadu district in Sindh province

48' E and latitudes 26 <sup>o</sup>

186 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

conditions.

longitudes 67 <sup>o</sup>

*3.1.1. Description of study area*

25' - 67 <sup>o</sup>

alluvial deposits along the Naig Nai stream.

small in the deep water table in the piedmont plain.

The data of subsoil properties, aquifer characteristics and existing groundwater conditions were collected through reconnaissance level field investigations including geophysical survey. The electrical resistivity soundings survey (ERSS) was conducted primarily to collect required input data for the modeling. The surface drainage and topographic information were extracted from the topo-sheets prepared by Survey of Pakistan on 1:50,000 scale. The climate data i.e. precipitation and temperature, of 1961-1976 period were acquired of nearby meteorological stations like Karachi, Nawabshah, Moenjodaro and Khuzdar from Meteorological Department and Water and Power Development Authority (WAPDA) Pakistan. The published literatures of the region i.e. see [29-32] had been used to firm up the study results.

## *3.2.1. Geophysical data analysis*

The geophysical/drill logs of the injection well field suggest that the subsurface material is composed of layers of sandstone, limestone, dolomite and clay-stone of different formations. One of the interpretative seismic sections of the Bhit concession area (shown in Figure 12) indicates the deeper Chiltan limestone formation beyond the depth of injection well. The hard rock aquifers are mainly composed of partially fractured limestone and sandstone belonging to Nari and Kirthar formations. Limestone, which is the dominant formation, has solution channels due to water action having secondary permeability characteristics. Further, chances of transport contamination could take place much significantly through fracture zone of limestone and dolomitic formations. The hydraulic properties of the underlying overburden and rocks evaluated through geophysical well log data are shown in Table 5.

*3.2.2. Model conceptualization and simulation*

The groundwater flow system was treated as a multi-layered system. The upper layer aquifer is mainly unconfined and at depths where silt and clay horizons are present, the sand could probably cause partial confinement in some areas. The chances of transport contamination could take place much significantly through fracture zone of limestone and dolomitic forma‐ tions. Seven aquifer layers were defined on the basis of physical characteristics of lithological formations (Table 1). The disposal of produced water in the injection well is set at 2,100 m depth in the Pab sandstone formation. The main source of recharge to groundwater is rainfall which is highly variable. During rains, different nullahs carry flows and infiltrations through the piedmonts and alluvial fans and cause subsequent lateral movement at depth. The recharge is higher in summer especially in the months of July and August. It may occur during rainy month of March in winter, a mean value for the limestone recharge may be taken as 200 mm/ year based on the recharge data of rainfall. The main discharge components are groundwater extraction from water wells/dug wells, evapotranspiration, and spring discharge. The abstrac‐ tion of groundwater becomes higher during months of little recharge to groundwater i.e. November to January, which may affect the storage of groundwater for a limited period.

The MODFLOW code and the MT3DMS code were used to solve the flow and transport equations. The model domain comprised of 40 x 30 grid network with total area of 1200 sq km (Figure 13). First, the model was run for steady state condition. The model took up ground‐ water extraction from the oil-field water wells (TW1 & TW2) and the community tube wells

wastewater from the injection well was considered during simulation. Once the flow model was completed and run was carried out, the contaminant transport model was set and

The steady-state hydraulic heads were used as initial condition in MT3DMS option available in PM5 to simulate the dispersion of plume. The MT3DMS model simulates the processes i.e. advection, mechanical dispersion, retardation, decay and molecular diffusion related to the fate of contaminant. Initial concentration was set to zero in all the layers. In Advection package, 3rd order TVD scheme [12, 33] was selected. This method is considered as a good compromise between the standard finite difference and particle tracking approaches. In dis‐ persion package, TRPT (Horizontal transverse dispersivity/Longitudinal dispersivity) was set to 0.3 for all the layers except in layer 3, where it was set to 0.1. In dispersion package, TRPV (Vertical transverse dispersivity/Longitudinal dispersivity) was set to 0.3 for all the layers except in layer 3. For this layer, it was set to 0.1. Longitudinal dispersivity was set to 10 m. There was no sorption selected in chemical reaction package. The injection well was set at layer 7 in the sink/source menu. The concentration of the injection well liquid was con‐ sidered to be 100 ppm. The model was then simulated for 1-, 10-, and 30-year period for

simulated to evaluate the groundwater contamination and movement of plume.

studying the behaviors produced wastewater in and around the injection well.

/sec and discharge rate of

/sec. The injection of the produced

Groundwater and Contaminant Hydrology http://dx.doi.org/10.5772/54732 189

(TW3 & TW4). The discharge rate of the water wells was 0.0083 m3

deep seated injection well (DW) as supplied was 0.00152 m3

*3.2.3. Simulation of contaminant transport of injection well*


The hydraulic conductivity of the unconsolidated deposits is about 19 m/day and the effective porosity is 0.25. Several hydrocarbon wells are producing gas in the Bhit concession area. The Ghazij (claystone) was found to be a cap rock (regional seal) over Pab sandstone - an enriched reservoir of hydrocarbon. The Ranikot formation is a prolific rock unit having good transmis‐ sivity and storativity to accumulate disposal waste of the producing gas Bhit concession.

**Figure 12.** Interpretative seismic depth section indicating several geological formations

## *3.2.2. Model conceptualization and simulation*

of transport contamination could take place much significantly through fracture zone of limestone and dolomitic formations. The hydraulic properties of the underlying overburden

> *K (m/d)*

113 5 130

616 24 2000

1,259 0.165 4

1,799 0.68 14.97

2,000 to onward 0.138 17.5

*Transmissivity T (m2/d)*

(S=.01)

(S=.04)

(S=0.00005)

(S=0.000007)

(S=0.000005)

and rocks evaluated through geophysical well log data are shown in Table 5.

*(m)*

3 Ghazij (claystone) 762 zero Regional Seal

The hydraulic conductivity of the unconsolidated deposits is about 19 m/day and the effective porosity is 0.25. Several hydrocarbon wells are producing gas in the Bhit concession area. The Ghazij (claystone) was found to be a cap rock (regional seal) over Pab sandstone - an enriched reservoir of hydrocarbon. The Ranikot formation is a prolific rock unit having good transmis‐ sivity and storativity to accumulate disposal waste of the producing gas Bhit concession.

5 Dunghan (dolomite and claystone) - - -

**Figure 12.** Interpretative seismic depth section indicating several geological formations

*Formation Depth*

188 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

Nari (sandstone)

Kirthar (limestone)

Laki (limestone & dolomite)

Ranikot (Lakhra+Bara+Khadro)

Pab (sandstone)

**Table 5.** Summary of the Aquifer characteristics of hard rock formations

*S.No.*

1

2

4

6

7

The groundwater flow system was treated as a multi-layered system. The upper layer aquifer is mainly unconfined and at depths where silt and clay horizons are present, the sand could probably cause partial confinement in some areas. The chances of transport contamination could take place much significantly through fracture zone of limestone and dolomitic forma‐ tions. Seven aquifer layers were defined on the basis of physical characteristics of lithological formations (Table 1). The disposal of produced water in the injection well is set at 2,100 m depth in the Pab sandstone formation. The main source of recharge to groundwater is rainfall which is highly variable. During rains, different nullahs carry flows and infiltrations through the piedmonts and alluvial fans and cause subsequent lateral movement at depth. The recharge is higher in summer especially in the months of July and August. It may occur during rainy month of March in winter, a mean value for the limestone recharge may be taken as 200 mm/ year based on the recharge data of rainfall. The main discharge components are groundwater extraction from water wells/dug wells, evapotranspiration, and spring discharge. The abstrac‐ tion of groundwater becomes higher during months of little recharge to groundwater i.e. November to January, which may affect the storage of groundwater for a limited period.

The MODFLOW code and the MT3DMS code were used to solve the flow and transport equations. The model domain comprised of 40 x 30 grid network with total area of 1200 sq km (Figure 13). First, the model was run for steady state condition. The model took up ground‐ water extraction from the oil-field water wells (TW1 & TW2) and the community tube wells (TW3 & TW4). The discharge rate of the water wells was 0.0083 m3 /sec and discharge rate of deep seated injection well (DW) as supplied was 0.00152 m3 /sec. The injection of the produced wastewater from the injection well was considered during simulation. Once the flow model was completed and run was carried out, the contaminant transport model was set and simulated to evaluate the groundwater contamination and movement of plume.

## *3.2.3. Simulation of contaminant transport of injection well*

The steady-state hydraulic heads were used as initial condition in MT3DMS option available in PM5 to simulate the dispersion of plume. The MT3DMS model simulates the processes i.e. advection, mechanical dispersion, retardation, decay and molecular diffusion related to the fate of contaminant. Initial concentration was set to zero in all the layers. In Advection package, 3rd order TVD scheme [12, 33] was selected. This method is considered as a good compromise between the standard finite difference and particle tracking approaches. In dis‐ persion package, TRPT (Horizontal transverse dispersivity/Longitudinal dispersivity) was set to 0.3 for all the layers except in layer 3, where it was set to 0.1. In dispersion package, TRPV (Vertical transverse dispersivity/Longitudinal dispersivity) was set to 0.3 for all the layers except in layer 3. For this layer, it was set to 0.1. Longitudinal dispersivity was set to 10 m. There was no sorption selected in chemical reaction package. The injection well was set at layer 7 in the sink/source menu. The concentration of the injection well liquid was con‐ sidered to be 100 ppm. The model was then simulated for 1-, 10-, and 30-year period for studying the behaviors produced wastewater in and around the injection well.

Manchhar Lake. The groundwater flow had shown decline in confining layers like 4 and 7.

community wells were used for the study. The results obtained from the 3-D transport model are shown in Table 6 and transport of contaminant plumes in three simulation periods in

**Layer 1 Year 10 Years 30 Years**

1 - 2.21 x 10-23 1.96 x 10-20

2 3.16 x 10-27 3.88 x 10-21 1.32 x 10-18

3 3.14 x 10-21 3.54 x 10-16 4.01 x 10-14

4 3.16 x 10-15 3.18 x 10-11 1.06 x 10-09

5 9.58 x 10-10 9.53 x 10-07 2.45 x 10-05

6 8.04 x 10-05 7.95 x 10-03 7.21 x 10-02

After 30 years of simulation period, only traces of contamination were found in Ghazij Formation. Moreover, it is found that after 1-year period of simulation the produced wastewater will reach upward in layer-5 (Ranikot Formation) emerging from layer 7 (Pab sandstone) as shown in Figure 14. In this period, no contamination was found in layer 1 and 2. In 10-year simulation a plume of produced water moved from layer 7 to layer 5 (Figure 15). Only traces of contamination were found in layer 3 (Ghazi For‐ mation). In Figure 6, plume of produced water contamination indicates movement from layer 7 to layer 4 after 30 years simulation. The layer 3 was found to be acting as a regional confining seal. In this layer, only traces of contamination were present. The movement of produced wastewater was found within a radius of 3 km at the bot‐ tom of injection well in the Pab sandstone. The upper aquifers in the alluvial deposit, Nari sandstone, and Kirthar limestone was remain safe from the effects of produced wastewater disposal from the deep seated injection well. The community wells tapping in the upper few tens of meters, naturally oozing springs and the Manchhar Lake lo‐ cated about 43 km from the injection well were also found to be safe from the effects of produced water injection even after contaminant transport simulation of 30-year pe‐ riod. The development of plume was significant in layer 7 and upward in the three

**Table 6.** Maximum concentration observed in different simulation periods (ppm)

6.28 x 10-02 0.626 1.861

/sec along with

191

Groundwater and Contaminant Hydrology http://dx.doi.org/10.5772/54732

The two tube wells of the oil-fields each discharging at the rate of 8.3 x 10-3 m3

Figures 14-16.

7 Pab (Sandstone)

cases (shown in Figures 14-16).

**Figure 13.** Manchhar Lake treated as constant head boundary; Bhit Range as impervious boundary. *DW* deep injec‐ tion well, *TW* tube well

#### **3.3. Results and discussion**

The computations were carried out for three cases i.e. in case-I, Injection well was continuously discharged for one-year period, in case-II it well was simulated for 10-year period and in case-III for 30-year period. The hydraulic heads and drawdown were computed in all three cases. The velocity vectors prominent in layer-1 tend to move in the northeast direction towards Manchhar Lake. The groundwater flow had shown decline in confining layers like 4 and 7. The two tube wells of the oil-fields each discharging at the rate of 8.3 x 10-3 m3 /sec along with community wells were used for the study. The results obtained from the 3-D transport model are shown in Table 6 and transport of contaminant plumes in three simulation periods in Figures 14-16.


**Table 6.** Maximum concentration observed in different simulation periods (ppm)

**Figure 13.** Manchhar Lake treated as constant head boundary; Bhit Range as impervious boundary. *DW* deep injec‐

190 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

The computations were carried out for three cases i.e. in case-I, Injection well was continuously discharged for one-year period, in case-II it well was simulated for 10-year period and in case-III for 30-year period. The hydraulic heads and drawdown were computed in all three cases. The velocity vectors prominent in layer-1 tend to move in the northeast direction towards

tion well, *TW* tube well

**3.3. Results and discussion**

After 30 years of simulation period, only traces of contamination were found in Ghazij Formation. Moreover, it is found that after 1-year period of simulation the produced wastewater will reach upward in layer-5 (Ranikot Formation) emerging from layer 7 (Pab sandstone) as shown in Figure 14. In this period, no contamination was found in layer 1 and 2. In 10-year simulation a plume of produced water moved from layer 7 to layer 5 (Figure 15). Only traces of contamination were found in layer 3 (Ghazi For‐ mation). In Figure 6, plume of produced water contamination indicates movement from layer 7 to layer 4 after 30 years simulation. The layer 3 was found to be acting as a regional confining seal. In this layer, only traces of contamination were present. The movement of produced wastewater was found within a radius of 3 km at the bot‐ tom of injection well in the Pab sandstone. The upper aquifers in the alluvial deposit, Nari sandstone, and Kirthar limestone was remain safe from the effects of produced wastewater disposal from the deep seated injection well. The community wells tapping in the upper few tens of meters, naturally oozing springs and the Manchhar Lake lo‐ cated about 43 km from the injection well were also found to be safe from the effects of produced water injection even after contaminant transport simulation of 30-year pe‐ riod. The development of plume was significant in layer 7 and upward in the three cases (shown in Figures 14-16).

(a) (b) (c)

Groundwater and Contaminant Hydrology http://dx.doi.org/10.5772/54732 193

**Figure 16.** Upward movement of the plume from layer 7 to layer 6 (a); from layer 6 to layer 5 (b) and from layer 5 to

The results of the three-dimensional groundwater modeling study highlighted the hydrogeo‐ logical characteristics and features of the contaminant transport of the deep injected tube well in the Bhit oil-field area. The groundwater contaminant transport modeling technique has proved to be effective in simulating produced wastewater plume from the deep seated injection well. The study would provide base for evaluating risks of contaminants on long term basis in similar conditions in future. Risk of expansion of plume regionally does not exist as the disposal of wastewater is made in the deeper horizon well below the aquifers and also the

Thorough understanding of surface hydrology, hydrogeological conditions and contaminant behavior in the aquifer system coupled with application of reliable modeling techniques could be helpful in dealing with water management issues related to contaminant hydrology.

, Gulraiz Akhter1

3 College of Earth and Environmental Sciences, Punjab University, Lahore, Pakistan

1 Department of Earth Sciences, Quaid-i-Azam University, Islamabad, Pakistan

and Iftikhar Ahmad3

layer 4 (c) after 30 year simulation period

**4. Conclusions**

quantity is quite limited.

**Author details**

Zulfiqar Ahmad1

, Arshad Ashraf2

2 National Agricultural Research Center,Islamabad, Pakistan

**Figure 14.** Upward movement of the plume from layer 7 to layer 6 (a); from layer 6 to layer 5 (b) and from layer 5 to layer 4 (c) after 1 year simulation period

**Figure 15.** Upward movement of the plume from layer 7 to layer 6 (a); from layer 6 to layer 5 (b) and from layer 5 to layer 4 (c) after 10 year simulation period

**Figure 16.** Upward movement of the plume from layer 7 to layer 6 (a); from layer 6 to layer 5 (b) and from layer 5 to layer 4 (c) after 30 year simulation period

## **4. Conclusions**

(a) (b) (c)

**Figure 14.** Upward movement of the plume from layer 7 to layer 6 (a); from layer 6 to layer 5 (b) and from layer 5 to

(a) (b) (c)

**Figure 15.** Upward movement of the plume from layer 7 to layer 6 (a); from layer 6 to layer 5 (b) and from layer 5 to

layer 4 (c) after 1 year simulation period

192 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

layer 4 (c) after 10 year simulation period

The results of the three-dimensional groundwater modeling study highlighted the hydrogeo‐ logical characteristics and features of the contaminant transport of the deep injected tube well in the Bhit oil-field area. The groundwater contaminant transport modeling technique has proved to be effective in simulating produced wastewater plume from the deep seated injection well. The study would provide base for evaluating risks of contaminants on long term basis in similar conditions in future. Risk of expansion of plume regionally does not exist as the disposal of wastewater is made in the deeper horizon well below the aquifers and also the quantity is quite limited.

Thorough understanding of surface hydrology, hydrogeological conditions and contaminant behavior in the aquifer system coupled with application of reliable modeling techniques could be helpful in dealing with water management issues related to contaminant hydrology.

## **Author details**

Zulfiqar Ahmad1 , Arshad Ashraf2 , Gulraiz Akhter1 and Iftikhar Ahmad3


## **References**

[1] Boulding, J. R, & Ginn, J. S. Practical handbook of soil, vadose zone and groundwater contamination: assessment, prevention and remediation, CRC Press; (2004).

[16] EPADetermining Soil Response Action Levels Based on Potential Contaminant Migration to Groundwater: A Compendium of Examples, EPA/540/(1989). , 2-89. [17] ASTMStandard Guide for Risk-Based Corrective Action Applied at Petroleum Release Sites, American Society for Testing and Materials, ASTM E-Philadelphia, PA; (1995). ,

Groundwater and Contaminant Hydrology http://dx.doi.org/10.5772/54732 195

[18] Newell, C. J, Mcleod, R. K, & Gonzales, J. R. BIOSCREEN: Natural Attenuation Decision Support System, User's Manual Version 1.3, EPA/600/R-96/087, National Risk Man‐ agement Research Laboratory, Office of Research and Development, U. S. Environ‐

[19] Chiang, W-H, & Kinzelbach, W. Processing MODFLOW. A simulation system for

[20] Pollock, D. W. Documentation of a computer program to compute and display path lines using results from the U.S Geological Survey modular three-dimensional finitedifference groundwater flow model: U.S Geological Survey, open-file report Denver;

[21] Zhou, Y. Z. Sampling frequency for monitoring the actual state of groundwater

[22] Plus, R. W, & Paul, C. J. Multi-layer sampling in conventional monitoring wells for improved estimation of vertical contaminant distribution and mass. *Jour of Contam.*

[23] Fisher, R. S, & Goodmann, P. T. Characterizing groundwater quality in Kentucky: from site selection to published information. Proceedings of 2002 National Monitoring Conference, national water Quality Monitoring council, May Madison, Wisconsin;

[24] Zeru, A. Investigations numériques sur l'inversion des courbes de concentration issues d'un pompage pour la quantification de la pollution de l'eau souterraine / Numerical investigations on the inversion of pumped concentrations for groundwater pollution

[25] Bauer, S, Bayer-raich, M, Holder, T, Kolesar, C, Muller, D, & Ptak, T. Quantification of groundwater contamination in an urban area using integral pumping tests *Jour of*

[26] Bockelmann, A, Zamfirescu, Z, Ptak, T, Grathwohl, P, & Teutsch, G. Quantification of mass fluxes and natural attenuation rates at an industrial site with a limiten monitoring

[27] NESPAKMaster Feasibility studies for flood management of hill-torrents of Pakistan, Supporting Vol-IV Sindh Province. National Engg. Services Pak. (Pvt) ltd. Lahore,

network: a case study site. *Jour of Contam. Hydrology* (2003). , 60, 97-121.

[28] http://enwikipedia.org/wiki/Manchar\_Lake (accessed 10 September (2012).

quantification. PhD Thesis, Université Louis Pasteur (France); (2004). , 192.

mental Protection Agency, Cincinnati, 45268, Ohio; (1996).

modeling groundwater flow and pollution; (2001).

systems, *Jour of Hydrology* (1996). , 180, 301-318.

1739-95.

(1989). , 89-381.

(2002). , 19-23.

Pakistan; (1998).

*Hydrology* (1997). , 25, 85-111.

*Contam. Hydrology* (2004). , 75, 183-213.


[16] EPADetermining Soil Response Action Levels Based on Potential Contaminant Migration to Groundwater: A Compendium of Examples, EPA/540/(1989). , 2-89.

**References**

03-4053.

A1, Washington D.C.; (1988). , 83-875.

194 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

(1999). , 25, 457-462.

*Research* (2000). , 34(3)

US Army Corps of Engineers, Washington, DC: (1999).

depots/ Installations/ airfields, Shell Pakistan Limited; (2000).

Indus Basin. *Geophysical Journal International* (2008). , 173, 17-24.

Facilities, NAVFAC Pacific, Oahu, Hawaii; (2007).

Contaminant Transport Model, *Groundwater* (1986). , 24(1)

Journal of Environmental Science and Health (2008). , 43(6)

Inc. for the U. S. Environmental Protection Agency; (1997).

[14] Fetter, C. W. Applied hydrogeology (4th edition), Prentice Hall; (2000).

[1] Boulding, J. R, & Ginn, J. S. Practical handbook of soil, vadose zone and groundwater contamination: assessment, prevention and remediation, CRC Press; (2004).

[2] Walter, D. A, & Masterson, J. P. Simulation of Advective Flow under Steady-State and Transient Recharge Conditions, Camp Edwards, Massachusetts Military Reservation, Cape Cod, Massachusetts, Water-Resources Investigations Report USGS; (2003). ,

[3] Anderson, M. P, & Woessner, W. W. Applied Groundwater Modeling, Simulation of Flow and Advective Transport. Academic Press, Inc San Diego, California; (1992). [4] Mcdonald, M. G, & Harbaugh, A. W. ModFlow, A modular three-dimensional finitedifference ground-water flow model, U.S. Geological Survey Open-File Report Chapter

[5] Zheng, C, & Wang, P. P. A Modular Three-Dimensional Multispecies Transport Model,

[6] Poeter, E. P, & Hill, M. C. UCODE, a computer modelling. *Computers & Geosciences*

[7] Kim, J, & Corapcioglu, M. Y. Modeling dissolution and volatilization of LNAPL sources migrating on the groundwater table, *Journal of Contaminant Hydrology* (2003). , 65(1-2)

[8] Mott McDonald Pakistan (pvt) LtdA report on soil and groundwater Assessment of

[9] Ashraf, A, & Ahmad, Z. Regional groundwater flow modeling of upper Chaj Doab,

[10] El-Kadi, A. Groundwater Modeling Services for Risk Assessment Red Hill Fuel Storage

[11] Periago, E. L, Delgado, A. N, & Diaz-fierros, F. F. Groundwater contamination due to cattle slurry: modeling infiltration on the basis of soil column experiments, *Water*

[12] Eric, W, Strecker, E. W, & Chu, W. Parameter Identification of a Ground-Water

[13] Jin, S, Fallgren, P, Cooper, J, Morris, J, & Urynowicz, M. Assessment of diesel contam‐ ination in groundwater using electromagnetic induction geophysical techniques,

[15] KeyGroundwater Fate and Transport Evaluation Report: South Cavalcade Superfund Site, Houston, Texas, prepared by KEY Environmental, Inc. on behalf of Beazer East,


[29] Ahmad, Z, Ahmad, I, & Akhtar, G. A report on groundwater reserve estimation of the Ahmad Khan well field and its safe utilization for the 4000 TPD Luck Cement plant Pezu, D.I.Khan; (1994).

**Section 3**

**Water Resources Sustainability**


**Water Resources Sustainability**

[29] Ahmad, Z, Ahmad, I, & Akhtar, G. A report on groundwater reserve estimation of the Ahmad Khan well field and its safe utilization for the 4000 TPD Luck Cement plant

[31] WAPDALower Indus Report. Physical Resources, groundwater, supl.6.1.6, Tube wells

[33] Harten, A. High resolution schemes for hyperbolic conservation laws. *Jour of Comput.*

[30] Ahmad, N. Groundwater resources of Pakistan, Ripon printing press; (1974).

[32] WAPDARohri Hydropower Project. Feasibility study, HEPO (1989). (97)

Pezu, D.I.Khan; (1994).

*Phys.* (1983). , 49, 357-393.

and Boreholes, Nara Command; (1965).

196 Current Perspectives in Contaminant Hydrology and Water Resources Sustainability

**Chapter 8**

**Geospatial Analysis of Water Resources for**

Matjaž Glavan, Rozalija Cvejić, Matjaž Tratnik and

Additional information is available at the end of the chapter

Marina Pintar

**1. Introduction**

http://dx.doi.org/10.5772/53528

a priority of any government in the world.

on fine temporal and spatial scale [4].

**Sustainable Agricultural Water Use in Slovenia**

Global population growth has greatly increased food demand. This, in turn, has intensified agricultural production, already the biggest consumer of water in the world [1]. Development of irrigation techniques has contributed to the global food production [2]. However, climate change simulations predict repeated droughts and deteriorating crop production, illustrating the critical need for sustainable irrigation [3]. Thus, a proactive water management strategy is

Globally, only 10% of estimated blue water (surface water, groundwater, and surface runoff) and 30% of estimated green water (evapotranspiration, soil water) resources are used for consumption. Nevertheless, water scarcity is a problem due to high variability of water resources availability in time and space [4]. Model results suggest that severe water scarcity occurs at least one month per year in almost one half of the world river basins [5]. One third of the water volume currently supplied to irrigated areas is supplied by locally stored runoff [6]. It is estimated that small reservoirs construction could increase global cereal production in low-yield regions (i.e. Africa, Asia) by approximately 35% [6]. Global water scarcity problems can now be, due to advances in hydrology science in the last decades, easily assessed

Irrigation development and management in Slovenia have completely stagnated in the last decade due to financial shortages. In 1994 the Slovenian government adopted a strategy for agricultural land irrigation (i.e. National Irrigation Programme) [7]. In 1999, the World Bank

> © 2013 Glavan et al.; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use,

© 2013 Glavan et al.; licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

distribution, and reproduction in any medium, provided the original work is properly cited.
