Analysis of aquifer mineralization by paleodrainage channels

Analysis of aquifer mineralization by paleodrainage channels

Journal of Hydrology 277 (2003) 280–304 www.elsevier.com/locate/jhydrol Analysis of aquifer mineralization by paleodrainage channels Hillel Rubina,*,...

833KB Sizes 0 Downloads 41 Views

Journal of Hydrology 277 (2003) 280–304 www.elsevier.com/locate/jhydrol

Analysis of aquifer mineralization by paleodrainage channels Hillel Rubina,*, Robert W. Buddemeierb a

Faculty of Civil Engineering, Technion—Israel Institute of Technology, Haifa 32000, Israel b Kansas Geological Survey, The University of Kansas, Lawrence, KS 66047, USA Received 13 May 2002; accepted 21 March 2003

Abstract Mineralization of groundwater resources is a problem in south-central Kansas, due to the penetration of saline water from Permian bedrock formations into the overlying alluvial aquifer. One of the mechanisms involved in the mineralization involves small bedrock features of high permeability located in places occupied by streams and rivers in past geological eras. These geological features are termed ‘paleodrainage channels’. The permeability of the overlying aquifer can be significantly smaller than that of the channel fill material. The comparatively fast migration of saline water through these channels of high permeability is associated with the transfer of minerals into the overlying freshwater aquifer. This study applies a set of boundary layer approaches to quantify the process of mineral transfer from the channels into the aquifer. The methods used in the present study provide quick estimation and evaluation of the dilution of the channel flow, as well as mineral concentration profile changes in the mineralized zone created in the overlying aquifer. More generally, the method can also be useful for the analysis and evaluation of various types of groundwater contamination in heterogeneous aquifers. The application of the method is exemplified by a complete set of calculations characterizing the possible mineralization process at a specific channel in south central Kansas. Sensitivity analyses are performed and provide information about the importance of the various parameters that affect the mineralization process. Some possible scenarios for the aquifer mineralization phenomena are described and evaluated. It is shown that the channel mineralization may create either several stream tubes of the aquifer with high mineral concentration, or many stream tubes mineralized to a lesser extent. Characteristics of these two patterns of aquifer mineralization are quantified and discussed. q 2003 Published by Elsevier Science B.V. Keywords: Aquifer mineralization; Paleodrainage channels; Boundary layers; Great Bend Prairie aquifer; Arkansas River; Groundwater salinity

1. Introduction This paper contributes to the understanding of major mechanisms and phenomena involved in * Corresponding author. Address: Department of Civil and Environmental Engineering, Stanford University, Stanford, CA 94305 4020, USA. Tel.: þ 972-4-829-2306; fax: þ 972-4-822-8898. E-mail address: [email protected], [email protected] (H. Rubin).

the mineralization of groundwater. Groundwater in the Great Bend Prairie alluvial aquifer of southcentral Kansas is subject to contamination by the intrusion of highly mineralized water originating from dissolution of minerals, mainly halite, in the underlying bedrock. Although the process in this case results in saline water, we use the term ‘mineralization’ in preference to ‘salinization’ to distinguish it from the very different processes of

0022-1694/03/$ - see front matter q 2003 Published by Elsevier Science B.V. doi:10.1016/S0022-1694(03)00123-9

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

Nomenclature a Ac

transverse dispersivity [L] cross section area of the paleodrainage channel [L2] b width of the channel [L] B thickness of the aquifer [L] BL boundary layer cr ratio between Ca and Cb at the boundary between the inner and outer BLs Ca mineral concentration of the overlying aquifer [ML23] Cav average mineral concentration of the overlying aquifer flowline [ML23] Cc mineral concentration of the channel flow [ML23] Cc0 mineral concentration of reference of the channel flow, mineral concentration at the channel entrance [ML23] Cb mineral concentration at the bottom of the aquifer [ML23] Ct mineral concentration at the top of the aquifer [ML23] CT mineral concentration at the top of the ROI, acceptable value of the mineral concentration [ML23] F function defined in Eq. (34) J number of channel segments comprising the dilution length Ka hydraulic conductivity of the aquifer [LT21] Kc hydraulic conductivity of the channel [LT21] Lc distance from the channel entrance to a cross section associated with the production of a particular value of Cav LD dilution length of the channel [L] LE establishment section length [L] Lj length of a channel segment [L] LR restructuring section length [L] m number of the term in a series expansion n; n1 ; n2 ; n3 ; n4 ; n5 ; n6 power coefficients of BL series expansions qa specific discharge of the aquifer flow [LT21]

281

flow-rate of the channel [L3T21] region of interest time time interval required for the advection of fluid particle from the downstream edge of the mineral penetration site to the attachment point TSBL top specified boundary layer Va flow velocity of the aquifer flow [LT21] Vc flow velocity of the channel flow [LT21] VR the velocity ratio x longitudinal coordinate [L] x0 coordinate extended along the channel centerline [L] X longitudinal coordinate of the flowline extended downstream of the attachment point [L] y transverse coordinate [L] z vertical coordinate [L] a1 coefficient defined in Eq. (15) a2 coefficient defined in Eq. (16) a3 coefficient defined in Eq. (31) bi ði ¼ 1; …; 5Þ coefficient defined in Appendix A g coefficient defined in Eq. (26) d thickness of the ROI [L] du thickness of the inner BL [L] du0 thickness of the inner BL at the right hand side of the attachment point [L] d0 thickness of the mineralized zone [L] Dx0 interval along the channel centerline [L] z outer BL coordinate zT value of z at the top of the ROI h BL coordinate hT value of h at the top of the ROI u angle of orientation of the channel l inner BL coordinate lT value of l at the top of the ROI j local coordinate of the interface of direct contact [L] fa porosity of the aquifer fc porosity of the channel v BL coordinate used at the mineral penetration site of a winding channel

Qc ROI t T

282

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

agricultural salinization and saltwater intrusion into coastal aquifers. The region of major groundwater mineralization considered is shown in Fig. 1. There have been numerous studies on the geology and hydrogeology of the region (e.g. Latta, 1950; Layton and Berry, 1973; Fader and Stullken, 1978; Cobb, 1980). Fig. 2 shows a cross section through a typical part of the salt-affected area of the aquifer, approximately 50 km north of the case study site used for the example calculations in this paper. Thickness of the alluvial aquifer ranges from about 25 to 60 m, and is about 40 m at the case study site.

During the past decade, significant efforts have been invested in measuring mineral concentration distribution in the region (Young, 1992; Whittemore, 1993; Buddemeier et al., 1992, 1994; Young and Rubin, 1998). These measurements have been used for the identification of sources of minerals (Whittemore, 1993), and of the mechanisms involved in the mineralization processes (Buddemeier et al., 1994). This paper focuses on mineralization mechanisms potentially associated with buried paleodrainage channels incised into the surface of the Cretaceous- and Permian-age bedrock formations (Fig. 1b).

Fig. 1. The region of major groundwater mineralization in south central Kansas. The irregular N–S line near Hwy 281 is the eastern limit of the Cretaceous bedrock confining layer beneath the alluvial Great Bend Prairie aquifer. (A) East of the confining layer boundary, groundwater salinity is a problem; contaminated and vulnerable areas are shown by shading. (B) Inferred locations of bedrock paleodrainage channels are depicted (from Sophocleous, 1991). The location of the sample calculation site used for this paper (Fig. 9) is shown in both A and B.

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

283

Fig. 2. Cross section illustration of the general hydrogeology of the Great Bend Prairie aquifer.

Sophocleous (1991) reviewed the geology of the region, and inferred that in the Great Bend Prairie aquifer a number of highly transmissive bedrock topographic features, termed ‘paleodrainage channels’, exist at the aquifer – bedrock interface (Fig. 1b).

By following changes of water level in a grid of monitoring wells, Sophocleous (1991) obtained substantial support for the hypothesized existence and effects of paleodrainage channels in the region. He also demonstrated that the channel sediments are

284

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

permeable, capable of rapid pressure pulse transmission, and may serve as water transmission conduits, and noted that the channels tend to have better quality water than the surrounding aquifer when they have good connection to a surface stream. Just as the bedrock channels can act as conduits of fresh water, they may act to transmit saline water through the aquifer when they are better connected with the bedrock than with surface water sources. The channels are incised into the bedrock formation, which consists of gently sloping near-horizontal layers of sandstone (Fig. 3), in which the saltwater (like the overlying freshwater) flows down-gradient in a generally easterly direction. Thus, permeable layers of sandstone exposed at the head and along the walls of the channels are sources of bedrock discharge into the alluvium contained in the channel. An additional source may be high-density brine discharged at upgradient locations that flows over the surface of the bedrock at the bedrock – alluvium interface. When this flow encounters a channel it will tend to flow into it with some sinking, thus making the channel a collection sink for brine discharged elsewhere as well as for local discharge. Statistical analysis of groundwater quality data shows a strong correlation between salinity and depth (Young et al., 1998), and a test drilling and monitoring program conducted by Groundwater Management District No. 5 has confirmed that bedrock lows and channels contain water of substantially higher salinity than the shallower aquifer (D. Zehr, pers. comm.). The saline water that enters the channel may flow along the channel course with a velocity an order of magnitude or more higher than that of the groundwater flowing through the overlying aquifer, based on the permeability estimates of Sophocleous (1991). This channel water of comparatively high mineral concentration and flow velocity will transfer minerals into the groundwater in the overlying aquifer groundwater, due to the direct contact between the two types of groundwater. The channel therefore acts as a secondary source of bedrock minerals by transferring them into the aquifer at locations distant from their primary discharge out of the bedrock. Thus, the detailed cross section, in terms of aquifer characteristics and the interactions of paleochannels with both bedrock and surface channels will be

important to assessing the mineralization effects in specific cases. The authors (Rubin and Buddemeier, 1996, 1998a – d) have developed conceptual and modeling approaches to provide quantitative estimates of phenomena and mechanisms leading to mineralization of the Great Bend Prairie aquifer. Mineral transfer into the freshwater aquifer may occur by either direct contact between the fresh and saline water, or by infiltration of saline water from the deep formation through semi-permeable strata into the freshwater aquifer. In most places, clay and shale layers effectively separate the freshwater aquifer from the saline water of the deep formations, so mineral penetrations are probably local phenomena. However, once the saline water penetrates the freshwater aquifer, the minerals are advected and disperse, and thereby contaminate the freshwater resources of the region. Rubin and Buddemeier (1996) suggested applying a boundary layer (BL) approach called the top specified boundary layer (TSBL) to analyze a variety of contaminant hydrology issues and problems associated with the mineral intrusion phenomenon that could not be analyzed by the traditional BL approaches. The major idea of the TSBL approach is to separate the definition of the region of similar normalized contaminant concentration profiles from the definition of the region of interest (ROI), in which contaminant concentration exceeds its acceptable value. Previous studies (Rubin and Buddemeier, 1998a – d), apply the TSBL approach to calculate and evaluate mineralization processes in an inland aquifer, either as a result of contact between saline and freshwater layers, or saltwater seepage into the freshwater aquifer. The present study applies the BL approach of the TSBL to a simplified quantitative analysis of mineral transport from the paleodrainage channels into and through the overlying aquifer, and suggests a possible set of mechanisms involved in the aquifer mineralization process as it is affected by the preferential transmission of highly mineralized water in bedrock channel features.

2. Basic conceptual approach Fig. 3 presents a schematic depiction of a cross section perpendicular to the channel centerline. It

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

285

† Curved channels, which cross a particular flowline once; and † Winding channels, which cross some flowlines more than once.

Fig. 3. Schematic cross section of the aquifer and the paleodrainage channels, illustrating the dimensions and hydraulic conductivities used in this study.

provides some information about dimensions and properties of the channel and overlying aquifer assumed for this analysis. The orientation of the channel centerline is commonly different from the direction of the flowlines in the overlying aquifer, and the orientation angle u between the aquifer flowlines and the channel centerline varies along the channel. As indicated by Fig. 4, the present study refers to two types of channels:

Fig. 4. Types of channels considered in the present study, and their relationship to the aquifer flow. (a) Curved channel. (b) Winding channel.

Fig. 5 shows the scheme for analyzing salt and water fluxes in an elementary volume of the channel. We calculate mineral concentration transport at the mineral penetration site and downstream of that location along flowlines of the overlying aquifer. Owing to the higher mineral concentration of the channel water, diffusion and transverse dispersion transfer minerals from the channel into the overlying aquifer water along the interface of contact between these two types of water, which is the mineral penetration site considered. We consider that the channel flow is completely mixed, and mineral concentration distribution in the channel is uniform. Mineral transfer in the overlying aquifer is affected by the build-up of the vertical mineral concentration profile. Effects of dispersion along the flowlines of the aquifer and in the horizontal transverse direction are usually much smaller, and therefore they are neglected. Due to the relatively small width of the channel and the assumed high permeability of its contents, we assume that the mineral concentration profile obtains its steady state shape very quickly along the mineral penetration site. We also neglect the minor flow component perpendicular to the channel centerline. This flow component originates from the difference between the hydraulic gradient direction and the direction of the channel centerline. Mass balance is a fundamental basis for the calculation approach. Fig. 5 shows how the balance between the different mineral fluxes is calculated. This figure incorporates the concept of a control section of the channel. There must be mass balance between the quantity of minerals entering that control section, and net fluxes of minerals leaving the control section by advection through the channel flow (flow through the high permeability portion of paleodrainage channel) and diffusion through the top of the control section into the major portion of the aquifer flow (flow through layers of lower permeability of the overlying aquifer). The mass balance described in preceding phrases is represented by the following

286

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

Fig. 5. An elementary volume of the channel with fluxes of channel water, groundwater, and mineral concentration.

expression:

›Cc 0 ›C dx þ Qc 0c dx0 ›t ›x   › ðB ¼ 2qa b cosu 0 Ca dz dx0 ›x 0

Ac f c

ð1Þ

where Ac is the cross section area of the channel; Cc is the mineral concentration of the channel flow; t is the time; Qc is the channel flow rate; Ca is the mineral concentration of the aquifer at a point with vertical distance z above the direct contact interface; z is the vertical coordinate; qa is the specific discharge of the overlying aquifer flow; x0 is a coordinate extended along the channel centerline; fc is the porosity of the channel material; B is the overlying aquifer thickness; b is the width of the channel; u is the angle of orientation of the channel, with regard to the flowlines of the aquifer flow; j is a local coordinate originating at the upstream (with regard to the overlying aquifer flowline) edge of the mineral penetration site extending along the interface of direct contact, and parallel to the x direction; the coordinate x represents the direction of the aquifer flowlines; y is the horizontal coordinate perpendicular to the flowlines.

The first term of Eq. (1) represents the time variation of the mineral mass incorporated with the control section. The second term represents the net mineral flux leaving the control section by the advection through the channel flow. The right hand side of Eq. (1) refers to the flux of minerals transferred by transverse diffusion –dispersion from the channel control section into the overlying aquifer flow. This term incorporates the net accumulation of mineral mass in the concentration profiles, which originates from the contact between the channel flow and the overlying aquifer flow. We assume that the channel flow rate is not significantly affected by minor changes of the channel orientation along the course of the control section. By applying a dimensional analysis calculation, we find that the first term of Eq. (1) can be neglected, provided that the component of the channel flow velocity in the direction of the overlying aquifer flowline is much larger than the aquifer flow velocity, namely: VR ¼

Va p1 Vc cos u

ð2Þ

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

287

where Va is the flow velocity in the overlying aquifer; Vc is the channel flow velocity, and VR is termed ‘the velocity ratio’. According to Fig. 5, there are the following relationships between the different coordinates: dx0 ¼ dx=cos u;

dy ¼ dx tan u

ð3Þ

Introducing the relationships of Eqs. (2) and (3) into Eq. (1), with some approximation, we obtain:

›Cc q sin u ¼2 a ›x 0 Qc  ðB  Ca dz 0

j¼b=sin u

2

ðB 0

 Ca dz



ð4Þ

j¼0

Besides the graphical description of Fig. 4, we may define curved channels as those channels that transfer minerals only to approaching fresh groundwater of the overlying aquifer—that is, they have only a single site of interaction with any overlying aquifer flowline. Winding channels transfer minerals to approaching groundwater that has been mineralized by a previous interaction with the channel, as well as to fresh groundwater. Winding channels are characterized by multiple mineral penetration sites that supply minerals into some of the overlying aquifer groundwater flowlines. The treatments of the previously mineralized water interactions developed for the winding channel case may be applied to any source of prior mineralization, and is not limited to single-channel applications. Fig. 6 shows some schematics of mineral concentration profile development in the overlying aquifer, due to mineral penetration from the paleodrainage channel. Case (a) refers to a flowline crossing a curved channel or those portions of a winding channel subject to initial contact with approaching fresh groundwater of the aquifer. Cases (b) and (c) refer to winding channels, illustrating possible mineral concentration profiles along flowlines where previously mineralized overlying aquifer groundwater interacts with the channel at a mineral penetration site subsequent to the initial contact. The mineral penetration site representing the channel is basically a line source of variable strength. By integrating the mineral flux over the interface of mineral penetration it is possible to obtain the total amount of minerals supplied from the channel into

Fig. 6. Build-up of the mineral concentration profile over the direct contact interface, which comprises the mineral penetration site. (a) Approaching aquifer groundwater is fresh. (b) and (c) Approaching groundwater is mineralized.

the overlying aquifer, and then to approximate the mineral penetration site as a point source. By using the image method (e.g. Fischer et al., 1979) to simulate the bottom and top of the aquifer, an estimate of mineral advection and diffusion in the aquifer could be obtained. However, this approach requires complete uniformity of the domain. Use of some of the available semi-analytical solutions of heat conduction in solids (Carslaw and Jaeger, 1959) suffers from similar difficulties. Therefore, to address the issues and conditions of the present study, we have developed a BL approach for the analysis and

288

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

calculation of mineral transport in the aquifer overlying paleodrainage channels that act as secondary sources of mineralized water. The following sections present specific derivations and formulations characterizing and quantifying mineral transfer and the build-up of mineral concentration profiles as a result of both curved and winding channel sources.

3. Curved channels In curved channels, the groundwater approaching the mineral penetration site is always freshwater, since at the upstream edge of the channel, which is represented by j ¼ 0; the aquifer flowline does not carry any mineral concentration. The curved channel calculations are also relevant to the portion of the winding channel that is crossed first by freshwater flowlines of the overlying aquifer. In these cases, the second integral on the right-hand side of Eq. (1) vanishes. Along the overlying aquifer flowlines, we identify several sections of development of the mineral concentration profile. Fig. 7 shows sections

of various stages of BL development along the flowline, with the symbols and terminology used in the development that follows (see nomenclature for definition of symbols). 3.1. The mineral penetration section The mineral penetration section comprises the overlying aquifer portion above the mineral penetration site, which is the interface of direct contact between the channel and overlying aquifer flowlines. Here, we apply the conceptual approach of the TSBL method (Rubin and Buddemeier, 1996, 1998a) to analyze the mineralization process. Owing to the mineral transfer from the channel into the overlying aquifer along the interface of direct contact, there is a build-up of a mineral concentration profile similar to a BL in the overlying aquifer, shown in Fig. 6a. The thickness of the mineralized portion of the aquifer given by d0 ; is similar to a BL. The mineral concentration effectively vanishes at the top of the BL. We define a ROI, at the top of which the mineral concentration is smaller than a defined acceptable

Fig. 7. Characteristics of the various sections of the mineral plume described in this study, with terminology and symbols shown.

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

value, CT : The thickness of the ROI is d; which is termed as TSBL. The ROI comprises a portion of the total mineralized zone. According to the BL approximation, the mineral concentration distribution in the overlying aquifer above the mineral penetration site is given by:   z n Ca ¼ Cc 1 2 ð5Þ d0 where d0 is the thickness of the mineralized zone, which is a BL; n is a power coefficient approximately equal to 3 (Rubin and Buddemeier, 1996, 1998a). Due to the small value of the channel width, b; we consider the build-up of the mineralized zone under steady state conditions offlow and mineral concentration distribution. By introducing Eq. (5) into the equation of advection–dispersion (Rubin and Buddemeier, 1998a), we obtain the following expression for the build-up of the BL along the mineral penetration section: ðd20 Þj ¼ 2anðn þ 1Þj

ð6Þ

where a is the transverse dispersivity of the aquifer. By introducing Eq. (5) into Eq. (4) and integrating over a finite interval Dx0 along the channel centerline, we obtain: ( ðDx0 q sin u a Cc ¼ Cc0 exp 2 Qc 0 )  ð1  n 0 ð7Þ £ d0 ð1 2 hÞ dh dx 0

j¼b=sin u

where Cc0 is the mineral concentration at a point of reference; h is the BL coordinate, which is given by: z h¼ ð8Þ d0 Eq. (6) was introduced into Eq. (7) to obtain: 8 sffiffiffiffiffiffiffiffiffiffiffiffiffiffi9 < q Dx0 2anb sin u = a Cc ¼ Cc0 exp 2 : Qc nþ1 ;

ð9Þ

Eq. (9) can be applied to calculate the variation of Cc along a channel segment in which values of b and u are kept constant. It should be noted that Eq. (9) can be applied provided that Eq. (2) is satisfied and the flow in the overlying aquifer is significantly different in direction from the both the channel centerline direction and its perpendicular. Eq. (9) is not

289

applicable to cases like u ¼ 0; p/2. If u < p=2 then such minor quantities of the groundwater flow-rate take place in the direction of the centerline of the narrow channel that the length of the penetration site becomes too small to be an effective source of mineral transfer. If u ¼ 0; then the channel course is parallel to the overlying aquifer flowlines. Under such conditions of parallel flows with different rates in the channel and overlying aquifer, a different type of mineral transport calculation should be applied. Eq. (9) indicates that the process of mineral transfer from the channel into the overlying aquifer occurs along a limited length of the channel and the overlying aquifer as the source strength (concentration) of the channel is progressively decreased by exchange with the overlying aquifer along its length. Therefore, the overlying aquifer flowlines located closer to the source of mineral input to the channel get larger amounts of minerals. We consider that a channel mineral concentration equal to 1% of the value at the channel entrance is equivalent to complete dilution of the mineralized water. At complete dilution, the transfer of minerals from the channel flow into the overlying aquifer effectively stops. If a curved channel had constant values of b and u then the length of complete dilution of that channel, LD would be given by: rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4:6Qc nþ1 LD ¼ ð10Þ 2anb sin u qa In the general case, the curved channel can be treated as consisting of various segments of variable length, width and orientation. We refer to a channel such that in each one of its segments, values of the width and orientation are kept constant. If the number of channel segments in which mineral concentration is greater than the dilution value is J; then J can be determined by applying the following relationships: 2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 J X q L 2anbj sin uj 5 a j 4 $ 4:6 ð11Þ Q nþ1 c j¼1 where j is the number of the particular channel segment; Lj is the length of the segment; uj is the orientation angle of the segment; and bj is the width of the segment. Any cross section of the curved channel, which is located upstream of the complete dilution point,

290

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

represents a mineral penetration site where minerals are transferred into the overlying aquifer through the interface of direct contact between saline and freshwater. Downstream of the mineral penetration site along the aquifer flowline, calculations of the variation of the mineral concentration distribution in the overlying aquifer are provided in following sections. Owing to the small length of the interface of direct contact, the mineralized zone thickness is generally smaller than the aquifer thickness at the downstream edge of the mineral penetration site, namely at j ¼ b=sin u: Downstream of the mineral penetration site, the mineral concentration plume is subject to extension and expansion, until it occupies the entire thickness of the aquifer, at the attachment point (Rubin and Buddemeier, 1998c,d). We refer to the section of the plume between the downgradient end of the penetration site and the attachment point as the ‘extension – expansion section’.

concentration was cr Cc : A particular normalized isohaline cr was chosen to represent the top of the inner BL. Conservation of mass was applied to adjust the mineral concentration profile parameters for the effects of a transition from a profile of constant bottom concentration to a profile of no diffusive flux at the aquifer bottom. Numerous numerical simulations of the BLs developed at the right hand side of the downstream edge of the mineral penetration site were performed (Rubin and Buddemeier, 1998c). These simulations showed that if the thickness of the mineralized zone, d0 is smaller than the thickness B; of the aquifer, and if cr ¼ 0:5; then there is almost no mineral transfer between the inner BL and the outer BL. The mineral concentration distribution in the inner BL is given by:

3.2. Section of mineral concentration plume extension– expansion



A mineralized zone with a thickness smaller than the aquifer thickness characterizes this section of the aquifer. Downstream of the mineral penetration site, the boundary condition at the bottom of the aquifer changes from a constant mineral concentration boundary to a vanishing diffusive mineral flux boundary. There is no more supply of minerals into the aquifer flowlines, but the mineral concentration profile is subject to expansion due to dispersion. The mineralized zone occupies a portion of the aquifer thickness, and the mineral concentration profile consists of two BLs (Rubin and Buddemeier, 1998c), shown in Fig. 7. One BL develops at the bottom of the aquifer and is termed the inner BL. The other one is termed the outer BL, and it develops on top of the inner BL. The thickness of the inner BL is du : In the extension– expansion section, normalized isohaline contours were calculated. Each normalized isohaline was associated with a particular value cr ; which represented the ratio between the mineral concentration of that particular isohaline and the mineral concentration at the bottom of the aquifer. At the left-hand side of the downstream edge of the mineral penetration site, where j ¼ b=sin u; we referred to the ordinate z ¼ du at which the mineral

Ca ¼ Cb ½1 2 ð1 2 cr Þð1 2 lÞn1 ; ð12Þ

du 2 z ; 0 # z # du du

where Cb is the contaminant concentration at the bottom of the overlying aquifer, n1 is a power coefficient; l is the inner BL coordinate. Mineral concentration distribution at the outer BL is represented by: Ca ¼ cr Cb ð1 2 zÞn2 ;



z 2 du ; d0 2 du

ð13Þ

du # z # d0 where n2 is a power coefficient; z is the outer BL coordinate, and d0 is the thickness of the portion of the aquifer subject to mineralization. Previous numerical experiments (Rubin and Buddemeier, 1998c) have shown that the following values are applicable: n1 ¼ 1:84;

n2 ¼ 4

ð14Þ

As long as the top of the aquifer does not bind the mineral concentration plume, the plume incorporates two sections: (a) the extension– expansion section, and (b) the spearhead section, as described by Rubin and Buddemeier (2002a,b). The interface between these two sections represents the greatest vertical thickness of the mineral concentration plume, which migrates with the aquifer flow at an advection velocity

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

equal to the interstitial flow velocity of the overlying aquifer. At the attachment point, this cross section of the mineral concentration plume reaches the top of the aquifer. Then the extension –expansion plume section reaches its final development. The spearhead section proceeds in its migration down gradient with the aquifer flow, and a new section, the restructuring section, starts to develop between the attachment point and the down gradient migrating spearhead section. In the extension –expansion section, the contaminant transport equation yielded values of du and d0 developments as given by the following expressions:

d2u ¼ ðd2u Þx¼0 þ a1

2að1 2 cr Þn1 ðn1 þ 1Þ x; n1 þ c r ð15Þ

0 # x # Va T ðd0 2 du Þ2 ¼ ðd 2 du Þ2x¼0 þ 2a2 an2 ðn2 þ 1Þx;

ð16Þ

0 # x # Va T

where Va is the aquifer interstitial flow velocity, T is the time period required for a fluid particle, emerging from the downstream edge of the mineral penetration site, to arrive at the attachment point, x is the longitudinal coordinate of the extension –expansion section, x ¼ 0 represents the downstream edge of the mineral penetration site, a1 and a2 are constant coefficients. Numerical experiments and the principle of mass conservation (Rubin and Buddemeier, 1998c) were used to determine their values:

a1 ¼ 0:935;

a2 ¼ 0:775

ð17Þ

According to the principle of mass conservation: Cb du ¼ ðCb du Þj¼b=sin u

ð18Þ

Cb d0 ¼ ðCb d0 Þj¼b=sin u

ð19Þ

These expressions provided the value of Cb ; and were also useful for the determination of the coefficients a1 and a2 : Only part of the mineralized zone thickness comprises the ROI, which is the TSBL. In the ROI the mineral concentration is higher than the acceptable value CT : Either one of the following

291

expressions might gives the thickness of the ROI:

d ¼ du ð1 2 lT Þ;   1=n1 C lT ¼ 1 2 2 1 2 T Cb d ¼ zT ðd0 2 du Þ þ du ;

zT ¼ 1 2

ð20Þ 

CT cr C b

1=n2

ð21Þ

The first expression is applicable if d , du : The second expression is applicable if d . du : 3.3. The restructuring section The restructuring section is located downstream of the attachment point and behind the spearhead region, where the thickness of the mineralized zone is the same as the saturated thickness of the overlying aquifer. Two BLs of different thickness characterize the restructuring section (Rubin and Buddemeier, 1998d). The inner BL develops in this section at the expense of the outer BL. The thickness of the inner BL is du : We defined two types of BL coordinates, referring to the inner and outer BL, respectively:



du 2 z ; 0 # z # du du

ð22Þ



z 2 du ; du # z # B B 2 du

ð23Þ

According to Rubin and Buddemeier (1998d), the mineral concentration profiles of the inner and outer BLs were given, respectively by: Ca 2 Ct ¼ 1 2 0:5ð1 2 lÞn3 ; 0 # z # du ð24Þ Cb 2 Ct Ca 2 Ct ¼ 0:5ð1 2 zÞn4 ; du # z # B ð25Þ Cb 2 Ct where Cb and Ct represent the mineral concentration at the bottom and top of the overlying aquifer, respectively. In order to comply with the conservation of mass, some adjustment of BL quantities was made at the attachment point. We defined a longitudinal coordinate X for the restructuring and establishment sections, where X ¼ 0 represented the attachment point. By applying the equation of mineral transfer in the restructuring section, the value of X associated with a given value

292

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

of du in the restructuring section was found to be: "      # an1 X du du0 du 2 du0 2 2 g 2 ¼ 2 b1 2 þ b2 B B B B B  3 2 du 6 b4 2 b5 B 7  7 2 b3 ln6 ð26Þ 4 du0 5 b4 2 b5 B where bi ði ¼ 1; …; 5Þ are coefficients, whose values are given in Appendix A, and g is a coefficient of calibration of the model. Comparison of numerical simulation results with the BL approach implied g < 2:5 (Rubin and Buddemeier, 1998d). At the downstream end of the restructuring section, the establishment point was identified as the location at which the inner BL thickness was equal to approximately half of the thickness of the aquifer. Introducing this value into Eq. (26) yielded the length of the restructuring section: 8     2 > B < 1 du0 LR ¼ 2 b1 2 2 gan3 > B : "   # 1 2 du0 2 2 þ b2 2 B   39 2 1 > > b 2 b 4 5 6 7= 2 6 7  5 þb3 ln4 ð27Þ > d > ; b4 2 b5 u0 B By the employment of mass conservation and previous expressions, the mineral concentrations at the bottom and top of the aquifer in the restructuring section were determined to be given, respectively, by: Cb ¼

Cb0 2ðn3 þ 0:5Þðn4 þ 0:5Þ 2 0:5   d B 2 du0  2ðn3 þ 0:5Þðn4 þ 0:5Þ u0 2 0:5 ð28Þ du B 2 du

  Cb0 ðn3 þ 0:5Þ B 2 du0 du0 2 Ct ¼ ð29Þ 2ðn3 þ 0:5Þðn4 þ 0:5Þ2 0:5 B 2 du du where Cb0 is the mineral concentration of the aquifer bottom at the right hand side of the attachment point. Results of numerical experiments (Rubin and Buddemeier, 1998d) suggested the following values

of the power coefficients: n3 ¼ 1:5;

n4 ¼ 2:5

ð30Þ

3.4. The establishment section In the establishment section, as in the restructuring section, there are two BLs, inner and outer. Unlike the case of the restructuring section the dimensions of the two BLs in the restructuring section do not change, but their concentrations do change. The inner BL occupies the lower half of the aquifer thickness, and the upper BL occupies the upper half. Mineral concentration profiles in the inner and outer BLs are given by Eqs. (24) and (25), respectively, while considering that du ¼ B=2: However, in the establishment section, both BLs have the identical power coefficient n5 : Therefore some adjustment was made to the values of Cb and Ct ; to comply with the principle of mass conservation. At the downstream end of the establishment section, the concentrations in the two BLs have effectively equilibrated, and the difference in mineral concentration between the top and bottom of the aquifer is very small. For an appropriate definition of ‘very small’, we assumed that the mineral concentration profile was uniform if either of two conditions was met: if mineral concentrations at all points in the cross section were smaller than 1% of the mineral concentration entering the paleochannel (entrance mineral concentration), or, if the difference between Cb and Ct was smaller than 1% of the entrance mineral concentration. The value of X; which represents the downstream end of the establishment section is termed the location of the ‘uniformity point’. The length LE of the establishment section was determined by the criteria of the preceding paragraph as: LE ¼

B2 ln½100ðCb 2 Ct Þes  8a3 an5

ð31Þ

where a3 is a constant coefficient and the subscript ‘es ’refers to the right hand side of the establishment point (see Fig. 7). The symbol n5 represents the power coefficient of the establishment section; its value is about n5 ¼ 1:8: The value of the coefficient a3 ; according to numerical experiments (Rubin and Buddemeier, 1998d), is about 0.6.

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

Values of Cb and Ct at the right hand side of the establishment point were obtained by applying the principle of mass conservation to adjust the BL quantities. The mineral concentrations at the bottom and top of the aquifer in the establishment section, obtained by using the equation of mineral concentration transport, are given, respectively, by: Cb ¼

1 2

½ðCb þ Ct Þes 2 ðCb 2 Ct Þes F

ð32Þ

Ct ¼

1 2

½ðCb þ Ct Þes 2 ðCb 2 Ct Þes F

ð33Þ

where   8an F ¼ exp 2a3 2 5 ðX 2 LR Þ : B

ð34Þ

4. Winding channel The equations and expressions developed for the curved channel case are also applicable to those portions of the winding channel that are in contact with approaching overlying aquifer groundwater free of minerals. Differences between calculations referring to curved and winding channels occur in locations where minerals penetrate from the winding channel into overlying aquifer groundwater that has been previously mineralized by contact with an upgradient section of the winding channel. Mineral transfer from the mineralized overlying aquifer into the winding channel is also possible, provided that the direct contact interface is located downstream of the complete dilution point of the channel. However, because the mineral concentration values of the approaching overlying aquifer groundwater were usually small, we ignore such phenomena, and consider only the portion of the channel located upstream of the point of complete dilution. Calculations for the curved channel case, which refer to sections of the overlying aquifer downstream of the initial mineral penetration site, are basically applicable to winding channel calculations, provided that adjustments of mineral concentration profile quantities are made at points of interfaces between adjacent sections of the mineral concentration plume. For the mineral penetration site of the winding channel, the calculation of Cc can be performed according to Eq. (4), without the simplification of Eqs. (7) and (9). Mineral penetration sites may occur in any

293

of the sections (extension – expansion, restructuring, establishment, or spearhead), which characterize the mineralized portion of the overlying aquifer. An approaching mineralized groundwater front has the shape of a spearhead whose length is equal to b=sin u; where b and u refer to the last mineral penetration site of the winding channel. However, since we are considering relatively narrow channel widths, the spearhead section is comparatively small and we neglect phenomena associated with this section. If the approaching mineral concentration plume occupies only a portion of the aquifer thickness then the upstream edge of the new mineral penetration site is characterized by a profile of the extension– expansion type. Other options for the mineral concentration profile at that edge are profiles characterizing either the restructuring or establishment section. The different mineral concentration profiles that may take place at the upstream edge of the mineral penetration site are shown in Figs. 6 –8. In all cases of mineralized groundwater approaching a new mineral penetration site, the approaching mineral concentration profile incorporates two BLs. As groundwater of the overlying aquifer moves over the new mineral penetration site, the mineral concentration profile is subject to an abrupt change in its bottom boundary condition, which leads to restructuring of the mineral concentration profile as a single BL. Calculations for this restructuring process adjust the BL quantities based on the principle of mass conservation. If the adaptation of the mineral concentration profile at the upstream edge of the new mineral penetration site ðj ¼ 0Þ creates a mineralized zone with a thickness smaller than the thickness of the aquifer then the profile above the mineral penetration site, is represented by Eq. (5). Then the equation of mineral transport can be applied to calculate the thickness of the mineralized zone at the downstream edge of the mineral penetration site, where j ¼ b= sin u : ðd20 Þj¼b=sin u ¼ ðd20 Þj¼0 þ 2anðn þ 1Þb=sin u

ð35Þ

By introducing the expressions for the mineral concentration profiles in Eq. (4), the following expression was obtained: dCc q sin u ½ðd Þ ¼2 a 2 ðd0 Þj¼0 dx0 Cc Qc ðn þ 1Þ 0 j¼b=sin u

ð36Þ

294

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

Fig. 8. Restructuring of the mineral concentration profile at the upstream edge of the mineral penetration site ðj ¼ 0Þ of a winding channel. (a) Thickness of mineral plume is smaller than the aquifer thickness, B: (b) Thickness of the mineral plume occupies the entire thickness of the aquifer.

Direct integration of Eq. (36) yielded:   C q sin u ln c ¼2 a Cc0 Qc ðnþ1Þ ðDx0 ½ðd0 Þj¼b=sin u 2ðd0 Þj¼0 dx0 £

the aquifer is smaller than in the curved channel case with identical values of b and u because the average concentration difference between channel and aquifer water is less. Therefore, the calculated complete dilution length LD of a winding channel is longer than that of a comparable curved channel. At the downstream edge of the mineral penetration site, the mineral concentration profile should accommodate to the abrupt change in the bottom boundary condition. Therefore, we implemented conservation of mass to adjust the BL quantities. Fig. 8 shows mineralized zone of the approaching mineral concentration plume that occupies the entire thickness of the overlying aquifer. In this case the mineral concentration profile of the approaching plume consists of two BLs. Therefore, the mineral concentration profile could be expressed by Eqs. (24) and (25), provided that the approaching mineral concentration plume profile was typical of the restructuring section (see Fig. 7). Eqs. (24) and (25) were also applied with n5 replacing n3 and n4 if the approaching mineral concentration profile was of the establishment section. As in the previous case, the abrupt change in the bottom boundary condition of the mineral concentration profile at the upstream edge of the mineral penetration site suggests adapting the profile of two BLs to a mineral concentration profile of a single BL. The mineral concentration profile of this single BL occupying the entire thickness of the aquifer was approximated by: Ca ¼ Ct þ ðCc 2 Ct Þð1 2 vÞn6

ð39Þ

where n6 is a power coefficient, and v is the BL coordinate, which is given by: ð37Þ

0

Where there was a minor difference between the thickness of the mineralized zone at the downstream and upstream edges of the mineral penetration site Eq. (37) was approximated by:   0 Cc q anb ðDx dx0 ln ð38Þ ¼2 a Qc 0 ðd0 Þj¼0 Cc0 It should be noted that in the winding channel case, the average intensity of mineral penetration into



z B

ð40Þ

The value of the power coefficient n6 of Eq. (39) was determined by applying the analytical solution for an analogous problem of heat conduction in solids. The analogy originates from the identity between the equation of mineral transport with its boundary conditions at the mineral penetration site to the equation of heat transfer and build-up of the temperature profile in a solid slab subject to initial temperature profile and bounded by two parallel planes, where one of them is insulated and the other

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

one is subject to constant temperature. This analogy implies: 1 X

"

Ca 4 1 ð2mþ1Þ2 p2 aj exp 2 ¼ 12 p m¼0 2mþ1 Cc 4B2   pz £sin ð2mþ1Þ 2B

#

ð41Þ

where j is the local longitudinal coordinate of the mineral penetration site, and the location of j ¼ 0 is chosen by curve fitting to the mineral concentration profile at the upstream edge of the mineral penetration site. The BL approximate mineral concentration profile represented by Eq. (39) was introduced into the equation of mineral transport. Then by using Leibniz’s theorem and performing integration of the expressions we obtained the following expression for the calculation of Ct :  ln

Cc 2 Ct Cc



¼2

aðn6 þ 1Þj B2

ð42Þ

In order to find the appropriate value of n6 ; calculations were performed with regard to the development of mineral concentration profiles for various values of j by applying Eqs. (39) and (42), as well as by using Eq. (41). It was found that about six terms were needed to obtain practically complete convergence of the series expansion shown in Eq. (41). The numerical experiments indicated that for small values of j the series expansion solution of Eq. (41) was not accurate, and could be replaced by the complementary error function solution, which is typical of semi-infinite domain. Under such conditions the value of n6 could be replaced by n; namely 3. For large values of j; smaller values of n6 were obtained. However, generally values between 1.2 and 1.8 were found to be satisfactory. Only minor errors were introduced into the calculation by considering a constant value of n6 for the entire domain. At the downstream edge of the mineral penetration site, due to the change of the bottom boundary condition of the mineral plume, the mass conservation principle again required some adjustment of the mineral concentration profile.

295

5. Practical use of the method As noted by one of the reviewers of this manuscript, the approach and method described in this paper can be useful for the quick analysis, evaluation, and possibly prediction, of various cases of contaminant supply and migration in aquifers containing minor portions of strata with very high permeability and major portions of strata with comparatively low permeability. Only very limited quantities of field data were available for such applications in the framework of this study. However, as presented and discussed in following paragraphs, various studies have already applied the approach of this paper to practical issues of sensitivity analysis, testing of numerical modeling results, and evaluation of mineralization processes at some specific field sites. All these cases indicate that in complex strata consisting of high and low permeability portions, very often aquifer contamination results mainly from high rates of contaminant advection in the high permeability portions, contaminant transfer from those portions into the low permeability portions by diffusion – dispersion through contact interfaces between the high and low permeability portions, and contaminant advection in the low permeability portions. The method developed in this study was used to analyze the mineralization induced in the Great Band Prairie aquifer of south central Kansas by a salinewater channel, and the potential effect of that mineralization on both groundwater and discharged surface water quality. Fig. 9 shows the particular paleochannel to which the calculations were applied as a test case; the location of this channel is indicated in Fig. 1. Water level measurements in the overlying aquifer during 1996 provided contours of the piezometric heads, shown with 20 ft intervals in Fig. 9. We used the piezometric head contours to depict flowlines of the aquifer flow in the region subject to mineralization by the saline water flowing through the channel. Calculations of the build-up and spatial variation of the mineral concentration profile along several of the flowlines were performed. The calculations provided details about variations of mineral concentration along the course of the channel flow up to the point of complete dilution or uniformity of the mineral concentration profile. Details obtained for the specific flowlines identified by A, B, and C in Fig. 9

296

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

Fig. 9. A specific channel (see Fig. 1 for location) with an assumed mineral source and flowlines of the overlying aquifer indicated. Distances in km are indicated along the channel and the flowline, marked by A and used for sample calculations. Points A, B, and C represent the contact locations for the sample calculation flow line (on which Cav ¼ 0:025Cc0 ) and the flowlines with Cav ¼ 0:01Cc0 and Cav ¼ 0:005Cc0 ; respectively.

are shown in Figs. 1 –13. The longitudinal coordinate x is emanating at the downstream edge of the bedrock mineral penetration site (Fig. 1), and is in an area known to have saline water in the aquifer just above the bedrock surface (Young et al., 2000). Point A is located where the paleodrainage channel intersects the region of saline bedrock exposure. Along the channel, as along the flowline, numbers show segments of 1 km length with an origin at point A. The average hydraulic gradient along the overlying aquifer flowlines is about 0.33%. For the first portion (about 12 km) of the channel, the hydraulic gradient along the channel course is about 0.2%. As there are no direct measurements of channel dimensions at this location, we assumed that the average width and depth of the channel were 30 and 5 m, respectively (see

the schematic of Fig. 3). Therefore the assumed cross section area of the channel was about 150 m2. It was also assumed that the channel width at the interface of direct contact ðb sin uÞ was about b ¼ 50 m. The channel hydraulic conductivity was taken as Kc ¼ 500 m/d, and the hydraulic conductivity of the aquifer was assumed to be about Ka ¼ 10 m/d. These values are consistent with the upper and lower ranges, respectively, of the values reported by Layton and Berry (1973), Fader and Stullken (1978), and Cederstrand and Becker (1998). Therefore, the calculated channel flow-rate was about 3 Qc ¼ 150 m /d, and the specific discharge of the aquifer flow was qa ¼ 3:3 £ 1022 m/d. Along the first 12 km of the channel, its assumed angle of orientation to the aquifer flowlines was u ¼ 508: Considering an

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

297

Fig. 10. Variations of the mineralization zone characteristics at the extension–expansion section along the sample calculation flowline, shown in Fig. 8, and marked by A. (a) Growth of d0 and du until d0 ¼ B: (b) Variation of Cb :

aquifer dispersivity of a ¼ 0:1 m, and introducing values of the appropriate parameters into Eq. (10), the obtained dilution length of the channel was about LD ¼ 8700 m. However, the study of Garneau (1995) suggested considerably smaller values of the aquifer dispersivity. If a ¼ 0:052 m, then the obtained dilution length was about 12 km. For smaller values of the dispersivity, variations of the channel

orientation were taken into account. If a ¼ 0:01 m was assumed, then the channel was defined as a winding channel. For such a small dispersivity value, the channel segments between points at 17 and 24 km interact with aquifer groundwater that was mineralized by upgradient segments of the channel. It was assumed that the South Fork Ninnescah River gets some saline water from the mineralized

298

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

Fig. 11. Variations of the mineral plume characteristics at the restructuring section along the sample calculation flowline shown in Fig. 8, and marked by A. (a) Variation of Cb and Ct : (b) Development of du :

portion of the aquifer. Therefore, we calculated mineral concentration profile variations in the flowline shown in Fig. 9 up to a distance of 13 km downstream from the mineral source, where the aquifer flowline reaches the river. The calculation results shown were obtained by using a ¼ 0:01 m; however, they can easily be modified with other values for dispersivity to study the sensitivity of the physical phenomenon to this parameter. The method of this paper could also be applied to calculations of

mineral transport downstream of the river, provided that some more information about mineral discharge into the river was available. By applying Eq. (6), we obtained d0 ¼ 3:95 m, du ¼ 0:82 m as values for the mineralized zone and inner BL thickness, respectively, at the downstream edge of the mineral penetration site ðj ¼ b=sin uÞ: By applying Eqs. (15) and (16), we found that the attachment point (see Fig. 7) was located at x ¼ 3200 m, were x is measured along flowline

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

299

Fig. 12. Variation of Cb and Ct in the establishment section of the mineral plume along the sample calculation flowline shown in Fig. 8, and marked by A.

A. The mineral concentration at the bottom of the overlying aquifer at that location was about 10% of the channel mineral concentration at its entrance. At the attachment point, the thickness of the inner BL was modified to comply with the mass conservation principle, changing it from du ¼ 8:2 to 11.4 m. By applying either Eq. (20) or Eq. (21), we obtained

the value of the ROI thickness. By considering that CT ¼ 0:01Cc0 ; where Cc0 is the mineral concentration at the entrance of the channel, we obtained d ¼ 32:1 m. Considering that the aquifer thickness was B ¼ 40 m, that value of d indicated that although the entire thickness of the aquifer was mineralized at the attachment point, the top 8 m of the aquifer

Fig. 13. Development of inner and outer BLs in the entire simulated domain and the ROIs ðdÞ for the flowlines B ðCav ¼ 0:01Þ and C ðCav ¼ 0:005Þ for the system shown in Fig. 8.

300

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

contained groundwater that was practically fresh. By applying Eqs. (28) and (29), we calculated values of the mineral concentration at the top and bottom of the aquifer at the establishment point. We found that Cb was about 4% of the channel entrance concentration, and Ct was about 1.1% of the channel entrance concentration. If an acceptable mineral concentration is 1% of the channel entrance concentration, our result indicated that the aquifer was completely mineralized at the establishment point. By applying Eq. (27) we obtained a restructuring length of about LR ¼ 6000 m. This result indicated that where the sample flowline intersected the South Fork Ninnescah River its mineral concentration profile was in the establishment section (see Fig. 7). However, only about 4000 m of that region could practically be considered, as the distance between the channel entrance and the river measured along the aquifer flowline, was about 13.2 km, and the attachment point was at 3200 m downstream of the channel. Nonetheless, we applied Eqs. (31) – (34) and calculated variations of the mineral concentration profile until achieving the hypothetical condition of uniformity, in which difference of mineral concentration between bottom and top of the aquifer was smaller than 0:01Cc0 : By applying Eq. (31), and ignoring any effect of the South Fork Ninnescah River, we obtained the length of the establishment section as LE ¼ 20 km. Eqs. (34) and (35) indicate that at the downstream end of the establishment section, the mineral concentration distribution was practically uniform, and its value was about 2.4% of the mineral concentration of the channel flow at its entrance. This was, of course, the value of the average mineral concentration of the flowline. However the establishment section calculations were relevant only up to the point of the South Fork Ninnescah River, at which values of the mineral concentration at the bottom and top of the aquifer were about 3.7 and 1.4%, respectively, of the mineral concentration at the entrance of the channel. At that point, substantial quantities of the groundwater probably discharge into the South Fork Ninnescah River. Simulation of mineralization characteristics of the flowline considered in Fig. 9 could easily be adapted for neighboring flowlines. In Figs. 10 – 13,

details of such simulations and their implications are provided. Fig. 10a describes the growth of the inner BL and the mineralized zone at the extension– expansion section. The coordinate x; refers to distance from the downstream edge of the mineral penetration site measured along the flowline. Fig. 10b shows the variation of the mineral concentration at the bottom of the overlying aquifer. The value of Cb =Cc0 in this figure represents the mineral concentration at the bottom of the aquifer normalized with regard to Cc0 : Fig. 11a refers to the restructuring section. It shows the variation of the normalized mineral concentration at the top and bottom of the aquifer. Fig. 11b also refers to the restructuring section. It shows the moderate rate of growth of the inner BL, until its thickness occupies the lower half of the overlying aquifer. Fig. 12 refers to the establishment section. It indicates that the normalized value of Cb decreases and that of Ct increases, until they become almost identical. However, owing to the very small vertical mineral concentration gradients in the overlying aquifer cross section, complete uniformity of the mineral concentration profile may be reached only at a very large distance downstream of the mineral penetration site. Calculations of the thickness of the various BLs, as well as calculations of values of Cb and Ct normalized with regard to Cc ; are not affected by the value of Cc : However, the thickness d of the ROI is determined according to the definition of an acceptable value of mineral concentration, CT : We assumed that CT ¼ 0:01Cc0 : We identified the overlying aquifer flowlines that carried mineral concentration larger than CT indefinitely, and those carrying it only a limited distance along their course. The selection of the different types of flowlines depends on the average mineral concentration of the flowline. This quantity is obtained by dividing the total mineral quantity present at the overlying aquifer cross section per unit width of the aquifer, by the total overlying aquifer thickness. This quantity is kept constant along the overlying aquifer flowline. We integrated the mineral concentration profile over the overlying aquifer thickness at the downstream edge of the mineral penetration site to obtain the following expression for the average

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

mineral concentration of the flowline: Cav ¼

1 d0 Cc B nþ1

ð43Þ

By applying Eqs. (43), (6) and (9), we identified the relationship between the location, Lc of the channel cross section (downstream of the channel entrance) and its average value of mineral concentration Cav : rffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffi! Qc nþ1 Cav B nþ1 Lc ¼ 2 ln ð44Þ Cc0 2anb=sin u qa 2anb sin u Flowlines with an average mineral concentration larger than CT convey water defined as indefinitely mineralized or saline indefinitely. A flowline with average mineral concentration smaller than CT ; conveys mineralized water along a portion of its course; although the mineral load remains constant, dispersion ultimately results in a concentration profile that does not exceed the defined level of mineral concentration at any subsequent point along the flowline. The average mineral concentration of the flowline marked by A in Fig. 9 was found to be about 0:025Cc0 ; so it conveyed saline water indefinitely. It should be noted that along the course of the flowline marked by A, there was no input or output of minerals, so the average mineral concentration along that aquifer flowline was constant. We considered two other flowlines of Fig. 9: the flowline with average mineral concentration Cav ¼ 0:01Cc0 ; and the flowline with average mineral concentration Cav ¼ 0:005Cc0 : The first flowline, marked by B in Fig. 9, crossed the channel at a distance of 5500 m downstream from the channel entrance. The other flowline, marked by C, crossed the channel at a distance of 9600 m downstream. In Fig. 13, the variation of the ROI thickness ðdÞ along these two flowlines is depicted. It should be noted that the value of d; associated with Cav ¼ 0:01Cc0 ; may represent the practical thickness of the mineral plume, which is bounded by the isohaline at which Ca ¼ Cav for all flowlines in proximity to the channel entrance. As Fig. 13 covers a significant portion of the course of the aquifer flowlines, the adaptation of values of du ; and thereby also d; at the attachment point is clearly shown in this figure. Fig. 13 shows that along the first 6.3 km of the flowline of Cav ¼ 0:005Cc0 saline water saturated a layer whose maximum thickness was

301

about 8.6 m. The flowline of Cav ¼ 0:01Cc0 represents the boundary between flowlines that were indefinitely saline, and those that became practically fresh in a finite distance downstream from the channel.

6. Discussion Figs. 9 –13 and the calculations and results of Section 5 demonstrate that the major parameters which determine the characteristics of the channel mineralization process are: (a) the ratio between the specific discharge of the aquifer and the channel flowrate, (b) the width of the channel, (c) the orientation of the channel, and (d) the transverse dispersivity of the aquifer. The effect of each of these parameters is quantified, and the sensitivity of the mineralization process to changes in each parameter is easily inferred from the relevant analytical expressions. In following paragraphs we use the concept of a stream tube to represent the three-dimensional volume of the aquifer subject to flow. This volume is comprised of, and bounded by, streamlines. The stream tube can be represented by a three-dimensional section of the aquifer; we consider the flowline and the mineral concentration profile as two-dimensional representatives of the stream tube and its transverse average mineral concentration profile. The width of the channel and its orientation determine the length of the mineral penetration site, which is the interface of direct contact along which salt is transferred from the channel into the overlying aquifer stream tubes. The dispersivity of the aquifer is a measure of the capability of the aquifer to obtain minerals from a unit length of the interface of direct contact. For small values of the channel orientation angle the interface of direct contact is long, and larger portions of the channel are subject to contact with mineralized groundwater. Therefore, as indicated by Eq. (10), small orientation angles lead to longer dilution length of the channel. In cases of large channel orientation angle, approaching 908, there is probably a relatively small hydraulic gradient along the channel, so that the channel flow rate is diminished, and saline water may not be transported through the channel over long distances. A large channel width increases the size of the interface of mineral penetration site and intensifies

302

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

the mineral transfer from the channel to the aquifer, so the dilution length of the channel is reduced. An increase of aquifer dispersivity intensifies the mineral transfer from the channel into the overlying aquifer, reducing the dilution length of the channel. For a channel of constant width and orientation, the channel dilution length is inversely proportional to the value of the dispersivity. Our calculations indicated that for the case-study channel (Fig. 9) and a ¼ 0:1 m, the dilution length was about 8700 m. In cases of high dispersivity the number of stream tubes was comparatively small, as was the total number of overlying aquifer stream tubes subject to mineralization. However, the majority of the mineralized stream tubes contained significant amounts of minerals, and in such stream tubes the mineral concentration was high along the entire course of the aquifer flowlines. On the other hand, for a ¼ 0:01 m, the channel became a winding channel, and its dilution length became larger than 30 km. Therefore, in Fig. 9, a comparatively large portion of the aquifer was potentially subject to mineralization if the dispersivity was small. However, in such a case, comparatively minute quantities of salt were transferred to each stream tube of the overlying aquifer, and the mineralized stream tubes had low average mineral concentrations upstream of the point of complete dilution of the channel flow. Such stream tubes conveyed water meeting the ROI definition of mineral concentration along comparatively short distances from the channel. Parameters of the channel and the aquifer, mentioned in the preceding paragraphs, determine how the flux of minerals is transferred to the overlying aquifer stream tubes, and how the mineral concentration is finally distributed in the aquifer. For a stream tube of specified width, one of two conditions occurs: (a) significant mineralization of a small number of stream tubes, or (b) minor mineralization of many stream tubes. The first option creates mineralized stream tubes, the majority of which carry minerals along their entire course. The latter option creates a larger number of mineralized stream tubes that carry minerals to a limited length of their course. Our calculations were based on the assumption that there were only minor changes in the streamline pattern. However, changes in head and flow pattern may occur due to seasonal variations originating from

natural causes, or because of artificial pumping. If such changes are significant they may affect the distribution and direction of advection of the minerals in the aquifer, and should be taken into account.

7. Summary Groundwater mineralization may occur as a result of the flow of saline water in bedrock channels with high permeability. Brines originating from deep saline bedrock formations mineralize these channel formations, which serve as conduits for rapid transmission of brine to parts of the aquifer distant from the initial mineral intrusion site. The mineralization phenomenon was analyzed by applying the conceptual approach of the TSBL, in which it is assumed that the mineral concentration profile developed in the aquifer may consist of one or two BLs. Minerals are transferred from the channel into the overlying aquifer through interfaces of direct contact between the channel flow and the overlying aquifer flowlines. A single BL can represent the mineral concentration profile developed on top of the interface of direct contact, which represents a mineral penetration site. Downstream of the edge of that site, treating the mineral concentration profile as consisting of two BLs, an inner and an outer BL, permits useful approximations and classifications. The profile of the mineral plume is subject to changes originating from the mineral dispersion and changes of the boundary conditions imposed by the top and bottom of the aquifer. Several different sections of the aquifer flowlines are identified, to characterize variations in the mineral concentration profile for particular regions of the aquifer. We do not consider the leading edge of the plume (call the spearhead section) because of its limited extent. In the first major section of plume (the extension– expansion section) the thickness of the mineral plume is smaller than the thickness of the aquifer. The following section is termed the restructuring section. Here the mineral plume occupies the entire thickness of the aquifer, and the section is characterized by an expanding inner BL that is less than half of the aquifer thickness. In the final section, termed the establishment section, each of the BLs employed occupies half of the aquifer thickness, but there is exchange of

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

mineral mass between them so that their mineral concentration values gradually become identical. Characteristic expressions were developed for the calculation of mineral transport and development of the BLs in the different sections of the mineral plume. We identify two general cases of paleodrainage channels: (a) curved channels, and (b) winding channels. Curved channels transfer minerals only to approaching freshwater of the overlying aquifer. Winding channels transfer minerals to both fresh and mineralized approaching groundwater. Appropriate expressions were developed for the quantification of the mineral transfer from each type of channel to the aquifer flowlines. The methods of analysis and calculation were exemplified by application to a particular scenario of channel mineralization in south-central Kansas. The quantitative results of this example were used to identify major parameters of the mineralization process. It was shown that the channel might supply significant quantities of minerals to a small number of stream tubes of the aquifer, or relatively inconsequential quantities of minerals to a comparatively large number of stream tubes, depending on the dimensions and hydrologic characteristics of

303

the aquifer-channel system. The outcomes of these two options, and their implications for water quality, were discussed.

Acknowledgements This work was supported under contracts between the Kansas Water Office and the Kansas Geological Survey. The authors thank Marios Sophocleous and David Young for valuable comments and information. The authors are grateful for the assistance of Mark Schoneweis in preparing the illustrations.

Appendix A. Values of coefficients of Eq. (26) By the employment of BL approach and some assumptions with regards to mechanisms leading to the development of the inner BL in the restructuring section Rubin and Buddemeier (1998d) showed that the values of the coefficients bi ði ¼ 1; …; 5Þ of Eq. (26) are given by:

   d du0 2ðn1 þ 0:5Þ 1 2 u0 ½n2 ðn1 þ 1Þ þ n1 ðn2 þ 0:5Þ B B b1 ¼      du0 du0 2 2ðn1 þ 0:5Þðn2 þ 0:5Þ þ ðn1 þ 1Þ 1 2 B B

ðA1Þ



   du0 du0 2n2 ðn1 þ 0:5Þ 2 n1 1 2 B B   b2 ¼  du0 d 2 2ðn1 þ 0:5Þðn2 þ 0:5Þ þ ðn1 þ 1Þ 1 2 u0 B B

ðA2Þ

   d du0 2 4ðn1 þ 0:5Þ2 ðn2 þ 0:5Þ 1 2 u0 ½n2 ðn1 þ 1Þ þ n1 ðn2 þ 0:5Þ B B b3 ¼     3 d d 2ðn1 þ 0:5Þðn2 þ 0:5Þ u0 þ ðn1 þ 1Þ 1 2 u0 B B

ðA3Þ

304

H. Rubin, R.W. Buddemeier / Journal of Hydrology 277 (2003) 280–304

  d b4 ¼ 2ðn1 þ 0:5Þðn2 þ 0:5Þ u0 B   d b5 ¼2ðn1 þ 0:5Þðn2 þ 0:5Þ u0 B   du0 þ ðn1 þ 1Þ 1 2 B

ðA4Þ

ðA5Þ

In Eqs. (A1) –(A5), du0 is the thickness of the inner BL at the right hand side of the attachment point.

References Buddemeier, R.W., Sophocleous, M.A., Whittemore, D.O., 1992. Mineral intrusion—investigation of salt contamination of groundwater in the Eastern Great Bend Prairie Aquifer, OpenFile Report 92-25, Geohydrology, Kansas Geological Survey, The University of Kansas, Lawrence, KS. Buddemeier, R.W., Garneau, G.W., Young, D.P., Whittemore, D.O., Zehr, D., Lanterman, J., Ma, T.S., Falk, S., 1994. The mineral intrusion project: progress and activities during fiscal year 1994, Open-File Report 94-28, Geohydrology, Kansas Geological Survey, The University of Kansas, Lawrence, KS. Carslaw, H.S., Jaeger, J.C., 1959. Conduction of Heat in Solids, Second ed, Clarendon Press, Oxford, p. 100. Cederstrand, J.R., Becker, M.F., 1998. Digital map of hydraulic conductivity for the High Plains aquifer in parts of Colorado, Kansas, Nebraska, New Mexico, Oklahoma, South Dakota, Texas, and Wyoming, US, Geological Open File Report 98548. Cobb, P.M., 1980. The distribution and mechanisms of salt water intrusion in the fresh water aquifer and in Rattlesnake Creek Stafford County Kansas. MS thesis, Department of Civil Engineering, The University of Kansas, Lawrence, KS. Fader, S.W., Stullken, L.E., 1978. Geohydrology of the Great Bend Prairie South-Central Kansas, Irrigation Series 4, Kansas Geological Survey. Fischer, H.B., List, E., Koh, R.C.Y., Imberger, J., Brooks, N.H., 1979. Mixing in Inland and Coastal Waters, Academic Press, New York. Garneau, G.W., 1995. Detection and characterization of the distribution of mineral intrusion in the Great Bend Prairie aquifer—south-central Kansas, Open-file Report 95-35, Kansas Geological Survey, The University of Kansas, Lawrence, KS. Latta, B., 1950. Geology and groundwater resources of Barton and Stafford Counties, Kansas Bulletin 88, Kansas Geological Survey. Layton, D.W., Berry, D.W., 1973. Geology and ground-water resources of Pratt County, South-Central Kansas, Bulletin 205, Kansas Geological Survey. Rubin, H., Buddemeier, R.W., 1996. A top specified boundary layer (TSBL) approximation approach for the simulation of ground-

water contamination processes. Journal of Contaminant Hydrology 22, 123 –144. Rubin, H., Buddemeier, R.W., 1998a. Application of the top specified boundary layer (TSBL) approximation to initial characterization of an inland aquifer mineralization. Part 1: direct contact between fresh and saltwater. Journal of Contaminant Hydrology 32, 353–376. Rubin, H., Buddemeier, R.W., 1998b. Application of the top specified boundary layer (TSBL) approximation to initial characterization of an inland aquifer mineralization. Part 2: seepage of saltwater through semi-confining layers. Journal of Contaminant Hydrology 32, 377–402. Rubin, H., Buddemeier, R.W., 1998c. Approximate analysis of groundwater mineralization due to local discontinuities in impermeable layer—Part 1. Direct contact between fresh and saltwater, Open-File Report 98-31, Kansas Geological Survey, The University of Kansas, Lawrence, KS. Rubin, H., Buddemeier, R.W., 1998d. Approximate analysis of aquifer mineralization due to horizontal penetration of salinity, Open-File Report 98-33, Kansas Geological Survey, The University of Kansas, Lawrence, KS. Rubin, H., Buddemeier, R.W., 2002a. Groundwater contamination downstream of a contaminant penetration site. Part 1: extension—expansion of the contaminant plume. Journal of Environmental Science and Health Part A A37 (10), 1781–1812. Rubin, H., Buddemeier, R.W., 2002b. Groundwater contamination downstream of a contaminant penetration site. Part 2: horizontal penetration of the contaminant plume. Journal of Environmental Science and Health Part A A37 (10), 1813–1839. Sophocleous, M.A., 1991. Stream-floodwater propagation through the Great Bend alluvial aquifer Kansas—field measurements and numerical simulations. Journal of Hydrology 124, 207–228. Whittemore, D.O., 1993. Ground water geochemistry in the mineral intrusion area of groundwater management district no. 5, SouthCentral Kansas, Kansas Geological Survey, The University of Kansas, Lawrence, KS. Young, D.P., 1992. Mineral intrusion: geohydrology of Permian bedrock underlying the Great Bend Prairie aquifer in southcentral Kansas, Open-File Report 92-44, Kansas Geological Survey, The University of Kansas, Lawrence, KS. Young, D., Rubin, H., 1998. Determining vertical flow in variable density groundwater—a hydrostatic approach applied in southcentral Kansas. Eighteenth Annual Hydrology Days, Colorado State University, Fort Collins, Colorado March–April, 1998. Young, D.P., Buddemeier, R.W., Whittemore, D.O., 1998. Equus Beds Mineral Intrusion Project Report, FY 1998, Open-File Report 98-24, Kansas Geological Survey, University of Kansas, Lawrence, KS. Young, D.P., Buddemeier, R.W., Whittemore, D.O., Rubin, H., 2000. Final Summary and Data Report: The Equus Beds Mineral Intrusion Project, Open-File Report 2000-30, Kansas Geological Survey, University of Kansas, Lawrence, KS, http:// www.kgs.ukans.edu/Hydro/Equus/OFR00_30/index.htm.